數(shù)值分析課程課程設(shè)計(jì)

1、課程設(shè)計(jì)設(shè)計(jì)題目數(shù)值分析 學(xué)生姓名李飛吾學(xué)號 xxxxxxxx專業(yè)班級 信息計(jì)xxxxx班 指導(dǎo)教師設(shè) 計(jì)題 目共15題如下成績課 程 設(shè) 計(jì) 主 要 內(nèi) 容設(shè)計(jì)目的：通過不同題目的理解，進(jìn)行算法分析。通過MATLAB軟件進(jìn)行編程對題目進(jìn)行解決。個別題目設(shè)計(jì)驗(yàn)證，加深對數(shù)值分析的理解。函數(shù)的圖像繪制的運(yùn)用設(shè)計(jì)題目：題11 利用逆向遞推的方法求解問題，通過條件終止地推題12 從某個初始值開始，利用遞推公式進(jìn)行積分估值題13 繪制Koch分形曲線，節(jié)點(diǎn)之間的關(guān)系與坐標(biāo)變換題21 用高斯消元法的消元過程作矩陣分解，LU分解題22 矩陣分解方法求上題中A的逆矩陣，針對不同的b,而重復(fù)利用已知的LU題2

2、3 驗(yàn)證希爾伯特矩陣的病態(tài)性,矩陣基本運(yùn)算題31 用泰勒級數(shù)的有限項(xiàng)逼近正弦函數(shù)，由圖像觀察逼近效果題32 繪制飛機(jī)的降落曲線，線性方程組求解，與繪圖題41 線性擬合的函數(shù)表達(dá)式的推導(dǎo)，使用了兩種代碼方法題51 用幾種不同的方法求積分，觀察數(shù)值積分的逼近效果題55 求空間曲線L弧長。求導(dǎo)后使用符號函數(shù)積分計(jì)算題61 用歐拉公式和四階龍格-庫塔法分別求解下列初值問題，代碼搜索內(nèi)容。題64 常微分方程的解，dsolve()函數(shù)使用題82 差分法解常微分方程邊值問題，ode函數(shù)無能為力，Matlab中提供bvp解算器。 solinit=bvpinit(x,yinit,params)sol= bvps

3、olver(odefun,bcfun,solinit,options) 題83 求解圓的半徑，圓心。線性方程組解參數(shù)設(shè)計(jì)總結(jié)：（1） 算法是題目的解題核心，好的算法可以使計(jì)算更加精確 （題5.1）（2） 圖形繪制在今后的課程設(shè)計(jì)，或者是論文中可以用到。（3） 無法解決的問題，需要請教室友，或者上網(wǎng)查閱。（4） MATLAB是一個很強(qiáng)大的軟件，提供了很多內(nèi)置的數(shù)學(xué)函數(shù)，直接進(jìn)行解題。查閱資料時了解到一些MATLAB論壇。通過帖子閱讀，了解到了MATLAB在科學(xué)計(jì)算方面，模擬，圖形，視頻等起到的作用。增加了對“計(jì)算科學(xué)“的理解。指 導(dǎo) 老 師 評 語建議：從學(xué)生的工作態(tài)度、工作量、設(shè)計(jì)(論文的)創(chuàng)

4、造性、學(xué)術(shù)性、使用性及書面表達(dá)能力等方面給出評價(jià)。 簽名： 20 年 月 日 數(shù)值分析課程設(shè)計(jì)11 水手、猴子和椰子問題：五個水手帶了一只猴子來到南太平洋的一個荒島上，發(fā)現(xiàn)那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊，很快就入睡了。第一個水手醒來后，把椰子平分成五堆，將多余的一只給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個水手也陸續(xù)起來，和第一個水手一樣，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再給猴子，試問原先共有幾只椰子？（15621）試分析椰子數(shù)目的變化規(guī)律，利用逆向遞推的方法求解這一問題解：算法

5、分析：解該問題主要使用遞推算法，關(guān)于椰子數(shù)目的變化規(guī)律可以設(shè)起初的椰子數(shù)為，第一至五次猴子在夜里藏椰子后，椰子的數(shù)目分別為再設(shè)最后每個人分得x個椰子,由題: （k=0,1,2,3,4）所以,利用逆向遞推方法求解 （k=0,1,2,3,4）輸出結(jié)果 ： 1023 15621MATLAB代碼：n=input('n= ');n= 15621for x=1:n p=5*x+1; for k=1:5p=5*p/4+1; end結(jié)果分析：此題的解題思想是在迭代法中，判斷p為整數(shù)時，輸出與p if p=fix(p),breakendenddisp(x,p)12 設(shè)，（1）從盡可能精確的近似值

6、出發(fā)，利用遞推公式：計(jì)算機(jī)從到的近似值；（2）從較粗糙的估計(jì)值出發(fā)，用遞推公式：計(jì)算從到的近似值；解：首先第一步，估計(jì)和的值：syms x n; int (x0/(5+x),0,1)ans=log(2)+log(3)-log(5)eval(ans)ans=0.1823則取為0.18syms x n;int(x30/(5+x),0,1)ans =eval(ans)ans = 0MATLAB中小數(shù)點(diǎn)后保留四位，由上面計(jì)算知道積分的值不為了零。所以的取值為MATLAB代碼：i=input('i=');i=0.18;>> if i>=0.1&&i<

7、;=0.2 for n=1:1:20 i=-5*i+1/n endelseif i>0&&i<=0.0001 for n=30:-1:2 i=(-1/5)*i+1/(5*n) endend i =1.1336e+005i = -5.6679e+005i = 2.8339e+006i = -1.4170e+007i = 7.0848e+007i = -3.5424e+008i = 1.7712e+009i = -8.8560e+009i = 4.4280e+010i = -2.2140e+011同理輸入積分初始值i=0時可以得 i=0.0884結(jié)果分析：第二種方法所得

8、的結(jié)果相對來說比較精確一些，也比較可靠因?yàn)榈谝环N方法每一迭代都將最初的誤差放大了五倍，使得最終的誤差越來越大；而第一種方法經(jīng)過每一次迭代都將誤差縮小為初始誤差的五分之一，使得最終的誤差越來越小，因此相對來說比較可靠，性能較好。13 繪制Koch分形曲線問題描述：從一條直線段開始，將線段中間的三分之一部分用一個等邊三角形的另兩條邊代替，形成具有5個結(jié)點(diǎn)的新的圖形（圖1-4）；在新的圖形中，又將圖中每一直線段中間的三分之一部分都用一個等邊三角形的另兩條邊代替，再次形成新的圖形（圖1-5），這時，圖形中共有17個結(jié)點(diǎn)。這種迭代繼續(xù)進(jìn)行下去可以形成Koch分形曲線。問題分析：考慮由直線段（2個點(diǎn)）產(chǎn)生

9、第一個圖形（5個點(diǎn)）的過程，設(shè)和分別為原始直線段的兩個端點(diǎn)?，F(xiàn)在需要在直線段的中間依次插入三個點(diǎn)產(chǎn)生第一次迭代的圖形（圖1-4）。顯然，位于點(diǎn)右端直線段的三分之一處，點(diǎn)繞旋轉(zhuǎn)60度（逆時針方向）而得到的，故可以處理為向量經(jīng)正交變換而得到向量，形成算法如下：（1）；（2）；（3）；在算法的第三步中，A為正交矩陣。這一算法將根據(jù)初始數(shù)據(jù)（和點(diǎn)的坐標(biāo)），產(chǎn)生圖1-4中5個結(jié)點(diǎn)的坐標(biāo)。這5個結(jié)點(diǎn)的坐標(biāo)數(shù)組，組成一個5×2矩陣。這一矩陣的第一行為為的坐標(biāo)，第二行為的坐標(biāo)，第二行為的坐標(biāo)第五行為的坐標(biāo)。矩陣的第一列元素分別為5個結(jié)點(diǎn)的x坐標(biāo)，第二列元素分別為5個結(jié)點(diǎn)的y坐標(biāo)。問題思考與實(shí)驗(yàn)：（1

10、）考慮在Koch分形曲線的形成過程中結(jié)點(diǎn)數(shù)目的變化規(guī)律。設(shè)第k次迭代產(chǎn)生結(jié)點(diǎn)數(shù)為，第迭代產(chǎn)生結(jié)點(diǎn)數(shù)為，試寫出和之間的遞推關(guān)系式；（2）參考問題分析中的算法，考慮圖1-4到圖1-5的過程，即由第一次迭代的5個結(jié)點(diǎn)的結(jié)點(diǎn)坐標(biāo)數(shù)組，產(chǎn)生第二次迭代的17個結(jié)點(diǎn)的結(jié)點(diǎn)坐標(biāo)數(shù)組的算法；（3）考慮由第k次迭代的個結(jié)點(diǎn)的結(jié)點(diǎn)坐標(biāo)數(shù)組，產(chǎn)生第次迭代的個結(jié)點(diǎn)的結(jié)點(diǎn)坐標(biāo)數(shù)組的算法；（4）設(shè)計(jì)算法用計(jì)算機(jī)繪制出如下的Koch分形曲線（圖1-6）解：(1) (2)(3)算法及(4)代碼分析：p=0 0;10 0; %P為初始兩個點(diǎn)的坐標(biāo),第一列為x坐標(biāo),第二列為y坐標(biāo)n=2; %n為結(jié)點(diǎn)數(shù)A=cos(pi/3) -si

11、n(pi/3);sin(pi/3) cos(pi/3); %旋轉(zhuǎn)矩陣for k=1:4 d=diff(p)/3; %diff計(jì)算相鄰兩個點(diǎn)的坐標(biāo)之差,得到相鄰兩點(diǎn)確定的向量 %則d就計(jì)算出每個向量長度的三分之一,與題中將線段三等分對應(yīng) m=4*n-3; %迭代公式 q=p(1:n-1,:); %以原點(diǎn)為起點(diǎn),前n-1個點(diǎn)的坐標(biāo)為終點(diǎn)形成向量 p(5:4:m,:)=p(2:n,:); %迭代后處于4k+1位置上的點(diǎn)的坐標(biāo)為迭代前的相應(yīng)坐標(biāo) p(2:4:m,:)=q+d; %用向量方法計(jì)算迭代后處于4k+2位置上的點(diǎn)的坐標(biāo) p(3:4:m,:)=q+d+d*A' %用向量方法計(jì)算迭代后處于

12、4k+3位置上的點(diǎn)的坐標(biāo) p(4:4:m,:)=q+2*d; %用向量方法計(jì)算迭代后處于4k位置上的點(diǎn)的坐標(biāo) n=m; %迭代后新的結(jié)點(diǎn)數(shù)目end plot(p(:,1),p(:,2) %繪出每相鄰兩個點(diǎn)的連線axis(0 10 0 3.5)21 用高斯消元法的消元過程作矩陣分解。設(shè)消元過程可將矩陣A化為上三角矩陣U，試求出消元過程所用的乘數(shù)、并以如下格式構(gòu)造下三角矩陣L和上三角矩陣U驗(yàn)證：矩陣A可以分解為L和U的乘積，即A=LU。矩陣LU分解MATLAB代碼:function hl=zhjLU(A)n,n=size(A);RA=rank(A);if RA=n disp('因?yàn)锳的n階

13、行列式hl等于零，所以A不能進(jìn)行LU分解.A的秩RA如下：'); RA,hl=det(A); returnendif RA=n for p=1:n h(p)=det(A(1:p,1:p); end hl=h(1:n); for i=1:n if h(1,i)=0 disp('因?yàn)锳的各階主子式等不等于零，所以A能進(jìn)行LU分解.A的秩RA和各階順序主子式值hl依次如下：');RA,hl return end end if h(1,i)=0 disp('因?yàn)锳的各階主子式都不等于零，所以A能進(jìn)行LU分解.A的秩RA和各階順序主子式值hl如下：'); for

14、j=1:n U(1,j)=A(1,j); end for k=2:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if i>jL(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)-L(i,1:k-1)*U(1:k-1,k)/U(k,k); else U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);end end end end RA,hl,U,L· endend以上上代碼保存為M文件，并在命令窗口輸入A=20 2 3;1 8 1; 2 -3

15、15;b=0 0 0'h=zhjLU(A)程序運(yùn)行結(jié)果：L = U= 1.0000 0 0 20.0000 2.0000 3.0000 0.0500 1.0000 0 0 7.9000 0.8500 0.1000 -0.4051 1.0000 0 0 15.0443實(shí)驗(yàn)驗(yàn)證： 可以直接使用MATLAB內(nèi)置LU分解A=20 2 3;1 8 1; 2 -3 15;L U=lu(A) 輸出結(jié)果與上程序輸出結(jié)果一致。22 用矩陣分解方法求上題中A的逆矩陣。記分別求解方程組由于三個方程組系數(shù)矩陣相同，可以將分解后的矩陣重復(fù)使用。對第一個方程組，由于A=LU，所以先求解下三角方程組,再求解上三角方

16、程組，則可得逆矩陣的第一列列向量；類似可解第二、第三方程組，得逆矩陣的第二列列向量的第三列列向量。由三個列向量拼裝可得逆矩陣。解:MATLAB代碼如下：b1=1;0;0; b2=0;1;0; b3=0;1;1;A=20,2,3;1,8,1;2,-3,15;L=1,0,0;0.05,1,0;0.1,-0.4051,1;U=20 2 3;0 7.9 0.85;0 0 15.0443;Y1=Lb1 X1=UY1 Y2=Lb2X2=UY2Y3=Lb3X3=UY3Y1 = 1.0000 -0.0500 -0.1203X1 = 0.0517 -0.0055 -0.0080Y2 =實(shí)驗(yàn)驗(yàn)證： X1 X2 X

17、3ans = 0.0517 -0.0164 -0.0257 -0.0055 0.1237 0.1165 -0.0080 0.0269 0.0934而：inv（A）=X1 X2 X3 得證 0 1.0000 0.4051X2 = -0.0164 0.1237 0.0269Y3 = 0 1.0000 1.4051X3 = -0.0257 0.11650.093423驗(yàn)證希爾伯特矩陣的病態(tài)性：對于三階矩陣取右端向量，驗(yàn)證：（1）向量是方程組的準(zhǔn)確解；（2）取右端向量b的三位有效數(shù)字得，求方程組的準(zhǔn)確解，并與X的數(shù)據(jù)作比較。說明矩陣的病態(tài)性。解：(1)H=1 1/2 1/3;1/2 1/3 1/4;1

18、/3 1/4 1/5;X=1;1;1;b=H*Xb = 1.8333 1.08330.7833 與題中相同(2) 先求出解X，與數(shù)據(jù)作比較。H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5;b=1.83;1.08;0.783;X=HbX = 1.0800 0.54001.4400與相差較大，矩陣為病態(tài)矩陣31用泰勒級數(shù)的有限項(xiàng)逼近正弦函數(shù)用計(jì)算機(jī)繪出上面四個函數(shù)的圖形。解：MATLAB代碼如下(1) syms x;taylor(sin(x)x=0:0.01*pi:piplot(x,sin(x)(2) syms x;taylor(x)x=0:0.01*pi:pi/2plot

19、(x)(3) syms x;taylor(x-x3/6)fplot('x-x3/6',0 pi/2)(4) syms x;taylor(x-x3/6+x5/120)fplot('x-x3/6+x5/120',0 pi/2)結(jié)果圖形右：x=0:0.1:pi;y=sin(x);plot(x,y,'-sk');hold onx=0:0.1:pi/2;y=x;plot(x,y,'-b*')hold onfplot('x-x.3/6',0,pi/2,0,2,2e-3,'-gx')hold onfplot(&#

20、39;x-x.3/6+x.5/120',0,pi/2,0,2,2e-3,'-r.')hold offlegend('sin(x)','x','x-x3/6','x-3/6+x5/120',2)xlabel('x')ylabel('y')title('Taylor approximation')結(jié)果分析：圖中紅色點(diǎn)線為正弦曲線，藍(lán)色的星線為一階泰勒逼近，綠色叉線為二階泰勒逼近，黑色正方形線為三階泰勒逼近。可見三階泰勒逼近效果最好，泰勒級數(shù)越高，逼近效果越好。32

21、 繪制飛機(jī)的降落曲線一架飛機(jī)飛臨北京國際機(jī)場上空時，其水平速度為540km/h，飛行高度為1 000m。飛機(jī)從距機(jī)場指揮塔的橫向距離12 000m處開始降落。根據(jù)經(jīng)驗(yàn)，一架水平飛行的飛機(jī)其降落曲線是一條三次曲線。建立直角坐標(biāo)系，設(shè)飛機(jī)著陸點(diǎn)為原點(diǎn)O，降落的飛機(jī)為動點(diǎn)，則表示飛機(jī)距指揮塔的距離，表示飛機(jī)的飛行高度，降落曲線為該函數(shù)滿足條件：（1）試?yán)脻M足的條件確定三次多項(xiàng)式中的四個系數(shù)；（2）用所求出的三次多項(xiàng)式函數(shù)繪制出飛機(jī)降落曲線。function s=f(x1,y1,x2,y2,x3,y3,x4,y4)format longa1=1 ,x1,x12,x13;a2=1 ,x2,x22,x2

22、3;a3=0,1,2*x3,3*x32;a4=0,1,2*x4,3*x42;a=a1; a2;a3;a4;b=y1; y2;y3;y4;s=ab;x=-12000:250:0;y=s(3)*x.2-s(4)*x.3plot(x,y,'-d')xlabel('x')ylabel('y')以上為M文件內(nèi)容，在命令窗口鍵入f(0,0,12000,1000,0,0,12000,0)運(yùn)行結(jié)果： ans為多項(xiàng)式系數(shù)ans = 1.0e-004 * 0 0 0.20833333333333 -0.0000115740740741 曾任英特爾公司董事長的摩爾先生

23、早在1965年時，就觀察到一件很有趣的現(xiàn)象：集成電路上可容納的零件數(shù)量，每隔一年半左右就會增長一倍，性能也提升一倍。因而發(fā)表論文，提出了大大有名的摩爾定律（Moores Law），并預(yù)測未來這種增長仍會延續(xù)下去。下面數(shù)據(jù)中，第二行數(shù)據(jù)為晶片上晶體數(shù)目在不同年代與1959年時數(shù)目比較的倍數(shù)。這些數(shù)據(jù)是推出摩爾定律的依據(jù)：年代19591962196319641965增加倍數(shù)13456試從表中數(shù)據(jù)出發(fā)，推導(dǎo)線性擬合的函數(shù)表達(dá)式。解：解法一MATLAB代碼：x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p1=polyfit(x,y,1)y1=polyval(p1,x)p

24、lot(x,y1,'-',x,y,'r*')xlabel('x'),ylabel('y');運(yùn)行結(jié)果：p1 = 1.0e+003 * 0.0008 -1.6255y1 =0.8113 3.3019 4.1321 4.9623 5.7925線性擬合的函數(shù)表達(dá)式：Y=0.8302x-1.6255e+003解法二：運(yùn)行結(jié)果：xs = 1.0e+003 * -1.625528301891238 0.000830188679248x=1959 1962 1963 1964 1965;y=1;3;4;5;6;for i=1:length(x)

25、 for j=1:2 A(i,j)=x(i)(j-1); endendL,U=lu(A'*A);xs=U(L(A'*y)從而年代y與增加倍數(shù)x之間的關(guān)系為：y=-1625.528301891238+0.830188679248x51 用幾種不同的方法求積分的值。（1）牛頓-萊布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）復(fù)合梯形公式。解：syms xi1=int(4/(1+x2),x,0,1)運(yùn)行結(jié)果：i1 =pii2 = 3i3 = 3.1333i4 = 3.1416a=0;b=1;h=b-a;i2=(4/(1+a2)+4/(1+b2)/2i3=h/6*(4/(1+a2

26、)+4*4/(1+(a+b/2)2)+4/(1+b2)M=100;h=(b-a)/M;i4=0;for k=1:(M-1) x=a+h*k; i4=i4+4/(1+x2);endi4=h*(4/(1+a2)+4/(1+b2)/2+h*i4結(jié)果分析：牛頓萊布尼茲公式得到精確結(jié)果，采用梯形公式得到的結(jié)果比采用Simpson公式的精確度要低，采用復(fù)化梯形公式在步長取得越來越小的狀態(tài)下可以提高精度。55 求空間曲線L：弧長公式為解：運(yùn)行程序?yàn)椋簊yms t;x=diff(cos(t);y=diff(sin(t);z=diff(2-cos(t)-sin(t);y=int(x2+y2+z2)0.5,t,0

27、,2*pi); %符號函數(shù)積分digits(14);i=vpa(y)運(yùn)行結(jié)果：i = 8.737752570989461用歐拉公式和四階龍格-庫塔法分別求解下列初值問題；解：（1）<1>歐拉公式：function t,x=Euler(fun,t0,tt,x0,N)h=(tt-t0)/N;t=t0+0:N'*h;x(1,:)=x0'運(yùn)行結(jié)果：ans = 0 1.0000 0.0500 1.0450 0.1000 1.0878 0.1500 1.1285 0.2000 1.1676 0.2500 1.2051 0.3000 1.2413 0.3500 1.2762 0.

28、4000 1.3100 0.4500 1.3427 0.5000 1.3745 0.5500 1.4055 0.6000 1.4356 0.6500 1.4649 0.7000 1.4936 0.7500 1.5216 0.8000 1.5490 0.8500 1.5758 0.9000 1.6021 0.9500 1.6278 1.0000 1.6531for k=1:N f=feval(fun,t(k),x(k,:); f=f' x(k+1,:)=x(k,:)+h*f;end以'Euler.m'保存function f=Euler_fun(t,x)f=0.9*x./

29、(1+2*t);以'Euler_fun.m'保存function main_Eulert,x=Euler('Euler_fun',0,1,1,20);fh=dsolve('Dx=0.9*x/(1+2*t)','x(0)=1');for k=1:8 ft(k)=t(k); fx(k)=subs(fh,ft(k);endt,x以'main_Euler.m'保存輸入：main_Euler運(yùn)行結(jié)果：ans = 0 1.0000 0.1250 1.0954 0.2500 1.1807 0.3750 1.2585 0.5000

30、 1.3306 0.6250 1.3981 0.7500 1.4617 0.8750 1.5223 1.0000 1.5801<2>四階龍格-庫塔法function R=rk4(f,a,b,ya,N)%y'=f(x,y)%a,b為左右端點(diǎn)%N為迭代步長%h為步長%ya為初值h=(b-a)/N;T=zeros(1,N+1);Y=zeros(1,N+1);T=a:h:b;Y(1)=ya;for j=1:N k1=h*feval(f,T(j),Y(j); k2=h*feval(f,T(j)+h/2,Y(j)+k1/2); k3=h*feval(f,T(j)+h/2,Y(j)+k2

31、/2);k4=h*feval(f,Y(j)+h,Y(j)+k3); Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;endans=T' Y' 以'rk4.m'保存function z=f(x,y)z=0.9*y./(1+2*x);以'f.m'保存輸入：rk4('f',0,1,1,8)64 列出函數(shù)在區(qū)間0，e上的函數(shù)值表并作出它的圖象。其中，是初值問題的解。解v=dsolve('Dv*log(x)=2*x','v(0)=0','x');f=(1-log(v)*v f

32、=-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x)ezplot('f')輸出結(jié)果：f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x)7.4一個10次項(xiàng)式的系數(shù)為1 a1 a2 a9 a10=1 55 1320 181 50 157 773 902 055 341 693 0 -840 950 0 127 535 76 -106 286 40 632 880 0試用多項(xiàng)式的求根指令roots求出該10次方程的10個根，然后修改9次項(xiàng)的系數(shù)-55為-56，得新的10次方程，求解新的方程，觀察根的變化是否很顯

33、著。解： p=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>> roots(p)ans = 10.6051 + 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127ip=1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>> roots(p)ans = 21.7335 7.3501 + 7.7973i 7.3501 - 7.7973i 4.3589 + 3.3285i 4.3589 - 3.3285i 5.1831 2.4378 + 2.7974i 2.4378 - 2.7974i 0.3949 + 1