《數(shù)值分析》課程設(shè)計(jì)—作業(yè)

1、 數(shù)值分析課程設(shè)計(jì)作業(yè)（本文檔內(nèi)的有些運(yùn)行結(jié)果，限于篇幅，使文檔結(jié)構(gòu)更和諧、緊湊，已做相關(guān)的改動(dòng)，程序代碼沒(méi)變） 實(shí)驗(yàn)一1.1 水手、猴子和椰子問(wèn)題：五個(gè)水手帶了一只猴子來(lái)到南太平洋的一個(gè)荒島上，發(fā)現(xiàn)那里有一大堆椰子。由于旅途的顛簸，大家都很疲憊，很快就入睡了。第一個(gè)水手醒來(lái)后，把椰子平分成五堆，將多余的一只給了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五個(gè)水手也陸續(xù)起來(lái)，和第一個(gè)水手一樣，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再給猴子，試問(wèn)原先共有幾只椰子？試分析椰子數(shù)目的變化規(guī)律，利用逆向遞推的方法求解這

2、一問(wèn)題。解：1、 問(wèn)題分析：對(duì)于本題，比較簡(jiǎn)單，我們只需要判斷原來(lái)椰子的個(gè)數(shù)及每個(gè)人私藏了一份之后剩下的是否能被5除余1，直到最后分完。function fentao(n)a=cat(1,7);for j=n:-1:1 a(1)=j;i=1; while i fentao(20001)a = 15621 12496 9996 7996 6396 5116 4092對(duì)于第一個(gè)程序，n取2000；對(duì)于第二個(gè)程序，n取20001，就能得到我們想要的結(jié)果，即原先一共有15621個(gè)椰子，最終平均每人得4092個(gè)椰子。 1.2 當(dāng)時(shí)，選擇穩(wěn)定的算法計(jì)算積分.解：一、問(wèn)題分析：由知： 以及： 得遞推關(guān)系：

3、， 但是通過(guò)仔細(xì)觀察就能知道上述遞推公式每一步都將誤差放大十倍，即使初始誤差很小，但是誤差的傳播會(huì)逐步擴(kuò)大，也就是說(shuō)用它構(gòu)造的算法是不穩(wěn)定的，因此我們改進(jìn)上述遞推公式（算法）如下： 通過(guò)比較不難得出該誤差是逐步縮小的，即算法是穩(wěn)定的。2、 問(wèn)題求解：為了利用上面穩(wěn)定的算法，需要我們估計(jì)初值的值。因?yàn)?所以當(dāng)n=100的時(shí)，我們有： 9.090909090909091e-0049.645173882220725e-004于是可取他們的平均值，有=9.368041486564908e-004，利用上述穩(wěn)定算法，可求得相應(yīng)的值和程序如下（限于篇幅，這里只給出了前后各十個(gè)連續(xù)的數(shù)據(jù)，詳細(xì)的數(shù)據(jù)會(huì)連同作

4、業(yè)一并上交）： 積分計(jì)算對(duì)照表 1n計(jì)算值準(zhǔn)確值誤差10.0591827610.059182761020.0402929470.040292947030.0294028930.0294028931.17997e-1640.0227087130.0227087131.5278e-1650.0183339560.018333956060.015312850.015312851.13285e-1670.0131250370.013125037080.01147650.0114765090.0101930630.0101930630910.000999990.000999992.68937e-1192

5、0.000989110.000989112.71896e-10930.0009784640.0009784642.74854e-09940.0009680450.0009680452.77813e-08950.0009578450.0009578462.80771e-07960.0009478620.0009478592.83729e-06970.0009380510.0009380782.86687e-05980.0009287660.0009284970.000289646990.0009164210.000919110.0029260391000.0009368040.000909911

6、0.029556214 i=cat(1,100);j=cat(1,100);k=cat(1,100);i(100)=9.368041486564908e-004;format long;%求近似值for n=99:-1:1 i(n)=(1-exp(-n)/n-i(n+1)/10;end%求精確值for n=1:100 syms x; k=n; j(k)=int(exp(-k*x)/(exp(-x)+10),x,0,1);end%求誤差for n=1:100 k(n)=abs(i(n)-j(n)/j(n);endn=1:100;a=n;i;j;k;b=axlswrite(1_2.xls,b) f

7、ormat long min=(1-exp(-100)/(11*100),max=(1-exp(-100)/(100*(exp(-1)+10)min = 9.090909090909091e-004max = 9.645173882220725e-004 i100=(min+max)/2i100 = 9.368041486564908e-0041.3 繪制靜態(tài)和動(dòng)態(tài)的koch分形曲線問(wèn)題描述：從一條直線段開(kāi)始，將線段中間的三分之一部分用一個(gè)等邊三角形的另兩條邊代替，形成具有5個(gè)結(jié)點(diǎn)的新的圖形；在新的圖形中，又將圖中每一直線段中間的三分之一部分都用一個(gè)等邊三角形的另兩條邊代替，再次形成新的圖形，

8、這時(shí)，圖形中共有17個(gè)結(jié)點(diǎn)。這種迭代繼續(xù)進(jìn)行下去可以形成koch分形曲線。在迭代過(guò)程中，圖形中的結(jié)點(diǎn)將越來(lái)越多，而曲線最終顯示細(xì)節(jié)的多少取決于所進(jìn)行的迭代次數(shù)和顯示系統(tǒng)的分辨率。koch分形曲線的繪制與算法設(shè)計(jì)和計(jì)算機(jī)實(shí)現(xiàn)相關(guān)。圖1.1 koch曲線的形成過(guò)程解：一、算法分析： 考慮由直線段（2個(gè)點(diǎn)）產(chǎn)生第一個(gè)圖形（5個(gè)點(diǎn)）的過(guò)程。上圖中，設(shè)和分別為原始直線段的兩個(gè)端點(diǎn)，現(xiàn)需要在直線段的中間依次插入三個(gè)點(diǎn)，。顯然位于線段三分之一處，位于線段三分之二處，點(diǎn)的位置可看成是由點(diǎn)以點(diǎn)為軸心，逆時(shí)針旋轉(zhuǎn)600而得。旋轉(zhuǎn)由正交矩陣實(shí)現(xiàn)。算法根據(jù)初始數(shù)據(jù)（和點(diǎn)的坐標(biāo)），產(chǎn)生圖1中5個(gè)結(jié)點(diǎn)的坐標(biāo)。結(jié)點(diǎn)的坐標(biāo)數(shù)

9、組形成一個(gè)矩陣，矩陣的第一行為的坐標(biāo)，第二行為的坐標(biāo)，第五行為的坐標(biāo)。矩陣的第一列元素分別為5個(gè)結(jié)點(diǎn)的坐標(biāo)，第二列元素分別為5個(gè)結(jié)點(diǎn)的坐標(biāo)。進(jìn)一步考慮該曲線形成過(guò)程中結(jié)點(diǎn)數(shù)目的變化規(guī)律。設(shè)第次迭代產(chǎn)生的結(jié)點(diǎn)數(shù)為，第次迭代產(chǎn)生的結(jié)點(diǎn)數(shù)為，則和中間的遞推關(guān)系為。2、 問(wèn)題求解: 、程序代碼p=0 0;10 0; %p為初始兩個(gè)點(diǎn)的坐標(biāo),第一列為x坐標(biāo),第二列為y坐標(biāo)n=2; %n為結(jié)點(diǎn)數(shù)a=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); %旋轉(zhuǎn)矩陣for k=1:6 d=diff(p)/3; %diff計(jì)算相鄰兩個(gè)點(diǎn)的坐標(biāo)之差,得到相鄰兩點(diǎn)確定的向量 %則d就計(jì)

10、算出每個(gè)向量長(zhǎng)度的三分之一,與題中將線段三等分對(duì)應(yīng) m=4*n-3; %迭代公式 q=p(1:n-1,:); %以原點(diǎn)為起點(diǎn),前n-1個(gè)點(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ì)算迭代后處于4k+3位置上的點(diǎn)的坐標(biāo) p(4:4:m,:)=q+2*d; %用向量方法計(jì)算迭代后處于4k位置上的點(diǎn)的坐標(biāo) n=m; subplot(2,3,k) plot(p(:,1),p(:,2),r)

11、 %繪出每相鄰兩個(gè)點(diǎn)的連線axis(0,3,0,3) %迭代后新的結(jié)點(diǎn)數(shù)目 %axis offend 、運(yùn)行結(jié)果： 實(shí)驗(yàn)二2.1 小行星軌道問(wèn)題：一天文學(xué)家要確定一顆小行星繞太陽(yáng)運(yùn)行的軌道，他在軌道平面內(nèi)建立以太陽(yáng)為原點(diǎn)的直角坐標(biāo)系，在五個(gè)不同的對(duì)小行星作了五次觀察，測(cè)得軌道上五個(gè)點(diǎn)的坐標(biāo)數(shù)據(jù)（單位：萬(wàn)公里）如下表所示：p1p2p3p4p5x坐標(biāo)5360558460628596666268894y坐標(biāo)602611179169542349268894由開(kāi)普勒第一定律知，小行星軌道為一橢圓，橢圓的一般方程可表示為：現(xiàn)需要建立橢圓的方程以供研究。(1) 分別將五個(gè)點(diǎn)的數(shù)據(jù)代入橢圓一般方程中，寫(xiě)出五個(gè)

12、待定系數(shù)滿足的等式，整理后寫(xiě)出線性方程組以及方程組的系數(shù)矩陣和右端項(xiàng)b；(2) 用marlab求低階方程的指令ab求出待定系數(shù);(3) 分別用直接法、jacobi迭代法、gauss-seidel迭代法求出待定系數(shù).解：（1）、程序代碼：x=53605 58460 62859 66662 68894;y=6026 11197 16954 23492 68894;b=x.*x;2*x.*y;y.*y;2*x;2*y; %b是系數(shù)矩陣的轉(zhuǎn)置a=b;%a是系數(shù)矩陣a=vpa(a,16)syms a1 a2 a3 a4 a5;a=a1;a2;a3;a4;a5;b=-1;-1;-1;-1;-1;b=sym

13、(a*a) 、運(yùn)行結(jié)果：a = 2873496025.0, 646047460.0, 36312676.0, 107210.0, 12052.0 3417571600.0, 1309153240.0, 125372809.0, 116920.0, 22394.0 3951253881.0, 2131422972.0, 287438116.0, 125718.0, 33908.0 4443822244.0, 3132047408.0, 551874064.0, 133324.0, 46984.0 4746383236.0, 9492766472.0, 4746383236.0, 137788.0

14、, 137788.0 b = 2873496025.0*a1 + 646047460.0*a2 + 36312676.0*a3 + 107210.0*a4 + 12052.0*a5 3417571600.0*a1 + 1309153240.0*a2 + 125372809.0*a3 + 116920.0*a4 + 22394.0*a5 3951253881.0*a1 + 2131422972.0*a2 + 287438116.0*a3 + 125718.0*a4 + 33908.0*a5 4443822244.0*a1 + 3132047408.0*a2 + 551874064.0*a3 +

15、133324.0*a4 + 46984.0*a5 4746383236.0*a1 + 9492766472.0*a2 + 4746383236.0*a3 + 137788.0*a4 + 137788.0*a5 即是方程組的形式如下:-1=2873496025.0*a1 + 646047460.0*a2 + 36312676.0*a3 + 107210.0*a4 + 12052.0*a5-1=3417571600.0*a1 + 1309153240.0*a2 + 125372809.0*a3 + 116920.0*a4 + 22394.0*a5-1=3951253881.0*a1 + 21314

16、22972.0*a2 + 287438116.0*a3 + 125718.0*a4 + 33908.0*a5-1=4443822244.0*a1 + 3132047408.0*a2 + 551874064.0*a3 + 133324.0*a4 + 46984.0*a5-1=4746383236.0*a1 + 9492766472.0*a2 + 4746383236.0*a3 + 137788.0*a4 + 137788.0*a5（2） 、用指令ab求出待定系數(shù)的值如下： b=-1;-1;-1;-1;-1;a;a=ab a = 0.00000000008015133713757207070503

17、4479362305 -0.000000000099919136582879577308894350837601 -0.00000000017639794583181340912897589334276 -0.000012556065177690492529375516352892 0.000015497775048603393791640753306649 （3） 、1）、直接法 、程序代碼：function ra,rb,n,x=bilizy(a,b)b=a,b;n=length(b);ra=rank(a);rb=rank(b);rank_cha=rb-ra;if rank_cha0, di

18、sp(因?yàn)閞a!=rb，所以此方程組無(wú)解.) returnendif ra=rb if ra=n disp(因?yàn)閞a=rb=n,所以此方程組有唯一解.) x=zeros(n,1);c=zeros(1,n+1); t=0;d=b;s=max(abs(d(1:n,:); s=s; for p=1:n-1 y,j=max(abs(b(p:n,p)./s(p:n); j=min(j); c=b(p,:);b(p,:)=b(j+p-1,:);b(j+p-1,:)=c; t=s(p);s(p)=s(j+p-1);s(j+p-1)=t; for k=p+1:n m=b(k,p)/b(p,p); b(k,p:

19、n+1)=b(k,p:n+1)-m*b(p,p:n+1); end end b=b(1:n,n+1); a=b(1:n,1:n);x(n)=b(n)/a(n,n); for q=n-1:-1:1 x(q)=(b(q)-sum(a(q,q+1:n)*x(q+1:n)/a(q,q); end else disp(請(qǐng)注意：因?yàn)閞a=rbb=-1;-1;-1;-1;-1;ra,rb,n,x=bilizy(a,b)因?yàn)閞a=rb=n,所以此方程組有唯一解.ra = 5rb = 5n = 5x = 1.0e-004 * 0.0000 -0.0000 -0.0000 -0.1256 0.1550 x=vpa

20、(x,16) x = 0.00000000008015133713756955 -0.0000000000999191365828799 -0.0000000001763979458318154 -0.00001255606517769044 0.00001549777504860351 2）、jacobi迭代法 、程序代碼：%p:范數(shù)的名稱，p=1，2，inf%error:近似解a的誤差%maxl:迭代的最大次數(shù)function a=jacobi(a,b,a0,p,error,maxl)n n=size(a);a=zeros(n,1);for k=1:maxl for j=1:n a(j)=

21、(b(j)-a(j,1:j-1,j+1:n) a0(1:j-1,j+1:n)/a(j,j); end erra=norm(a-a0,p); a0=a; if(erra=error) disp(請(qǐng)注意：jacobi迭代次數(shù)已經(jīng)超過(guò)最大迭代次數(shù)maxl。)end 、運(yùn)行結(jié)果：x=53605 58460 62859 66662 68894;y=6026 11197 16954 23492 68894;b=x.*x;2*x.*y;y.*y;2*x;2*y; %b是系數(shù)矩陣的轉(zhuǎn)置a=b;%a是系數(shù)矩陣b=-1;-1;-1;-1;-1;a0=0;0;0;0;0; a=jacobi(a,b,a0,1,0.0

22、001,100)近似解a分別是：a = -0.0000000003480081375786834 -0.0000000007638525188999265 -0.000000003479009721870011 -0.000007500525036752573 -0.000007257526054518536 3）、gauss-seidel迭代法 、程序代碼：%p:范數(shù)的名稱，p=1,2,inf;error：近似解x的誤差;maxl：迭代的最大次數(shù)function a=gs(a,b,a0,p,error,maxl)n n=size(a);a=zeros(n,1);for k=1:maxl fo

23、r j=1:n aa=0; for i=1:n if ij aa=aa+a(j,i)*a0(i); end end a(j)=(b(i)-aa)/a(j,j); end erra=norm(a-a0,p); a0=a; if (erra=error) disp(請(qǐng)注意：gauss-seidel迭代次數(shù)已經(jīng)超過(guò)最大迭代次數(shù)maxl.)end 、運(yùn)行結(jié)果： x=53605 58460 62859 66662 68894;y=6026 11197 16954 23492 68894;b=x.*x;2*x.*y;y.*y;2*x;2*y; a=b;b=-1;-1;-1;-1;-1;a0=0;0;0;0

24、;0; a=gs(a,b,a0,inf,0.00001,100)近似解a分別是：a = 0.0000000001678142444990526 0.000000002024408102712696 0.000000000369101766105352 -0.000009750490974607461 -0.00015547175272548952.3 設(shè) (1) 將矩陣a進(jìn)行l(wèi)u分解a=lu，其中u是上三角矩陣，l是主對(duì)角線上的元素都是1的下三角矩陣。解： a=20 2 3;1 8 1;2 -3 15; l,u=lu(a)l = 1.0000 0 0 0.0500 1.0000 0 0.100

25、0 -0.4051 1.0000u = 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443(2) 利用上述分解分別求解方程組并由此求出逆矩陣。 b1=1;0;0;b2=0;1;0;b3=0;0;1; x1=u(lb1), x2=u(lb2),x3=u(lb3),inv(a)x1 = 0.0517 -0.0055 -0.0080解：inv(a) = 0.0517 -0.0164 -0.0093 -0.0055 0.1237 -0.0072 -0.0080 0.0269 0.0665x3 = -0.0093 -0.0072 0.0665x2 = -0

26、.0164 0.1237 0.0269(3) 用lu分解求下列線性方程組的解（方程組的精確解是） 解： a=4 2 -3 -1 2 1 0 0 0 0; 8 6 -5 -3 6 5 0 1 0 0; 4 2 -2 -1 3 2 -1 0 3 1; 0 -2 1 5 -1 3 -1 1 9 4; -4 2 6 -1 6 7 -3 3 2 3; 8 6 -8 5 7 17 2 6 -3 5; 0 2 -1 3 -4 2 5 3 0 1; 16 10 -11 -9 17 34 2 -1 2 2; 4 6 2 -7 13 9 2 0 12 4; 0 0 -1 8 -3 -24 -8 6 3 -1; b

27、=5;12;3;2;3;46;13;38;19;-21; l u=lu(a); x=u(lb)x = 1.0000 -1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000 -1.0000 2.0000結(jié)果是：2.4 (1) 用追趕法求解方程組(a) ，(b) ，解： (a)、程序代碼：%追趕法function x=zhuiganfa(a,b)n=rank(a);for i=1:n if a(i,i)=0 return; endendd=ones(n,1);a=ones(n-1,1);c=ones(n-1);for i=1:n-1 a(i,1)=a(i

28、+1,i);c(i,1)=a(i,i+1);d(i,1)=a(i,i);endd(n,1)=a(n,n);%解ly=b的解，解保存在b中for i=2:n d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1)*b(i-1,1);end%求解uy=x的解x(n,1)=b(n,1)/d(n,1);for i=n-1:-1:1 x(i,1)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1);end 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.166

29、6 0.1668 0.1662 0.1683 0.1607 0.1891 0.0829 0.4793x = 0.4793 0.0829 0.1891 0.1607 0.1683 0.1662 0.1668 0.1666 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 0.1667 運(yùn)行結(jié)果： clear all a=zeros(30,30);a(30,30)=4;b=ones(30,1);b(1,1)=2;b(30,1)=2;for i=1:29 a(i,i+1)=1;a(i+1,i)=1; a(i,i)=4;end x=zhuiganfa(a,b) a=z

30、eros(20,20);a(19,19)=5;a(20,19)=-2;a(19,20)=-2;a(20,20)=5;b=zeros(20,1);b(1,1)=1;for i=1:18 a(i,i+2)=1;a(i+2,i)=1; a(i,i+1)=-2;a(i+1,i)=-2; a(i,i)=5;end l,u=lu(a); x=vpa(u(lb),6) 0.0000504376 0.000106309 0.0000292324 -0.0000143431 -0.0000126677 -0.00000146436 0.00000249126 0.00000131198 -9.33542*10(

31、-8) -2.99738*10(-7)x = 0.242121 0.0943914 -0.0218242 -0.0313623 -0.00694256 0.00488693 0.00358612 0.000214823 -0.000784527 -0.000357862 (b)、 (2) 設(shè)計(jì)實(shí)驗(yàn)驗(yàn)證hilbert矩陣的病態(tài)性，其中format rat;for n=1:10 a=hilb(n); tiaojianshua=norm(a,1)*norm(inv(a),1)end解：tiaojianshua = 1 tiaojianshua = 27 tiaojianshua = 748 tiao

32、jianshua = 28375 tiaojianshua = 943656 tiaojianshua = 29070279 tiaojianshua =985194890 tiaojianshua =33872791987 tiaojianshua =1099651961846 tiaojianshua =35353690673203 從上面的結(jié)果可以看出，hilbert矩陣的條件數(shù)的增長(zhǎng)速度驚人，才十階，矩陣的條件數(shù)就為35353690673203，這是一個(gè)，可怕的數(shù)字，由它組成的方程組肯定是病態(tài)的。實(shí)驗(yàn)三3.1 用taylor級(jí)數(shù)的第項(xiàng)多項(xiàng)式，分別在區(qū)間和上逼近正弦函數(shù)，并用計(jì)算機(jī)繪出上

33、面31個(gè)函數(shù)的圖形。解： 一、程序代碼：syms x;y=sin(x);ezplot(y,0,pi);grid on;axis(0,pi,0,1.5);hold on;for m=1:2:30 p=taylor(y,x,m); ezplot(p,0,pi); axis(0,pi,0,1.5);end 二、函數(shù)圖形： 3.2 已知四個(gè)點(diǎn)a，b，c，d的具體位置a(0, 0)，b(0, 3), c(8, 1), d(10, 5),求兩個(gè)點(diǎn)h1，h2的具體位置，使ah1+bh1+h1h2+h2c+h2d為最短。解：function d =ju_li(a,b)d=sqrt(sum(a-b).2);re

34、turnendsyms x1 y1 x2 y2;a=0 0;b=0 3;c=8 1;d=10 5;h1=x1 y1;h2=x2 y2;sum_d=ju_li(a,h1)+ju_li(b,h1)+ju_li(h1,h2)+ju_li(c,h2)+ju_li(d,h2)x1 y1 x2 y2=solve(diff(sum_d,x1),diff(sum_d,y1),diff(sum_d,x2),diff(sum_d,y2),x1,y1,x2,y2)sum_d =3*x12 - 2*x1*x2 + 3*x22 - 36*x2 + 3*y12 - 2*y1*y2 - 6*y1 + 3*y22 - 12*

35、y2 + 199x1 =9/4y1 =27/4x2 =15/8y2 =21/8實(shí)驗(yàn)四4.1 方程的根實(shí)際上是兩個(gè)函數(shù)的交點(diǎn)的橫坐標(biāo)，用計(jì)算機(jī)繪出兩個(gè)函數(shù)在區(qū)間的圖形，觀察圖形，分析它們的交點(diǎn)分布規(guī)律及特點(diǎn)，試寫(xiě)出方程的全部實(shí)根所在的隔根區(qū)間；并用二分法求出每一個(gè)根的近似值。 解： 一、畫(huà)函數(shù)圖像如下： syms x y1 y2; y1=sin(x);y2=1/x; ezplot(y1,-6,6);grid on; hold on; ezplot(y2,-6,6) gtext(y1=sin(x) gtext(y2=1/x) 2、 用二分法求近似解: 通過(guò)上面的函數(shù)圖像可以知道，y1=sin(x)

36、與y2=1/x在區(qū)間-6,6上共有四個(gè)交點(diǎn)，即是函數(shù)y=xsin(x)在區(qū)間-6,6上共有四個(gè)零點(diǎn)，調(diào)用matlab中的函數(shù)【 ginput(4) 】 可以大概得到四個(gè)根的大概區(qū)間為：-4,-2、-2,-0.5、0.5,2、2,4，通過(guò)二分法求四個(gè)點(diǎn)的近似值的過(guò)程如下：%用二分法求非線性方程f(x)=0的根；fun為函數(shù)的表達(dá)式%a,b為左右端點(diǎn)，eps為精度，x為近似根，k為二分次數(shù)function x,k=eff(fun,a,b,eps)if nargin4 eps=1e-5;%如果輸入自變量數(shù)目0 disp(a,b不是有根區(qū)間，請(qǐng)重新調(diào)整); return;endk=0;while ab

37、s(b-a)/2eps x=(a+b)/2;fx=feval(fun,x); if fx*fa fun=(x)(x*sin(x)-1) ;eps=1e-5; x1,k1=eff(fun,-4,-2,eps),x2,k2=eff(fun,-2,-0.5,eps), x3,k3=eff(fun,0.5,2,eps),x4,k4=eff(fun,2,4,eps)x1 = -2.7726k1 = 17 x4 = 2.7726k4 = 17x3 = 1.1142k3 = 17x2 = -1.1142k2 = 174.2 設(shè)方程有3個(gè)實(shí)根盡可能多地采用下面六種不同計(jì)算格式，求的根或或的近似值，并觀察相應(yīng)格

38、式的收斂性。(1) (2) (3) (4) (5) (6) 解： syms x fun; fun=x3-3*x-1; ezplot(fun,-2.5,2.5),grid1、 函數(shù)圖像： 得到該函數(shù)的零點(diǎn)大概在：-1.5、-0.3、1.8的附近。2、 迭代格式： 程序代碼:%迭代法%用迭代法求非線性方程f(x)=0的根，fun為函數(shù)（x）的表達(dá)式%x0為初值，eps為精度（默認(rèn)1e-5），n為最大迭代次數(shù)（默認(rèn)500），x為近似根function x,k=ddf(fun,x0,eps,nmax)if nargin4 nmax=500;endif nargineps&k fun=inline(3*

39、x+1)/x2); x1,k1=ddf(fun,-1.5,1e-5),x2,k2=ddf(fun,-0.3,1e-5),x3,k3=ddf(fun,1.8,1e-5), -1.5321k1 = 28x2 = nank2 = 31x3 = nank3 = 48 x2 = nank2 = 31x3 = nank3 = 48x1 = -1.5321k1 = 28(2) fun=inline(x3-1)/3); x1,k1=ddf(fun,-1.5,1e-5),x2,k2=ddf(fun,-0.3,1e-5),x3,k3=ddf(fun,1.8,1e-5),x2 = -0.3473k2 = 5x3 =

40、 -0.3473k3 = 9x1 = -0.3473k1 = 12(3) fun=inline(3*x+1)(1/3); x1,k1=ddf(fun,-1.5,1e-5),x2,k2=ddf(fun,-0.3,1e-5),x3,k3=ddf(fun,1.8,1e-5),x3 = 1.8794k3 = 8x2 = 1.8794k2 = 12x1 =1.8794 + 0.0000ik1 = 12 fun=inline(1/(x2-3); x1,k1=ddf(fun,-1.5,1e-5),x2,k2=ddf(fun,-0.3,1e-5),x3,k3=ddf(fun,1.8,1e-5),(4)x2 = -0.3473k2 = 5x3 = -0.3473k3 = 7x1 = -0.3473k1 = 8(5) fun=inline(3+1/x)(1/3); x1,k1=ddf(fun,-1.5,1e-5),x2,k2=ddf(fun,-0.3,1e-5),x3,k3=ddf(fun,1.8,1e-5),x1 =