非線性方程組的數(shù)值解法的matlab程序

1、q第二章非線性方程（組）的數(shù)值解法本章主要介紹方程根的有關(guān)概念，求方程根的步驟，確定根的初始近似值的方法（作圖法，逐步搜索法等），求根的方法（二分法，迭代法，牛頓法，割線法，米勒（Mdler ）法和迭代法的加速等）及其 MATLA典序，求解非線性方程組的方法及其MATLA叁序.方程（組）的根及其 MATLABt令求解方程（組）的solve命令求方程f(x)=q(x)的根可以用MATLAB令: x=solve(方程f(x)=q(x),彳寺求符號變量x)求方程組fi(xi，,xn)=qi(xi，,xn) (i=1,2,n)的根可以用 MATLAB令:E1=sym(方程 f1(x1,xn)=q1(x

2、1,xn);En=sym(方程fn(x1，,xn)=qn(x1,xn);x1,x2, ,xn=solve(E1,E2,，En, x1,xn)求解多項式方程（組）的roots命令如果f（x）為多項式，則可分別用如下命令求方程f（x）=0的根，或求導數(shù) f （x）（見表 2-1 ）.表2-1求解多項式方程（組）的 roots 命令命令功能xk =roots(fa)輸入多項式f（x）的系數(shù)fa （按降哥排列），運行后輸 出xk為f （x） = 0的全部根.dfa=polyder(fa)輸入多項式f（x）的系數(shù)fa （按降哥排列），運行后輸 出df a為多項式f（x）的導數(shù)f（x）的系數(shù).dfx=po

3、ly2sym(dfa)輸入多項式f（x）的導致f （x）的系數(shù)dfa （按降哥排 一. . 1 .、 列），運行后輸出dfx為多項式f（x）的導數(shù)f （x）.求解方程（組）的fsolve命令如果非線性方程（組）是多項式形式，求這樣方程（組）的數(shù)值解可以直接調(diào)用上面已 經(jīng)介紹過的roots命令.如果非線性方程（組）是含有超越函數(shù)，則無法使用roots命令，需要調(diào)用MATLA廉統(tǒng)中提供的另一個程序 fsoke來求解.當然，程序fsoke也可以用于多項式方程 （組），但是它的計算量明顯比roots命令白勺大.fsolve命令使用最小二乘法（least squares method ）解非線性方程（組

4、）F（ X ） =0的數(shù)值解，其中X和F （X）可以是向量或矩陣.此種方法需要嘗試著輸入解 X的初始值（向 量或矩陣）X。，即使程序中的迭代序列收斂， 也不一定收斂到F（ X ） = 0的根（見例2.1.8 ）.fsolve 的調(diào)用格式：X=fsolve（F,X0）輸入函數(shù)F （x ）的“文件名和解X的初始值（向量或矩陣）X。，嘗試著解方程（組）F( X ) =0,運行后輸出F( X ) =0解的估計值(向量或矩陣)X.要了解更多的調(diào)用格式和功能請輸入：help fsolve，查看說明.搜索根的方法及其 MATLA能序求解非線性方程根的近似值時，首先需要判斷方程有沒有根？如果有根，有幾個根？如

5、果有根，需要搜索根所在的區(qū)間或確定根的初始近似值(簡稱初始值).搜索根的近似位置的常用方法有三種：作圖法、逐步搜索法和二分法等，使用這些方法的前提是高等數(shù)學中的 零點定理.作圖法及其MATLA醒序作函數(shù)的圖形的方法很多，如用計算機軟件的圖形功能畫圖，或用高等數(shù)學中應(yīng)用導數(shù)作圖，或用初等數(shù)學的函數(shù)疊加法作圖等.下面介紹兩種作圖程序.作函數(shù)y = f(x)在區(qū)間a,b的圖形的MATLA叁序一x=a:h:b; % h 是步長 y=f(x); plot(x,y) grid, gtext(y=f(x)說明： 此程序在MATLAB勺工作區(qū)輸入，運行后即可出現(xiàn)函數(shù)y=f(x)的圖形.此圖形與x軸交點的橫坐標

6、即為所要求的根的近似值. 區(qū)間a,b的兩個端點的距離 b-a和步長h的絕對值越小，圖形越精確.作函數(shù)y = f (x)在區(qū)間a,b上的圖形的 MATLA翼序二將y = f (x)化為h(x) = g(x),其中h(x)和g(x)是兩個相等的簡單函數(shù)x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2) grid,gtext( y1=h(x),y2=g(x) TOC o 1-5 h z 說明：此程序在MATLAB勺工作區(qū)輸入，運行后即可出現(xiàn)函數(shù) y1 = h(x)和y2 = g(x)的 圖形.兩圖形交點的橫坐標即為所要求的根的近似值下面舉例說明如何用計算機軟件

7、MATLABJ圖形功能作圖.逐步搜索法及其MATLA醒序逐步搜索法也稱試算法.它是求方程f (x) =0根的近似值位置的一種常用的方法.逐步搜索法依賴于尋找連續(xù)函數(shù)f(x)滿足f(a)與f(b)異號的區(qū)間a,b. 一旦找到區(qū)間，無論區(qū)間多大，通過某種方法總會找到一個根.MATLAB的庫函數(shù)中沒有逐步搜索法的程序，現(xiàn)提供根據(jù)逐步搜索法的計算步驟和它的 收斂判定準則編寫其主程序，命名為zhubuss.m .逐步搜索法的 MATLA在程序輸入?yún)^(qū)間端點a和b的值，步長h和精度tol ,運行后輸出迭代次數(shù)k=( b-a)/ h+1,方程f(x)=0根的近似值r.function k,r=zhubuss(

8、a,b,h,tol)%輸入的量-a和b是閉區(qū)間a,b的左、右端點；%h是步長；%-tol是預(yù)先給定的精度.%運行后輸出的量-k是搜索點的個數(shù)；% -r是方程 在a,b上的實根的近似值，其精度是tol;X=a:h:b;Y=funs(X);n=(b-a)/h+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n); for k=2:nX(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);endxielv=(Y(k+1)-Y(k)*(Y(k)-Y(k-1);if (abs(Y(k)t

9、ol)&( xielv k,r=zhubuss(-2,2,0.001,0.0001)運行后輸出的結(jié)果k=4001r = -1.2240-1.0000-1.0000-0.99901.2250即搜索點的個數(shù)為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.*x.A3+2.* x.V-3.* *-3用丫=$所(8$(2.* x.A3)代替，可得到方程在區(qū)間 -2,2上的根的近似值如下r = -1.9190-1.7640-1.5770-1.3300-0.

10、92200.92301.33101.57801.76501.9200如果讀者分別將方程的結(jié)果與例2.2.3比較，方程的結(jié)果與例2.1.2比較，將會發(fā)現(xiàn)逐步搜索法的MATLAB序的優(yōu)點.如果精度要求比較高，用這種逐步搜索法是不合算的二分法及其 MATLABg序二分法二分法也稱 逐次分半法.它的基本思想是：先確定方程 f(x) = 0含根的區(qū)間(a,b),再 把區(qū)間逐次二等分.我們可以根據(jù)式(2.3b)、區(qū)間a,b和誤差8 ,編寫二分法求方程根的迭代次數(shù)的MATLAB命令.已知閉區(qū)間a,b和誤差 名.用二分法求方程誤差不大于一的根的迭代次數(shù)k!TMATLA腌令k=-1+ceil(log(b-a)-

11、 log(abtol)/ log(2)% ceil 是向+ 方向取整，abtol 是誤差二分法的 MATLAB?序二分法需自行編制程序，現(xiàn)提供用二分法求萬程 f(x)=0的根x的近似值xk的步驟和式(2.3a)編寫一個名為 erfen.m的二分法的MATLAB主程序如下.二分法的MATLA在程序求解方程f(x)=0在開區(qū)間(a,b)內(nèi)的一個根的前提條件是f (x)在閉區(qū)間a,b上連續(xù)，且 f(a) f (b):二 0.輸入的量：a和b是閉區(qū)間a,b的左、右端點，abtol是預(yù)先給定的精度.運行后輸出的量：k是使用二分法的次數(shù).x是方程 在(a,b)內(nèi)的實根x的近似值，其精度是abtol.wuc

12、a=|bk-ak|/2是使用k次二分法所得到的小區(qū)間的長度的一半，即實根x的近似值 x的絕對精度限，滿足wucaw abtol .yx=f(xk),即方程 f(x)=0在實根 x*的近似值 x處的函數(shù)值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* yb0,disp(注意：ya*yb0,請重新調(diào)整區(qū)間端點 a和b.), return endmax1=-1+ceil(log(b-a)- 10g(abtol)/log(2);% ceil 是向十g

13、 方向取整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*yx0 b=x;yb=yx;else a=x; ya=yx;endif b-ax=-4:0.1:4;y=x.A3-x +4; plot(x,y) grid,gtext(y=xA3-x+4) 畫出函數(shù)f(x)=x3-x+4的圖像，如圖2-7.從圖一 、3像可以看出，此曲線有兩個駐點都在x3軸的上方，在(-2,-1)內(nèi)曲線與x軸

14、只有一個交點，則該方程有唯一一個實根，且在(-2,-1)內(nèi).方法2試算法.在MATLAB：作窗口輸入程序x=-4: 1:4,y=x.A3-x +4運行后輸出結(jié)果圖2-734y =由于連續(xù)函數(shù)(2)-56f( x)滿足-20-2f(-2) f(-1) k,x,wuca,yx=erfen (-2,-1,0.001)2-3 ,其余結(jié)果為運行后屏幕顯示用二分法計算過程被列入表k = 9 , x=-1.7959 ,wuca = 9.7656e-004 , yx = 0.0037迭代法及其MATLA程序確定根的近似位置以后，接下來的工作就是將根精確化，即按某種方法將初始近似值逐 TOC o 1-5 h z

15、 步精確化，直到滿足所要求的精確度為止 .上節(jié)介紹的二分法是將根精確化的方法之一，但是它的收斂速度較慢，且不能求出偶重根.迭代法可以克服這種缺陷.迭代法是求解方程的根、線性和非線性方程組的解的基本而重要的方法迭代法的MATLAB?序1迭代法需要自行編制程序.下面提供的迭代法的 MATLAB序1使用時只需輸入迭代初始 值X。、迭代次數(shù)k、迭代公式xk+i=中（xk）和一條命令，運行后就可以輸出求迭代序列xk、迭代k次得到的迭代值 xk和相鄰兩次迭代的偏差 piancha =| xk -xk_1| （簡稱偏差）和偏 差的相對誤差 xdpiancha =| xk - x /岡 的值，并且具有警報功能

16、（若迭代序列發(fā)散， 則提示用戶“請重新輸入新的迭代公式”；若迭代序列收斂，則屏幕會出現(xiàn)“祝賀您！此迭 代序列收斂，且收斂速度較快”）.我們可以用這個程序來判斷迭代序列的斂散性，也可以用于比較由一個方程得到的幾個迭代公式的斂散性的優(yōu)劣迭代法的MATLA典序1輸入的量：初始值x。、迭代次數(shù)k和迭代公式xk4=(xk)(k = 0,1,2,)；運行后輸出的量：迭代序列xJ、迭代k次得到的迭代值xk、相鄰兩次迭代的偏差 piancha =| xk xk|和它的偏差的相對誤差xdpiancha = |xk - xk1/xk|的值.根據(jù)迭代公式(2.4 )和已知條件，現(xiàn)提供名為_diedai1.m的M文件

17、如下function k,piancha,xdpiancha,xk=diedai1(x0,k)%輸入的量-x0是初始值，k是迭代次數(shù)x(1)=x0; for i=1:kx(i+1)=fun1(x(i); %程序中調(diào)用的 fun1.m 為函數(shù) y= 4 (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 xk end if (piancha 1)&(xdpiancha0.5)&(k3)disp(請用戶注意：此迭代序列發(fā)散，請重新輸入新的迭代公

18、式,)return ;endif (piancha 0.001)&(xdpiancha3)disp(祝賀您！此迭代序列收斂，且收斂速度較快,)return ; endp=(i-1) piancha xdpiancha xk;2例2.4.1 求萬程f(x) =x +2x-10的一個正根.解 在MATLAB：作窗口輸入程序 k,piancha,xdpiancha,xk= diedai1(2,5)運行后輸出用迭代公式xk中二(10 x2)/2的結(jié)果k,piancha,xdpiancha,xk=3.000000000000000.500000000000004.875000000000001.0000

19、00000000001.000000000000000.333333333333332.000000000000002.500000000000005.000000000000003.000000000000004.375000000000000.897435897435904.00000000000000 11.757812500000001.70828603859251-6.882812500000005.00000000000000 11.803741455078130.63167031671297 -18.68655395507813請用戶注意：此迭代序列發(fā)散，請重新輸入新的迭代公式k

20、=5 , piancha = 11.80374145507813, xdpiancha = 0.63167031671297,xk = -18.68655395507813由以上足仃后輸出的迭代序列與根x =2.316 624 790 355 40 相差越來越大，即迭代序列xk發(fā)散，此迭代法就失敗.這時偏差piancha逐漸增大且偏差的相對誤差xdpiancha的值大于0.5.用迭代公式xk書=10/(xk +2)運行后輸出的結(jié)果 k,piancha,xdpiancha,xk=1.000000000000002.000000000000003.000000000000004.000000000

21、000005.000000000000000.500000000000000.277777777777780.146198830409360.079264426125550.042304047651280.200000000000000.125000000000000.061728395061730.034626038781160.018144868631152.500000000000002.222222222222222.368421052631582.289156626506022.33146067415730k =5,piancha =0.04230404765128,xdpianch

22、a = 0.01814486863115,xk = 2.33146067415730可見，偏差piancha和偏差的偏差的相對誤差xdpiancha的值逐漸變小，且第5次的迭接近，則迭代序列xj收*代值 xk =2.331 460 674 157 30 與根 x =2.316 624 790 355 40斂，但收斂速度較慢，此迭代法較為成功用迭代公式xk書=xk -(xk2 +2xk 10)/(2xk +2)運行后輸出的結(jié)果k,piancha,xdpiancha,1.000000000000002.000000000000003.000000000000004.000000000000005.

23、00000000000000祝賀您！此迭代序列收斂,xk=0.333333333333330.016666666666670.000041876046900.000000000264370且收斂速度較快0.142857142857140.00719424460432 0.00001807631822 0.00000000011412 02.333333333333332.316666666666672.316624790619772.316624790355402.31662479035540k = 5; piancha = 0; xdpiancha = 0; y = 2.3166247903

24、5540.可見，偏差 piancha和偏差的相對誤差 xdpiancha的值越來越小，且第 5次的迭代值, ，一 *，、，一一 、 一xk=2.316 624 790 355 40 與根 x =2.316 624 790 355 40 相差無幾，則迭代序列 xk收斂,且收斂速度很快，此迭代法成功迭代法的MATLAB?序2迭代法的MATLA典序2輸入的量：初初始值小、精度tol和和迭代公式xk書=9(xk)(k = 0,1,2,)；運行后輸出的量：迭代次數(shù) k,滿足精度tol的根的近似根xk的迭代序列xk,相 鄰兩次迭代 的偏差pianch=| xk xk|和它的偏差的 相對誤差xdpianch

25、a =卜卜一xk /lx/的值，yk= 4 (xk).根據(jù)迭代公式(2.4 )和已知條件，現(xiàn)提供名為_diedai2.m 的M文件如下function k,piancha,xdpiancha,xk,yk=diedai2(x0,tol,ddmax) x(1)=x0;for i=1: ddmaxx(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 (pianchatol)|(xdpi

26、anchaddmaxdisp(迭代次數(shù)超過給定的最大值ddmax)return ;endP=(i-1),piancha,xdpiancha,xk,yk;例2.4.5求x5-3x+1=0在0.3附近的根，精確到4位小數(shù).解構(gòu)造迭代公式 xk =中(xk ).改寫原方程為等價方程x = (x5+1)/3.這時中(x) =(x5+1)/3,初始值為xo=0.5,迭代公xk1 =(x5 1)/3 (k=0,1,2,).利用迭代法的MATLAB序2計算精確到4位小數(shù)的根的近似值y和迭代次數(shù)k.在MATLAB：作窗口輸入程序 k,piancha,xdpiancha,xk,yk=diedai2(0.3,1e

27、-4,100)運行后輸出的結(jié)果k,piancha,xdpiancha,xk,yk=00.034140.102180.300000.3341410.034140.102180.334140.3347220.000580.001730.334720.3347330.000010.000040.334730.33473k =3;piancha =1.206089525390697e-005;xdpiancha =3.603129477781680e-005;xk =0.3347; yk =0.3347.2.5 迭代過程的加速方法及其MATLA能序?qū)τ谑諗康牡^程，只要迭代足夠多次，就可以使結(jié)果達到

28、任意的精度，但有時迭代過程收斂緩慢，從而使計算量變得很大，因此，迭代過程的加速是個重要的問題.2.5.2 加權(quán)迭代法的MATLA抑序加權(quán)迭代法的MATLA在程序已知初始值乂0、精度tol和迭代公式xk+=中(xk),求滿足精度tol的根的近似根xk和 它的函數(shù)值yk及迭代次數(shù)k.根據(jù)加權(quán)迭代的公式(2. 10)和已知條件，現(xiàn)提供名為jasudd.m的M文件如下：function k,xk,yk=jasudd (x0,tol,L,ddmax)x1(1)=x0;for i=1: ddmaxx(i+1)=fun(x1(i);x1(i+1)= x(i+1)+L*( x(i+1)- x1(i)/(1-L

29、);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 ykif (pianchatol)|(xdpianchaddmaxdisp(迭代次數(shù)超過給定的最大值ddmax)return ;endP=(i-1),xk,yk;例2.5.1 用(2.10)式求2ex=0在0.85附近的一個近似根，要求精度8=106.解 在MATLAB：作窗口輸入程序 k,xk,yk=jasudd (0.85,1e-6,-0.855,100)運行后輸出結(jié)果k,xk,yk=1

30、.000000000000002.000000000000003.00000000000000k =3 ; xk =0.852606 ;0.852603700280410.852605499754910.85260550207699yk = 0.852606.0.852607038305610.852605504062360.852605502082552.5.4艾特肯(Aitken)加速方法的MATLAB?序艾特肯加速方法的MATLA典序已知初始值x0、精度tol和迭代公式xk+仁中(xk),求滿足精度tol的根的近似根xk和 它的函數(shù)值 yk及迭代次數(shù) k, p= k, x1, x2, x

31、.根據(jù)艾特肯加速方法的公式(2.11)和已知條件，現(xiàn)提供名為Aitken.m 的M文件如下:function k,xk,yk,p= Aitken (x0,tol, ddmax)x(1)=x0;for i=1: ddmaxx1(i+1)=fun(x(i); x2(i+1)=fun(x1(i+1);x(i+1)=x2(i+1)-(x2(i+1)-x1(i+1)42/(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 (p

32、ianchatol)|(xdpianchaddmaxdisp(迭代次數(shù)超過給定的最大值ddmax)return ;endm=0,1:i-1; p=m,x1,x2,x;例2.5.3用艾特肯加速方法求2e=-x = 0在0.85附近的一個近似根，要求精度；=10*.解 在MATLAB：作窗口輸入程序 k,xk,yk,p= Aitken(0.85,1e-6, 100)運行后輸出結(jié)果k=3, xk=0.85260550201343 p =01.00000000000000 2.00000000000000 3.00000000000000,yk=0.8526055020137300.854829863

33、897450.852604364918110.8526055020134300.850711106524840.852606471508260.852605502013980.850000000000000.852606835686070.852605502014070.852605502013732.6 牛頓(Newton)切線法及其MATLAgg序牛頓切線法的收斂性及其MATLA勰序判別牛頓切線法的局部收斂性的MATLABE程序根據(jù)牛頓切線法局部收斂的條件f (x)2，.人、，- r-* . -_輸入的量：x是萬程f(x)=0根x或滿足1% - x* 1的初始值x0,函數(shù)f(x)及 其一階

34、導數(shù)f(x)和二階導數(shù)f(x)；輸出的量：W(x或()的值.如果憎“ 1或*(%) 1時，運行后輸出信息 恭喜您！此迭彳序列收斂 ，4 (x)的導數(shù)值的絕對值 y=|d 4 (x)/dx| 和方程f(x)=0 的函 數(shù)f(x)值f如下；否則，運行后輸出信息請注意觀察下面顯示的4 (x)的導數(shù)值的絕對值y=|d 4 (x)/dx| 和方程f(x)=0 的函數(shù)f(x) 值.根據(jù)牛頓切線法局部收斂的條件和方程f(x)=0的函數(shù)f(x)及其一、二階導數(shù)，現(xiàn)提 供程序:function y,f=newjushou(x)f=fnq(x); fz=fnq(x)*ddfnq(x)/(dfnq(x)A2+eps

35、);y=abs(fz);if (y y,f=newjushou(-1)運行后輸出結(jié)集請注意觀察下面顯示的 小(x)的導數(shù)值的絕對值 y=|d小(x)/dx|和方程f(x)=0的函數(shù)f(x)值y =139.5644 f =4.3096.說明：實際上，這個初始值所產(chǎn)生的迭代序列xk 發(fā)散，見例2.6.6.(2)輸入程序 y,f=newjushou(0)運行后輸出結(jié)果請注意觀察下面顯示的 小(x)的導數(shù)值的絕對值 y=|d小(x)/dx|和方程f(x)=0的函數(shù)f(x)值y =8.0000f =4.,說明：實際上，這個初始值所產(chǎn)生的迭代序列xk收斂到方程的根2.925 32 ,而不是收斂到離它最近的

36、根1.400 81,即不是局部收斂的，見例 2.6.6.(3)輸入程序 y,f=newjushou(1)運行后輸出結(jié)果恭喜您！此迭代序列收斂，(j)(x)導數(shù)值的絕對值y=|d(j)(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下 y =0.3566f =1.7126說明：實際上，這個初始值所產(chǎn)生的迭代序列xk收斂到離它最近的根 1.400 81,即xk是局部收斂的，見例 2.6.6.(4)輸入程序： y,f=newjushou(2)運行后輸出結(jié)果請注意觀察下面顯示的小(x)的導數(shù)值的絕對值 y=d小(x)/dx|和方程f(x)=0的函數(shù)f(x)值y =1.2593 f =-2.7188

37、說明：雖然y =1.259 31,不滿足牛頓切線法局部收斂的條件 p(x*)1 .但是，實 際上，這個初始值所產(chǎn)生的迭代序列xk 收斂到離它最近的根 1.400 81,即xk是局部收斂的.這說明世() y,f=newjushou(5.5)運行后輸出結(jié)果請注意觀察下面顯示的小(x)的導數(shù)值的絕對值 y=|d小(x)/dx|和方程f(x)=0的函數(shù)f(x)值y =1.0447e+005 f =176.6400.說明：實際上，這個初始值所產(chǎn)生的迭代序列&k收斂到235.619,但不是此方程的根，比較初始值5.5與迭代值235.619 ,可見不是局部收斂的，見例 2.6.6.(6)輸入程序 y,f=n

38、ewjushou(8)運行后輸出結(jié)果恭喜您！此迭代序列收斂，小(x)導數(shù)值的絕對值y=|d小(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下y =0.4038f =-2.9452e+003說明：實際上，這個初始值所產(chǎn)生的迭代序列xk收斂到方程的根6.290 60 ,而不是收斂到離它最近的根 9.424 46,即不是局部收斂的，見例 2.6.6.牛頓切線法的 MATLAB?序牛頓切線法求方程f(x)=0的近似根xk (精度為tol )需要自行編制程序.牛頓切線法的 MATLA在程序輸入的量：初始值x。、近似根xk的誤差限tol ,近似根xk的函數(shù)值f(xk)的誤差限 f tol ,迭代次數(shù)

39、的最大值 gxmax函數(shù)fnq (x) = f(x)及其導數(shù)dfnq (x) = f (x).輸出的量：迭代序列xk、迭代k次得到的迭代值 xk (迭代次數(shù)超過gxmax時，運 行后輸出信息 請注意：迭代次數(shù)超過給定的最大值 gxmax)、相鄰兩次迭代的偏差 piancha = |xk - xk1和它的偏差的相對誤差 xdpiancha = |xk - xk_11/l xk |的值.根據(jù)式(2. 12)和已知條件，現(xiàn)提供名為 newtonqx.m 的M文件：function k,xk,yk,piancha,xdpiancha=newtonqx(x0,tol,ftol,gxmax) x(1)=x

40、0;for i=1: gxmaxx(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+1;xk=x(i);yk=fnq(x(i); (i-1) xk yk piancha xdpianchaif (abs(yk)ftol)&(pianchatol)|(xdpianchagxmaxdisp(請注意：迭代次數(shù)超過給定的最大值gxmax。)return ;end(i-1),xk,yk,piancha,xdpiancha;例2.6.3用牛頓切線法求方

41、程 2x3 -3x2 +1=0在x0 = -0.4和0.9的近似根，要求精3度名=10 .然后判斷此方法分別在方程的根x =-0.5和1處的收斂速度.解 (1)先求方程的近似根.在MATLAB：作窗口輸入程序 k,xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001,0.001,100)運行后輸出初始值Xo =.4的迭代結(jié)果如表2-10所示.表 2-10kXkykXk -Xk_1xk -Xk/Ixk1-0.516 7-0.076 70.116 70.225 82-0.500 4-0.001 60.016 30.032 63P -0.500 0-0.000 00

42、.000 40.000 7迭代次數(shù)k = 3,根的近似值Xk =-0.500，它的函數(shù)值yk =-1.759e-013,偏差piancha=1.712e-007 和相對誤差 xdpiancha =3.423e-007.如果將步驟命令中的-0.4改作0.9,則運行后輸出初始值x0 = -0.4的迭代結(jié)果為k =7; piancha =7.290e-004; xdpiancha =7.295e-004; xk = 0.999; yk =1.590e-006.(2)再討論收斂速度比較初始值分別是 x0 = -0.4和0.9的兩組結(jié)果，在 x0 = -0.4處迭代3次就得到單根* 一一一 一 * .一

43、 一 x =-0.5的近似值 Xk =-0.500;而在X0 =0.9處迭代7次才得到二重根 x = 1的近似值 xk=0.999.可見，牛頓切線迭代法在單根處的迭代速度比二重根處的迭代速度快很多.這正如前面討論的結(jié)果一樣，即若 x是f(x)=0的單根，則牛頓切線法是平方收斂的；若 x是 f (x) = 0的二重根，則牛頓切線法是線性收斂的；我們也可以用階的定義判斷它們的收斂速度.留給讀者證明.求VC的方法及其MATLAB?序用c的n次根點的迭代公式求 我的近似值xk (精度為tol )需要自行編制程序.求c的n次根n/C (當n是偶數(shù)時，c0)的matlaBE程序輸入的量：初始值 x。、根指

44、數(shù)n、被開方數(shù)c、nG的誤差限tol ,迭代次數(shù)的最大 值 gxmax輸出的量：迭代序列xk、迭代k次得到的迭代值Xk (迭代次數(shù)達到最大值gxmax時，運行后輸出信息,請注意：迭代次數(shù)超過給定的最大值gxmax)、相鄰兩次迭代的偏差piancha =| Xk - Xk 1| 和它的偏差的相對誤差xdpiancha = |xk - xk|/|xk| 的值.根據(jù)求C的n次根C的迭代公式(2.25)和已知條件，現(xiàn)提供名為kai2fang .m的M文件: function k,xk,yk,piancha,xdpiancha,P=kainfang(x0,c,n,tol, gxmax)x(1)=x0;f

45、or i=1: gxmaxu(i)= (x(i)An-c)/(n*x(i)A(n-1); 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,xdpianchaif (pianchatol)|(xdpianchagxmaxdisp(請注意：迭代次數(shù)超過給定的最大值gxmax.)(i-1),xk,yk,piancha,xdpianchareturn ;endP=(i-1),xk,yk,piancha,

46、xdpiancha;例2.6.5求Ji13,要求精度為10-6.解本題介紹四種解法.方法1用求c的n次卞nc (當n是偶數(shù)時， O0)的MATLA翼序at算.取初始值x0=10,根指數(shù)n=2,被開方數(shù)c=113,近似根的精度tol =10,，迭代 的最大次數(shù)gxmax=100.在工作區(qū)間輸入程序 k,xk,yk,piancha,xdpiancha=kainfang(10,113,2,1e-5,100)運行后輸出結(jié)果k=4,piancha=1.610800381968147e-011,xdpiancha=1.515313534117706e-012xk =10.63014581273465,yk

47、 =1.710910332060509e+009可見，、HT3定10.630 15 ,滿足精度10”.方法2 用牛頓迭代公式(2.12 )計算.設(shè) &=113,貝Ux2 113 = 0,記 f (x) =x2 -113 , f(x) = 2x.由牛頓迭代公式得，xk書=xk -k，即 xM = ：(xk + )(其中k =0,1,2,3,)2xk2 xk取初始值x0 =10,計算結(jié)果列入表 2-12.因為，迭代次數(shù)k =4時，偏差 x4 -x3 k,xk,yk,piancha,xdpiancha=newtonqx(10, 1e-5, 1e-5,100)運行后，將輸出的結(jié)果列入下表2-13.迭代

48、k =4次，得到精度為103的結(jié)果J而定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：作空間輸入程序 113A(1/2)運行后輸出ans = 10.63014581273465經(jīng)過四舍五入后，得到精度為10”的結(jié)果J113=10.630 15.2

49、.6.6牛頓切線法的加速及其兩種MATLA醒序. 當x是方程f(x) =0的m (m之2)重根時，由牛頓切線公式(2.12)產(chǎn)生的迭代序 列xk 收斂到x*的速度較慢，是線性收斂的.我們要提高收斂速度，可以有選擇地采用已 知根的重數(shù)和未知根的重數(shù)兩種求重根的修正牛頓切線公式(一) 已知方程根的重數(shù) m已知方程f (x) = 0根的重數(shù) m ,求重根的修正牛頓切線法的MATLABfc程序輸入的量：初始值 小，方程 “*)=0根乂的重數(shù)m,近似根xk的誤差限toi ,近似 根xk的函數(shù)值f(xk)誤差限ftol ,迭代次數(shù)的最大值 gxmax函數(shù)fnq(x)= f(x)及其導數(shù) dfnq (x)

50、= f (x);輸出的量：迭代序列xj、迭代k次得到的迭代值xk (迭代次數(shù)達到最大值 gxmax時運 行后輸出信息 請注意：迭代次數(shù)超過給定的最大值gxmax)、相鄰兩次迭代的偏差 piancha=| xk xk|和它的偏差的相對誤差 xdpiancha = |xk xkJ/|xk|的值.根據(jù)式(2.26)和已知條件，現(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

51、(i)+eps);piancha=abs(x(i+1)-x(i);xdpiancha=piancha/( abs(x(i+1)+eps);i=i+1;xk=x(i);yk=fnq(x(i);if (pianchatol)|(xdpiancha tol)&(abs(yk)gxmaxdisp(請注意：迭代次數(shù)超過給定的最大值 gxmax.) return ; endP=(i-1),piancha,xdpiancha,xk,yk;例2.6.7 判斷x* =0是方程f (x) =2e3x 9x2 6x 2 = 0的幾重根？在區(qū)間0,1上，分別用牛頓切線法和求重根的修正牛頓切線公式求此根的近似值xk ,

52、使精確到;=10”.一 ，一 一 ，一. * 一 . . . . . .解 經(jīng)判斷知，x =0是方程f(x)=0的三重根.在MATLAB作窗口輸入程序 k,piancha,xdpiancha,xk,yk= newtonqx (0.5,1e-4, 1e-4,100)運行后整理結(jié)果列入表2-20.表 2-20pianchaxdpianchaxkyk10.144 100.404 890.355 900.541 9820.107 410.432 250.248 490.168 2030.077 450.452 800.171 040.051 4640.054 500.467 600.116 550.0

53、15 58 J50.037 690.477 980.078 850.004 6960.025 760.485 140.053 100.001 40 ：70.017 460.490 010.035 630.000 4280.011 770.493 300.023 860.000 12190.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 660.004 750.000 00130.001 580.499 110.003

54、 170.000 00140.001 050.499 400.002 110.000 00150.000 700.499 580.001 410.000 00160.000 470.499 690.000 940.000 00170.000 310.499 730.000 630.000 00180.000 210.499 670.000 420.000 00190.000 140.499 440.000 280.000 00200.000 090.498 860.000 190.000 00迭代次數(shù)k=20,精確到6 = 10的根的近似值是xk =0.000 2,其函數(shù)值是 yk =5.75

55、58e-011 ,是一階(線性)收斂速度 .根據(jù)求重根的修正牛頓切線公式，在MATLAB：作窗口輸入程序 k,piancha,xdpiancha,xk,yk= newtonxz(3,0.5,1e-4, 1e-4,100)運行后整理結(jié)果得表 2-21.表 2-21pianchaxdpianchaxkyk10.432 306.385 790.067 700.002 9420.066 5457.310 770.001 160.000 0030.001 162673.312 010.000 00-0.000 0040.000 130.996 700.000 130.000 0050.000 13151

56、.535 850.000 00-0.000 0060.000 100.991 440.000 100.000 0070.000 1088.974 840.000 00-8.881 78e-016迭代次數(shù) k=7,精確到君=104的根的近似值是 xk =1.121 1e-006, piancha = 9.975 3e-005, xdpiancha=88.974 8,其函數(shù)值是 yk = -8.881 8e-016,是二階收斂(平方收斂)速度 . 可見，求重根的修正牛頓切線公式比牛頓切線法收斂速度快得多(二)未知方程根的重數(shù)未知方程f(x)=0根的重數(shù)為 m ,求重根的修正牛頓切線法的MATLA在

57、程序輸入的量：初始值 x0,近似根Xk的誤差限tol ,近似根Xk的函數(shù)值f(xk)誤差限 ftol ，迭代次數(shù)的最大值 gxmax函數(shù)fnq(x)= f(x)及其一、二階導數(shù)dfnq (x) = f(x) 和 ddfnq (x) = f (x).運行后輸出的量：迭代序列x、迭代k次得到的迭代值xk (迭代次數(shù)達到最大 值gxmax時，運行后輸出信息請注意：迭代次數(shù)超過給定的最大值gxmax)、相鄰兩次迭代的偏差 piancha = %xy網(wǎng)它的偏差的相對誤差xdpiancha = |xk 乂/區(qū)|的值.根據(jù)式(2.27)和已知條件，現(xiàn)提供名為newtonxz1.m 的M文件function

58、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)A2+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(i);yk=fnq(x(i); if (pianchatol)|(xdpiancha tol)&(abs(yk)

59、gxmax disp(請注意：迭代次數(shù)超過給定的最大值gxmax.)return ; endP=(i-1),piancha,xdpiancha,xk,yk;例2.6.8用未知重數(shù)的求重根的修正牛頓切線法，求方程在區(qū)間0,1上的根，使精確到6=10”，并且與例2.6.7的結(jié)果比較.解 在MATLAB作窗口輸入程序 k,piancha,xdpiancha,xk,yk=newtonxz1(0.5,1e-4, 1e-4,100)運行后整理結(jié)果得表2-22.表 2-22kpianchaxdpianchaxkyk10.599 236.038 57-0.099 23-0.008 1820.096 9843.

60、054 79-0.002 25-0.000 0030.002 251590.394 68-0.000 000.000 0040.000 020.945 62-0.000 03-0.000 00即迭代 k =4 次，piancha =2.461 55e-005 , xdpiancha =0.945 62 , xk =-2.603 09e-005 , yk =-1.325 61e-013.這些結(jié)果與例 2.6.7 的結(jié)果比較見下表 2-23.表 2-23牛頓切線法已知重數(shù)的修正牛頓切線法未知重數(shù)的修正牛頓切線法迭代次數(shù)k2074根的近似值xk1.858 1e-0041.121 1e-006-2.6