第二章非線性方程(組)的數(shù)值解法

1、第二章 非線性方程（組）的數(shù)值解法2.1 方程（組）的根及其MATLAB命令 求解方程（組）的solve命令求方程f(x)=q(x)的根可以用MATLAB命令:>> x=solve('方程f(x)=q(x)','待求符號(hào)變量x')求方程組fi(x1,xn)=qi(x1,xn) (i=1,2,n)的根可以用MATLAB命令:>>E1=sym('方程f1(x1,xn)=q1(x1,xn)');.En=sym('方程fn(x1,xn)=qn(x1,xn)');x1,x2,xn=solve(E1,E2,En, x1

2、,xn) 求解方程（組）的fsolve命令fsolve的調(diào)用格式： X=fsolve(F,X0)2.2 搜索根的方法及其MATLAB程序 作圖法及其MATLAB程序作函數(shù)在區(qū)間 a,b 的圖形的MATLAB程序一x=a:h:b; % h是步長(zhǎng)y=f(x); plot(x,y)grid, gtext('y=f(x)')說(shuō)明： 此程序在MATLAB的工作區(qū)輸入，運(yùn)行后即可出現(xiàn)函數(shù)的圖形.此圖形與軸交點(diǎn)的橫坐標(biāo)即為所要求的根的近似值. 區(qū)間a,b 的兩個(gè)端點(diǎn)的距離 b-a 和步長(zhǎng)h的絕對(duì)值越小，圖形越精確.作函數(shù)在區(qū)間 a,b上的圖形的MATLAB程序二將化為，其中是兩個(gè)相等的簡(jiǎn)單函

3、數(shù).x=a:h:b; y1=h(x); y2=g(x);plot(x, y1, x, y2)grid,gtext(' y1=h(x),y2=g(x)')說(shuō)明：此程序在MATLAB的工作區(qū)輸入，運(yùn)行后即可出現(xiàn)函數(shù)的圖形.兩圖形交點(diǎn)的橫坐標(biāo)即為所要求的根的近似值. 逐步搜索法及其MATLAB程序逐步搜索法的MATLAB主程序function k,r=zhubuss(a,b,h,tol)% 輸入的量- a和b是閉區(qū)間a,b的左、右端點(diǎn)；%-h是步長(zhǎng)；%-tol是預(yù)先給定的精度.% 運(yùn)行后輸出的量-k是搜索點(diǎn)的個(gè)數(shù)；% - r是方程 在a,b上的實(shí)根的近似值，其精度是tol；X=a:h

4、:b;Y=funs(X);n=(b-a)/h+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n);for k=2:n X(k)=a+k*h;Y(k)=funs(X(k); %程序中調(diào)用的funs.m為函數(shù)sk=Y(k)*Y(k-1); if sk<=0, m=m+1;r(m)=X(k); end xielv=(Y(k+1)-Y(k)*(Y(k)-Y(k-1); if (abs(Y(k)<tol)&( xielv<=0) m=m+1;r(m)=X(k); endend例 用逐步搜索法的MATLAB程序分別求方程和在區(qū)間上的根的近似值，要求精度是0.000 1.解

5、 將逐步搜索法的MATLAB程序保存名為zhubuss.m的M文件.建立M文件funs.mfunction y=funs(x) y=2.*x.3+2.*x.2-3.*x-3在MATLAB工作窗口輸入如下程序>> k,r=zhubuss(-2,2,0.001,0.0001)運(yùn)行后輸出的結(jié)果k =4001r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250即搜索點(diǎn)的個(gè)數(shù)為k =4 001，其中有5個(gè)是方程的近似根，即r = -1.224 0，-1.000 0，-1.000 0，-0.999 0，1.225 0，其精度為0.000 1.在程序中將y=2.*

6、x.3+2.*x.2-3.*x-3用y=sin(cos(2.*x.3) 代替，可得到方程在區(qū)間上的根的近似值如下r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.9230 1.3310 1.5780 1.7650 1.92002.3 二分法及其MATLAB程序2.3.2 二分法的MATLAB程序二分法的MATLAB主程序function k,x,wuca,yx=erfen(a,b,abtol)a(1)=a; b(1)=b; ya=fun(a(1); yb=fun(b(1); %程序中調(diào)用的fun.m 為函數(shù) if ya* yb>0, disp(&

7、#39;注意：ya*yb>0,請(qǐng)重新調(diào)整區(qū)間端點(diǎn)a和b.'), returnendmax1=-1+ceil(log(b-a)- log(abtol)/ log(2); % ceil是向 方向取整for k=1: max1+1a;ya=fun(a); b;yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2; k=k-1;k,a,b,x,wuca,ya,yb,yxif yx=0a=x; b=x;elseif yb*yx>0b=x;yb=yx;elsea=x; ya=yx;endif b-a< abtol , return, e

8、ndendk=max1; x; wuca; yx=fun(x);例2.3.3 確定方程x3-x+4=0的實(shí)根的分布情況，并用二分法求在開(kāi)區(qū)間 (-2,-1)內(nèi)的實(shí)根的近似值，要求精度為0.001.解 （1）先用兩種方法確定方程x3-x+4=0 的實(shí)根的分布情況。方法1 作圖法.在MATLAB工作窗口輸入如下程序>>x=-4:0.1:4; y=x.3-x +4; plot(x,y) grid,gtext('y=x3-x+4') 畫出函數(shù)f(x)=x3-x+4的圖像.從圖像可以看出，此曲線有兩個(gè)駐點(diǎn)都在x軸的上方，在(-2,-1)內(nèi)曲線與x軸只有一個(gè)交點(diǎn)，則該方程有唯一

9、一個(gè)實(shí)根，且在 (-2,-1)內(nèi).方法2 試算法.在MATLAB工作窗口輸入程序>>x=-4: 1:4,y=x.3-x +4運(yùn)行后輸出結(jié)果x = -4 -3 -2 -1 0 1 2 3 4 y = -56 -20 -2 4 4 4 10 28 64 由于連續(xù)函數(shù)f(x)滿足，所以此方程在(-2,-1)內(nèi)有一個(gè)實(shí)根.（2） 用二分法的主程序計(jì)算.在MATLAB工作窗口輸入程序>> k,x,wuca,yx=erfen (-2,-1,0.001)運(yùn)行后屏幕顯示用二分法計(jì)算過(guò)程被列入表 2-3，其余結(jié)果為k = 9，x=-1.7959， wuca = 9.7656e-004，y

10、x = 0.0037表 2-3次數(shù)k左端點(diǎn)ak右端點(diǎn)bk中點(diǎn)xk函數(shù)值f(ak)函數(shù)值f(bk)函數(shù)值f(xk)0-2.000 0-1.000 0-1.500 00.500 0-2.000 04.000 02.125 01-2.000 0-1.500 0-1.750 00.250 0-2.000 02.125 00.390 62-2.000 0-1.750 0-1.875 00.125 0-2.000 00.390 6-0.716 83-1.875 0-1.750 0-1.812 50.062 5-0.716 80.390 6-0.141 84-1.812 5-1.750 0-1.781 30

11、.031 3-0.141 80.390 60.129 65-1.812 5-1.781 3-1.796 90.015 6-0.141 80.129 6-0.004 86-1.796 9-1.781 3-1.789 10.007 8-0.004 80.129 60.062 77-1.796 9-1.789 1-1.793 00.003 9-0.004 80.062 70.029 08-1.796 9-1.793 0-1.794 90.002 0-0.004 80.029 00.012 19-1.796 9-1.794 9-1.795 90.001 0-0.004 80.012 10.003 72

12、.4 迭代法及其MATLAB程序2.4.2 迭代法的MATLAB程序1迭代法的MATLAB主程序1function k,piancha,xdpiancha,xk=diedai1(x0,k)% 輸入的量-x0是初始值,k是迭代次數(shù)x(1)=x0; for i=1:k x(i+1)=fun1(x(i);%程序中調(diào)用的fun1.m為函數(shù)y=(x) piancha= abs(x(i+1)-x(i); xdpiancha=piancha/( abs(x(i+1)+eps); i=i+1;xk=x(i);(i-1) piancha xdpiancha xkendif (piancha >1)&

13、;(xdpiancha>0.5)&(k>3) disp('請(qǐng)用戶注意：此迭代序列發(fā)散，請(qǐng)重新輸入新的迭代公式') return; end if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝賀您！此迭代序列收斂，且收斂速度較快') return; endp=(i-1) piancha xdpiancha xk'例 求方程的一個(gè)正根.解 在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk= dieda

14、i1(2,5)運(yùn)行后輸出用迭代公式的結(jié)果k,piancha,xdpiancha,xk=1.00000000000000 1.00000000000000 0.33333333333333 3.000000000000002.00000000000000 2.50000000000000 5.00000000000000 0.500000000000003.00000000000000 4.37500000000000 0.89743589743590 4.875000000000004.00000000000000 11.75781250000000 1.70828603859251 -6.8

15、8281250000000請(qǐng)用戶注意：此迭代序列發(fā)散，請(qǐng)重新輸入新的迭代公式k=5，piancha = 11.80374145507813，xk = -18.68655395507813由以上運(yùn)行后輸出的迭代序列與根=2.316 624 790 355 40相差越來(lái)越大，即迭代序列發(fā)散，此迭代法就失敗.這時(shí)偏差piancha逐漸增大且偏差的相對(duì)誤差xdpiancha的值大于0.5.用迭代公式運(yùn)行后輸出的結(jié)果k,piancha,xdpiancha,xk= 1.00000000000000 0.50000000000000 0.20000000000000 2.50000000000000可見(jiàn)，偏

16、差piancha和偏差的相對(duì)誤差xdpiancha的值逐漸變小，且第5次的迭代值xk =2.331 460 674 157 30與根=2.316 624 790 355 40接近，則迭代序列收斂，但收斂速度較慢，此迭代法較為成功.用迭代公式運(yùn)行后輸出的結(jié)果k,piancha,xdpiancha,xk= 祝賀您！此迭代序列收斂，且收斂速度較快.可見(jiàn)，偏差piancha和偏差的相對(duì)誤差xdpiancha的值越來(lái)越小，且第5次的迭代值xk=2.316 624 790 355 40與根=2.316 624 790 355 40相差無(wú)幾，則迭代序列收斂，且收斂速度很快，此迭代法成功.2.4.5 迭代法的

17、MATLAB程序2迭代法的MATLAB主程序2function k,piancha,xdpiancha,xk,yk=diedai2(x0,tol,ddmax)x(1)=x0;for i=1: ddmax x(i+1)=fun(x(i);piancha=abs(x(i+1)-x(i); xdpiancha=piancha/( abs(x(i+1)+eps);i=i+1; xk=x(i);yk=fun(x(i); (i-1) piancha xdpiancha xk yk if (piancha<tol)|(xdpiancha< tol) k=i-1; xk=x(i); return;

18、endend if i>ddmax disp('迭代次數(shù)超過(guò)給定的最大值ddmax') k=i-1; xk=x(i);yk=fun(x(i); (i-1) piancha xdpiancha xk yk;return; endP=(i-1),piancha,xdpiancha,xk,yk'例 求x5-3x+1=0在0.3附近的根，精確到4位小數(shù).解 構(gòu)造迭代公式   .利用迭代法的MATLAB程序2計(jì)算精確到4位小數(shù)的根的近似值 y 和迭代次數(shù)k.在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk

19、=diedai2(0.3,1e-4,100)運(yùn)行后輸出的結(jié)果k,piancha,xdpiancha,xk,yk =0 0.03414 0.10218 0.30000 0.334141 0.03414 0.10218 0.33414 0.334722 0.00058 0.00173 0.33472 0.334733 0.00001 0.00004 0.33473 0.33473k =3;piancha =1.206089525390697e-005; xk =0.3347; yk =0.3347.2.5 迭代過(guò)程的加速方法及其MATLAB程序 2.5.2 加權(quán)迭代法的MATLAB程序加權(quán)迭代法的

20、MATLAB主程序現(xiàn)提供名為jasudd.m的M文件如下：function k,xk,yk=jasudd (x0,tol,L,ddmax)x1(1)=x0;for i=1: ddmax x(i+1)=fun(x1(i);x1(i+1)= x(i+1)+L*( x(i+1)- x1(i)/(1-L); piancha=abs(x1(i+1)-x1(i); xdpiancha= piancha/( abs(x1(i+1)+eps); i=i+1;xk=x1(i);yk=fun(x1(i);(i-1) xk yk if (piancha<tol)|(xdpiancha< tol) k=i

21、-1; xk=x1(i); return;endendif i>ddmax disp('迭代次數(shù)超過(guò)給定的最大值ddmax') k=i-1; xk=x1(i); (i-1) xk yk; return;endP=(i-1),xk,yk'例 求2e在0.85附近的一個(gè)近似根，要求精度.解 在MATLAB工作窗口輸入程序 >> k,xk,yk=jasudd (0.85,1e-6,-0.855,100)運(yùn)行后輸出結(jié)果k,xk,yk = 1.00000000000000 0.85260370028041 0.85260703830561 2.000000000

22、00000 0.85260549975491 0.85260550406236 3.00000000000000 0.85260550207699 0.85260550208255k =3； xk =0.852606； yk = 0.852606. 艾特肯(Aitken)加速方法的MATLAB程序艾特肯(Aitken)加速方法的MATLAB主程序現(xiàn)提供名為Aitken.m的M文件如下：function k,xk,yk,p= Aitken (x0,tol, ddmax)x(1)=x0;for i=1: ddmax x1(i+1)=fun(x(i); x2(i+1)=fun(x1(i+1); x(

23、i+1)=x2(i+1)-(x2(i+1)-x1(i+1)2/(x2(i+1)-2*x1(i+1)+ x(i); piancha=abs(x(i+1)-x(i); xdpiancha= piancha/( abs(x(i+1)+eps); i=i+1; xk=x(i);yk=fun(x(i); if (piancha<tol)|(xdpiancha<tol) k=i-1; xk=x(i); yk=fun(x(i); m=0,1:i-1; p=m',x1',x2',x' return;endendif i>ddmax disp('迭代次數(shù)

24、超過(guò)給定的最大值ddmax') k=i-1; xk=x(i); yk=fun(x(i); m=0,1:i-1; p=m',x1',x2',x' return;endm=0,1:i-1; p=m',x1',x2',x'例 用艾特肯加速方法求2e在0.85附近的一個(gè)近似根，要求精度.解 在MATLAB工作窗口輸入程序>> k,xk,yk,p= Aitken(0.85,1e-6, 100)運(yùn)行后輸出結(jié)果k=3，xk=0.85260550201343，yk=0.85260550201373p = 0 0 0 0.850

25、00000000000 1.00000000000000 0.85482986389745 0.85071110652484 0.85260683568607 2.00000000000000 0.85260436491811 0.85260647150826 0.85260550201407 3.00000000000000 0.85260550201343 0.85260550201398 0.852605502013732.6 牛頓(Newton)切線法及其MATLAB程序 牛頓切線法的收斂性及其MATLAB程序牛頓切線法的收斂性及其MATLAB主程序function y,f=newju

26、shou(x)f=fnq(x); fz=fnq(x)*ddfnq(x)/(dfnq(x)2+eps); y=abs(fz);if (y<1) disp('恭喜您！此迭代序列收斂,(x)導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下')else disp('請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值')end P=y,f'例 用牛頓切線法的局部收斂性判別方程 e的近似根時(shí)，由下列初始值產(chǎn)生的迭代序列是否收斂？ ；.解 在MATLAB工作窗口輸入程序>>

27、y,f=newjushou(-1)運(yùn)行后輸出結(jié)果請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值y = 139.5644f = 4.3096（2）輸入程序 >> y,f=newjushou(0)運(yùn)行后輸出結(jié)果請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)的值fy = 8.0000f = 4（3）輸入程序>> y,f=newjushou(1)運(yùn)行后輸出結(jié)果恭喜您！此迭代序列收斂,(x)導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下y = 0.3

28、566f = 1.7126（4）輸入程序：>> y,f=newjushou(2)運(yùn)行后輸出結(jié)果請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值y = 1.2593f = -2.7188（5）輸入程序 >> y,f=newjushou(5.5)運(yùn)行后輸出結(jié)果請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值y = 1.0447e+005f = 176.6400（6）輸入程序 >> y,f=newjushou(8)運(yùn)行后輸出結(jié)果恭喜您！此迭代序列收斂,(x)導(dǎo)數(shù)值的

29、絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下y = 0.4038f = -2.9452e+003 牛頓切線法的MATLAB程序牛頓切線法的MATLAB主程序現(xiàn)提供名為newtonqx.m的M文件：function k,xk,yk,piancha,xdpiancha=newtonqx(x0,tol,ftol,gxmax)x(1)=x0; for i=1: gxmax x(i+1)=x(i)-fnq(x(i)/(dfnq(x(i)+eps); piancha=abs(x(i+1)-x(i); xdpiancha= piancha/( abs(x(i+1)+eps); i=i

30、+1;xk=x(i);yk=fnq(x(i); (i-1) xk yk piancha xdpianchaif (abs(yk)<ftol)&(piancha<tol)|(xdpiancha< tol) k=i-1; xk=x(i);(i-1) xk yk piancha xdpiancha return;endend if i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax。') k=i-1; xk=x(i);(i-1) xk yk piancha xdpiancha return;end (i-1),xk,yk,pianc

31、ha,xdpiancha'例 用牛頓切線法求方程在的近似根，要求精度.解 在MATLAB工作窗口輸入程序>> k,xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001, 0.001,100)運(yùn)行后輸出初始值的迭代結(jié)果. 求的方法及其MATLAB程序求的方法及其MATLAB主程序現(xiàn)提供名為kai2fang.m的M文件：function k,xk,yk,piancha,xdpiancha,P=kainfang(x0,c,n,tol, gxmax)x(1)=x0;for i=1: gxmaxu(i)= (x(i)n-c)/(n*x(i)(n-1)

32、; x(i+1)= x(i)-u(i); piancha=abs(x(i+1)-x(i); xdpiancha=piancha/( abs(x(i+1)+eps);i=i+1; xk=x(i);yk=fnq(x(i);(i-1),xk,yk,piancha,xdpiancha if (piancha<tol)|(xdpiancha< tol) k=i-1;xk=x(i); yk=fnq(x(i); return;endendif i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax.') k=i-1;xk=x(i); yk=fnq(x(i);

33、(i-1),xk,yk,piancha,xdpiancha return;endP=(i-1),xk,yk,piancha,xdpiancha'例 求，要求精度為.解 本題介紹四種解法.方法1 用求的次根（當(dāng)n是偶數(shù)時(shí)，）的MATLAB程序計(jì)算.取初始值，根指數(shù)，被開(kāi)方數(shù)，近似根的精度tol=，迭代的最大次數(shù)gxmax=100.在工作區(qū)間輸入程序>> k,xk,yk,piancha,xdpiancha=kainfang(10,113,2,1e-5,100)運(yùn)行后輸出結(jié)果可見(jiàn),10.630 15，滿足精度.方法2 用牛頓迭代公式（2.12）計(jì)算.設(shè)，則，記，.由牛頓迭代公式得

34、，即取初始值，計(jì)算結(jié)果列入表 2-12. 表 2-12迭代次數(shù)偏差根的近似值10.650 00010.650 00020.019 83610.630 16430.000 01910.630 14640.000 00010.630 146因?yàn)?，迭代次?shù)=4時(shí)，偏差，滿足精度，所以，10.630 15.方法3 用牛頓切線法的MATLAB主程序計(jì)算.分別建立名為fnq.m和dfnq.m的M文件function y=fnq(x) y=x2-113; function y=dfnq(x) y=2*x;在MATLAB工作窗口輸入程序>> k,xk,yk,piancha,xdpiancha=ne

35、wtonqx(10, 1e-5, 1e-5,100)運(yùn)行后，將輸出的結(jié)果列入下表 2-13. 迭代k=4次，得到精度為的結(jié)果 10.630 15.表 2-13kpianchaxdpianchaxkyk10.650 0000.061 03310.650 0000.422 50020.019 8360.001 86610.630 1640.000 39330.000 0190.000 00210.630 1460.000 00040.000 0000.000 00010.630 1460.000 000方法4 在MATLAB工作空間輸入程序>> 113(1/2) 運(yùn)行后輸出經(jīng)過(guò)四舍五入

36、后，得到精度為的結(jié)果10.630 15. 牛頓切線法的加速及其兩種MATLAB程序（一） 已知方程根的重?cái)?shù)已知方程根的重?cái)?shù),求重根的修正牛頓切線法的MATLAB主程序現(xiàn)提供名為newtonxz.m的M文件：function k,piancha,xdpiancha,xk,yk=newtonxz(m,x0,tol,ftol,gxmax)x(1)=x0;for i=1: gxmaxx(i+1)=x(i)-m*fnq(x(i)/(dfnq(x(i)+eps); piancha=abs(x(i+1)-x(i);xdpiancha=piancha/( abs(x(i+1)+eps); i=i+1; xk=

37、x(i);yk=fnq(x(i); (i-1) piancha xdpiancha xk yk; if (piancha<tol)|(xdpiancha< tol)&(abs(yk)<ftol) k=i-1; xk=x(i); yk=fnq(x(i); (i-1) piancha xdpiancha xk yk; return;endend if i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax.') k=i-1; xk=x(i); yk=fnq(x(i); (i-1) piancha xdpiancha xk yk;ret

38、urn;end例 判斷是方程e的幾重根？在區(qū)間上，分別用牛頓切線法和求重根的修正牛頓切線公式求此根的近似值，使精確到.解 經(jīng)判斷知，是方程的三重根.在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk= newtonqx (0.5,1e-4, 1e-4,100)運(yùn)行后整理結(jié)果列入表 2-20.表 2-20kpianchaxdpianchaxkyk10.144 100.404 890.355 900.541 9820.107 410.432 250.248 490.168 2030.077 450.452 800.171 040.051 4640.0

39、54 500.467 600.116 550.015 5850.037 690.477 980.078 850.004 6960.025 760.485 140.053 100.001 4070.017 460.490 010.035 630.000 4280.011 770.493 300.023 860.000 1290.007 910.495 520.015 960.000 04100.005 300.497 000.010 660.000 01110.003 540.498 000.007 120.000 00120.002 370.498 670.004 750.000 00130.

40、001 580.499 110.003 170.000 00140.001 050.499 410.002 110.000 00150.000 700.499 600.001 410.000 00160.000 470.499 740.000 940.000 00170.000 310.499 820.000 630.000 00180.000 210.499 880.000 420.000 00190.000 140.499 920.000 280.000 00200.000 090.499 950.000 190.000 00迭代次數(shù)k=20，精確到的根的近似值是xk =0.000 2，其

41、函數(shù)值是yk =5.752 0e-011，收斂速度較慢.根據(jù)求重根的修正牛頓切線公式，在MATLAB工作窗口輸入程序 >>k,piancha,xdpiancha,xk,yk=newtonxz(3,0.5,1e-4, 1e-4,100)運(yùn)行后整理結(jié)果得表 2-21.表 2-21kpianchaxdpianchaxkyk10.432 306.385 790.067 700.002 9420.066 5457.310 850.001 160.000 0030.001 163443.447 270.000 000.000 0040.000 431.000 78-0.000 43-0.000

42、 0050.000 439228.792 130.000 00-0.000 0060.011 151.000 000.011 150.000 0170.011 12356.869 170.000 030.000 0080.000 031638.927 030.000 000.000 00迭代次數(shù)k=8，精確到的根的近似值是xk =1.900 3e-008，其函數(shù)值是yk = 4.440 9e-016，piancha = 3.114 5e-005，xdpiancha=1638.927 0，是二階收斂（平方收斂）.可見(jiàn)，求重根的修正牛頓切線公式比牛頓切線法收斂速度快得多.（二） 未知方程根的重?cái)?shù)未

43、知方程根的重?cái)?shù)為，求重根的修正牛頓切線法的MATLAB主程序現(xiàn)提供名為newtonxz1.m的M文件function k,piancha,xdpiancha,xk,yk=newtonxz1(x0,tol,ftol,gxmax)x(1)=x0;for i=1: gxmaxu(i)=fnq(x(i)/dfnq(x(i);du(i)=1-fnq(x(i)*ddfnq(x(i)/(dfnq(x(i)2+eps);x(i+1)=x(i)-u(i)/du(i); piancha=abs(x(i+1)-x(i);xdpiancha=piancha/( abs(x(i+1)+eps); i=i+1; xk=x

44、(i);yk=fnq(x(i); if (piancha<tol)|(xdpiancha< tol)&(abs(yk)<ftol) k=i-1; xk=x(i); (i-1) piancha xdpiancha xk yk return;endendif i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax.')k=i-1; xk=x(i); yk=fnq(x(i); (i-1) piancha xdpiancha xk yk ;return;end例 用未知重?cái)?shù)的求重根的修正牛頓切線法，求方程 在區(qū)間上的根，使精確到.解 在M

45、ATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk=newtonxz1(0.5,1e-4, 1e-4,100)運(yùn)行后整理結(jié)果得表 2-22. 表 2-22kpianchaxdpianchaxkyk10.599 236.038 57-0.099 23-0.008 1820.096 9843.054 85-0.002 25-0.000 0030.002 251 778.327 84-0.000 000400-0.000 0002.7 割線法及其MATLAB程序 割線法的MATLAB程序割線法的MATLAB主程序現(xiàn)提供名為gexian.m的M文件：func

46、tion k,piancha,xdpiancha,xk,yk=gexian (x01,x02,tol,ftol,gxmax)x(1)=x01;x(2)=x02;for i=2: gxmax u(i)= fnq(x(i)*(x(i)-x(i-1); v(i)= fnq(x(i)-fnq(x(i-1); x(i+1)=x(i)- u(i)/( v(i); piancha=abs(x(i+1)-x(i); xdpiancha= piancha/( abs(x(i+1)+eps); i=i+1; xk= x(i);yk=fnq(x(i); (i-2) piancha xdpiancha xk yk i

47、f (abs(yk)<ftol)&( piancha <tol)|(xdpiancha< tol) k=i-2; xk=x(i);yk=fnq(x(i); (i-2) piancha xdpiancha xk yk； return;endend if i>gxmaxdisp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax.')k=i-2; xk=x(i);yk=fnq(x(i); return;end2.8 拋物線法及其MATLAB程序 拋物線法的MATLAB程序拋物線法的MATLAB主程序function k,piancha,xdpiancha,

48、xk,yk=paowu(x1,x2, x3,tol,ftol,gxmax)X=x1, x2, x3;for i=1: gxmaxh0= X(1)- X (3); h1= X (2)- X (3); u0= fnq (X (1)- fnq (X (3) ); u1= fnq (X (2)- fnq (X (3); c= fnq (X (3);fenmu= h02* h1- h12* h0; afenzi= u0* h1- u1* h0;bfenzi= u1*h02-u0*h12; a= afenzi/ fenmu; b= bfenzi/ fenmu; deta=b2-4*a*c; cha=abs

49、(X(3)-X(1),abs(X(3)-X(2),abs(X(2)- X(1);piancha =min(cha); xdpiancha= piancha/( abs (X(3)+eps);if deta<0disp('提示：根的判別式值為負(fù)數(shù).'),detakf=sqrt(deta);elseif deta<=0disp('提示：根的判別式值為零.'),detakf=0;elsedisp('提示：根的判別式值為正數(shù).'),detakf=sqrt(deta);endif b<0 detakf=- detakf;endz=-2*c

50、/(b+ detakf);xk= X(3)+ z;q1= xk - X (1); q2= xk - X (2); q3= xk - X (3);if abs(q2)< abs(q1)X1=X(2), X(1), X(3); X= X1;Y= fnq(X);endif abs(q3)< abs(q1)X2=X(1), X(3), X(2); X= X2;Y= fnq(X);endX(3)= xk; Y(3)=fnq(X(3);yk= Y(3); i piancha xdpiancha xk yk if (abs(yk)<ftol)&( piancha <tol)|

51、(xdpiancha< tol) k=i; X(3)= xk; Y(3)=fnq(X(3);yk= Y(3); return;endendif i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax,請(qǐng)重新輸入初始值x1,x2,x3') return;endP=i,piancha,xdpiancha,xk,yk' 例 用拋物線法求方程e的一個(gè)實(shí)根的近似值，使精確到.解 先用作圖法確定初始值x01，x02和x03，然后在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk= paowu (-0.4,-

52、0.3, -0.5,1e-4, 1e-4,100)運(yùn)行后輸出結(jié)果 提示：根的判別式值為正數(shù).ans = Columns 1 through 4 Column 5 -0.00005581900299提示：根的判別式值為正數(shù).ans = Columns 1 through 4 2.00000000000000 0.00933649437761 0.02389906977038 -0.39064638365394 Column 5 -0.00000000923532提示：根的判別式值為正數(shù).ans = Columns 1 through 4 3.00000000000000 0.00001712196845 0.00004382983990 -0.39064638082059 Column 5 0.00000000000000k = 3piancha =xdpianc