1、第二章 非線性方程（組）的數值解法2.1 方程（組）的根及其MATLAB命令 求解方程（組）的solve命令求方程f(x)=q(x)的根可以用MATLAB命令:>> x=solve('方程f(x)=q(x)','待求符號變量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的調用格式： X=fsolve(F,X0)2.2 搜索根的方法及其MATLAB程序 作圖法及其MATLAB程序作函數在區間 a,b 的圖形的MATLAB程序一x=a:h:b; % h是步長y=f(x); plot(x,y)grid, gtext('y=f(x)')說明： 此程序在MATLAB的工作區輸入，運行后即可出現函數的圖形.此圖形與軸交點的橫坐標即為所要求的根的近似值. 區間a,b 的兩個端點的距離 b-a 和步長h的絕對值越小，圖形越精確.作函數在區間 a,b上的圖形的MATLAB程序二將化為，其中是兩個相等的簡單函

3、數.x=a:h:b; y1=h(x); y2=g(x);plot(x, y1, x, y2)grid,gtext(' y1=h(x),y2=g(x)')說明：此程序在MATLAB的工作區輸入，運行后即可出現函數的圖形.兩圖形交點的橫坐標即為所要求的根的近似值. 逐步搜索法及其MATLAB程序逐步搜索法的MATLAB主程序function k,r=zhubuss(a,b,h,tol)% 輸入的量- a和b是閉區間a,b的左、右端點；%-h是步長；%-tol是預先給定的精度.% 運行后輸出的量-k是搜索點的個數；% - r是方程 在a,b上的實根的近似值，其精度是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); %程序中調用的funs.m為函數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程序分別求方程和在區間上的根的近似值，要求精度是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)運行后輸出的結果k =4001r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250即搜索點的個數為k =4 001，其中有5個是方程的近似根，即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) 代替，可得到方程在區間上的根的近似值如下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); %程序中調用的fun.m 為函數 if ya* yb>0, disp(&

7、#39;注意：ya*yb>0,請重新調整區間端點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的實根的分布情況，并用二分法求在開區間 (-2,-1)內的實根的近似值，要求精度為0.001.解 （1）先用兩種方法確定方程x3-x+4=0 的實根的分布情況。方法1 作圖法.在MATLAB工作窗口輸入如下程序>>x=-4:0.1:4; y=x.3-x +4; plot(x,y) grid,gtext('y=x3-x+4') 畫出函數f(x)=x3-x+4的圖像.從圖像可以看出，此曲線有兩個駐點都在x軸的上方，在(-2,-1)內曲線與x軸只有一個交點，則該方程有唯一

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

10、x = 0.0037表 2-3次數k左端點ak右端點bk中點xk函數值f(ak)函數值f(bk)函數值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是迭代次數x(1)=x0; for i=1:k x(i+1)=fun1(x(i);%程序中調用的fun1.m為函數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('請用戶注意：此迭代序列發散，請重新輸入新的迭代公式') return; end if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝賀您！此迭代序列收斂，且收斂速度較快') return; endp=(i-1) piancha xdpiancha xk'例 求方程的一個正根.解 在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk= dieda

14、i1(2,5)運行后輸出用迭代公式的結果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請用戶注意：此迭代序列發散，請重新輸入新的迭代公式k=5，piancha = 11.80374145507813，xk = -18.68655395507813由以上運行后輸出的迭代序列與根=2.316 624 790 355 40相差越來越大，即迭代序列發散，此迭代法就失敗.這時偏差piancha逐漸增大且偏差的相對誤差xdpiancha的值大于0.5.用迭代公式運行后輸出的結果k,piancha,xdpiancha,xk= 1.00000000000000 0.50000000000000 0.20000000000000 2.50000000000000可見，偏

16、差piancha和偏差的相對誤差xdpiancha的值逐漸變小，且第5次的迭代值xk =2.331 460 674 157 30與根=2.316 624 790 355 40接近，則迭代序列收斂，但收斂速度較慢，此迭代法較為成功.用迭代公式運行后輸出的結果k,piancha,xdpiancha,xk= 祝賀您！此迭代序列收斂，且收斂速度較快.可見，偏差piancha和偏差的相對誤差xdpiancha的值越來越小，且第5次的迭代值xk=2.316 624 790 355 40與根=2.316 624 790 355 40相差無幾，則迭代序列收斂，且收斂速度很快，此迭代法成功.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('迭代次數超過給定的最大值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位小數.解 構造迭代公式   .利用迭代法的MATLAB程序2計算精確到4位小數的根的近似值 y 和迭代次數k.在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk

19、=diedai2(0.3,1e-4,100)運行后輸出的結果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 迭代過程的加速方法及其MATLAB程序 2.5.2 加權迭代法的MATLAB程序加權迭代法的

20、MATLAB主程序現提供名為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('迭代次數超過給定的最大值ddmax') k=i-1; xk=x1(i); (i-1) xk yk; return;endP=(i-1),xk,yk'例 求2e在0.85附近的一個近似根，要求精度.解 在MATLAB工作窗口輸入程序 >> k,xk,yk=jasudd (0.85,1e-6,-0.855,100)運行后輸出結果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主程序現提供名為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('迭代次數

24、超過給定的最大值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附近的一個近似根，要求精度.解 在MATLAB工作窗口輸入程序>> k,xk,yk,p= Aitken(0.85,1e-6, 100)運行后輸出結果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)導數值的絕對值y=|d(x)/dx|和方程f(x)=0的函數f(x)值f如下')else disp('請注意觀察下面顯示的(x)的導數值的絕對值y=|d(x)/dx|和方程f(x)=0的函數f(x)值')end P=y,f'例 用牛頓切線法的局部收斂性判別方程 e的近似根時，由下列初始值產生的迭代序列是否收斂？ ；.解 在MATLAB工作窗口輸入程序>>

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

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

29、絕對值y=|d(x)/dx|和方程f(x)=0的函數f(x)值f如下y = 0.4038f = -2.9452e+003 牛頓切線法的MATLAB程序牛頓切線法的MATLAB主程序現提供名為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('請注意：迭代次數超過給定的最大值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)運行后輸出初始值的迭代結果. 求的方法及其MATLAB程序求的方法及其MATLAB主程序現提供名為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('請注意：迭代次數超過給定的最大值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 用求的次根（當n是偶數時，）的MATLAB程序計算.取初始值，根指數，被開方數，近似根的精度tol=，迭代的最大次數gxmax=100.在工作區間輸入程序>> k,xk,yk,piancha,xdpiancha=kainfang(10,113,2,1e-5,100)運行后輸出結果可見,10.630 15，滿足精度.方法2 用牛頓迭代公式（2.12）計算.設，則，記，.由牛頓迭代公式得

34、，即取初始值，計算結果列入表 2-12. 表 2-12迭代次數偏差根的近似值10.650 00010.650 00020.019 83610.630 16430.000 01910.630 14640.000 00010.630 146因為，迭代次數=4時，偏差，滿足精度，所以，10.630 15.方法3 用牛頓切線法的MATLAB主程序計算.分別建立名為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)運行后，將輸出的結果列入下表 2-13. 迭代k=4次，得到精度為的結果 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) 運行后輸出經過四舍五入

36、后，得到精度為的結果10.630 15. 牛頓切線法的加速及其兩種MATLAB程序（一） 已知方程根的重數已知方程根的重數,求重根的修正牛頓切線法的MATLAB主程序現提供名為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('請注意：迭代次數超過給定的最大值gxmax.') k=i-1; xk=x(i); yk=fnq(x(i); (i-1) piancha xdpiancha xk yk;ret

38、urn;end例 判斷是方程e的幾重根？在區間上，分別用牛頓切線法和求重根的修正牛頓切線公式求此根的近似值，使精確到.解 經判斷知，是方程的三重根.在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk= newtonqx (0.5,1e-4, 1e-4,100)運行后整理結果列入表 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迭代次數k=20，精確到的根的近似值是xk =0.000 2，其

41、函數值是yk =5.752 0e-011，收斂速度較慢.根據求重根的修正牛頓切線公式，在MATLAB工作窗口輸入程序 >>k,piancha,xdpiancha,xk,yk=newtonxz(3,0.5,1e-4, 1e-4,100)運行后整理結果得表 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迭代次數k=8，精確到的根的近似值是xk =1.900 3e-008，其函數值是yk = 4.440 9e-016，piancha = 3.114 5e-005，xdpiancha=1638.927 0，是二階收斂（平方收斂）.可見，求重根的修正牛頓切線公式比牛頓切線法收斂速度快得多.（二） 未知方程根的重數未

43、知方程根的重數為，求重根的修正牛頓切線法的MATLAB主程序現提供名為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('請注意：迭代次數超過給定的最大值gxmax.')k=i-1; xk=x(i); yk=fnq(x(i); (i-1) piancha xdpiancha xk yk ;return;end例 用未知重數的求重根的修正牛頓切線法，求方程 在區間上的根，使精確到.解 在M

45、ATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk=newtonxz1(0.5,1e-4, 1e-4,100)運行后整理結果得表 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主程序現提供名為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('請注意：迭代次數超過給定的最大值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('提示：根的判別式值為負數.'),detakf=sqrt(deta);elseif deta<=0disp('提示：根的判別式值為零.'),detakf=0;elsedisp('提示：根的判別式值為正數.'),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('請注意：迭代次數超過給定的最大值gxmax,請重新輸入初始值x1,x2,x3') return;endP=i,piancha,xdpiancha,xk,yk' 例 用拋物線法求方程e的一個實根的近似值，使精確到.解 先用作圖法確定初始值x01，x02和x03，然后在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk= paowu (-0.4,-

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