本章主要介紹方程根的有關(guān)概念

1、第二章 非線性方程（組）的數(shù)值解法本章主要介紹方程根的有關(guān)概念，求方程根的步驟，確定根的初始近似值的方法（作圖法，逐步搜索法等），求根的方法（二分法，迭代法，牛頓法，割線法，米勒（Müller）法和迭代法的加速等）及其MATLAB程序，求解非線性方程組的方法及其MATLAB程序.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命令

2、:>>E1=sym('方程f1(x1,xn)=q1(x1,xn)');.En=sym('方程fn(x1,xn)=qn(x1,xn)');x1,x2,xn=solve(E1,E2,En, x1,xn)求解多項(xiàng)式方程（組）的roots命令如果為多項(xiàng)式，則可分別用如下命令求方程的根，或求導(dǎo)數(shù)（見表 2-1）.表 2-1 求解多項(xiàng)式方程（組）的roots命令命 令功 能xk =roots(fa)輸入多項(xiàng)式的系數(shù)fa（按降冪排列），運(yùn)行后輸出xk為的全部根.dfa=polyder(fa)輸入多項(xiàng)式的系數(shù)fa（按降冪排列），運(yùn)行后輸出dfa為多項(xiàng)式的導(dǎo)數(shù)的系數(shù).

3、dfx=poly2sym(dfa)輸入多項(xiàng)式的導(dǎo)數(shù)的系數(shù)dfa（按降冪排列），運(yùn)行后輸出dfx為多項(xiàng)式的導(dǎo)數(shù). 求解方程（組）的fsolve命令如果非線性方程（組）是多項(xiàng)式形式，求這樣方程（組）的數(shù)值解可以直接調(diào)用上面已經(jīng)介紹過(guò)的roots命令.如果非線性方程（組）是含有超越函數(shù)，則無(wú)法使用roots命令，需要調(diào)用MATLAB系統(tǒng)中提供的另一個(gè)程序fsolve來(lái)求解.當(dāng)然，程序fsolve也可以用于多項(xiàng)式方程（組），但是它的計(jì)算量明顯比roots命令的大. fsolve命令使用最小二乘法（least squares method）解非線性方程（組）的數(shù)值解，其中X和F（X）可以是向量或矩陣.此

4、種方法需要嘗試著輸入解X的初始值（向量或矩陣）X0，即使程序中的迭代序列收斂，也不一定收斂到的根（見例2.1.8）.fsolve的調(diào)用格式： X=fsolve(F,X0)輸入函數(shù)的M文件名和解X的初始值（向量或矩陣）X0，嘗試著解方程（組），運(yùn)行后輸出解的估計(jì)值（向量或矩陣）X.要了解更多的調(diào)用格式和功能請(qǐng)輸入:help fsolve,查看說(shuō)明.2.2 搜索根的方法及其MATLAB程序求解非線性方程根的近似值時(shí)，首先需要判斷方程有沒有根？如果有根，有幾個(gè)根？如果有根，需要搜索根所在的區(qū)間或確定根的初始近似值（簡(jiǎn)稱初始值）.搜索根的近似位置的常用方法有三種:作圖法、逐步搜索法和二分法等，使用這些

5、方法的前提是高等數(shù)學(xué)中的零點(diǎn)定理. 作圖法及其MATLAB程序作函數(shù)的圖形的方法很多，如用計(jì)算機(jī)軟件的圖形功能畫圖，或用高等數(shù)學(xué)中應(yīng)用導(dǎo)數(shù)作圖，或用初等數(shù)學(xué)的函數(shù)疊加法作圖等.下面介紹兩種作圖程序.作函數(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程

6、序二將化為，其中是兩個(gè)相等的簡(jiǎn)單函數(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)即為所要求的根的近似值.下面舉例說(shuō)明如何用計(jì)算機(jī)軟件MATLAB的圖形功能作圖. 逐步搜索法及其MATLAB程序逐步搜索法也稱試算法.它是求方程根的近似值位置的一種常用的方法. 逐步搜索法依賴于尋找連續(xù)函數(shù)滿足與異號(hào)的區(qū)間.一旦找到區(qū)間，無(wú)論區(qū)間多大，通過(guò)某種方法總會(huì)找到一個(gè)根.MATLAB的庫(kù)函數(shù)中沒有逐步搜索法

7、的程序，現(xiàn)提供根據(jù)逐步搜索法的計(jì)算步驟和它的收斂判定準(zhǔn)則編寫其主程序，命名為zhubuss.m.逐步搜索法的MATLAB主程序 輸入?yún)^(qū)間端點(diǎn)a和b的值，步長(zhǎng)和精度tol，運(yùn)行后輸出迭代次數(shù) k=(b-a)/h+1,方程 根的近似值r.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:b;Y=funs(X);n=(b-a)/h+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n)

8、;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.解 在MATLAB工作窗口輸入如下程序>> k,r=zhubuss(-2,2,0.001,0.

9、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.*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.9200如果讀

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

11、B程序二分法需自行編制程序，現(xiàn)提供用二分法求方程f(x)=0的根的近似值的步驟和式（2.3a）編寫一個(gè)名為erfen.m的二分法的MATLAB主程序如下.二分法的MATLAB主程序 求解方程在開區(qū)間（a,b）內(nèi)的一個(gè)根的前提條件是在閉區(qū)間a,b上連續(xù)，且.輸入的量：a和b是閉區(qū)間a,b的左、右端點(diǎn)，abtol是預(yù)先給定的精度.運(yùn)行后輸出的量：k是使用二分法的次數(shù).x是方程 在（a,b）內(nèi)的實(shí)根x*的近似值，其精度是abtol.wuca=|bk-ak|/2是使用k次二分法所得到的小區(qū)間的長(zhǎng)度的一半，即實(shí)根x*的近似值x的絕對(duì)精度限,滿足wuca abtol.yx=f(xk) ，即方程f(x)=0

12、在實(shí)根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* yb>0, disp('注意：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);

13、 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, endendk=max1; x; wuca; yx=fun(x);例2.3.3 確定方程x3-x+4=0的實(shí)根的分布情況，并用二分法求在開區(qū)間 (-2,-1)內(nèi)的實(shí)根的近似值，要求精度為0.001.圖2-7解 （1）先用兩種方法確定方程x3-x+4=0 的實(shí)根的分布情況。方法1 作圖法.在MATLAB工作窗口輸入如下程序>&g

14、t;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的圖像，如圖2-7.從圖像可以看出，此曲線有兩個(gè)駐點(diǎn)都在x軸的上方，在(-2,-1)內(nèi)曲線與x軸只有一個(gè)交點(diǎn)，則該方程有唯一一個(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

15、)內(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，yx = 0.00372.4 迭代法及其MATLAB程序確定根的近似位置以后，接下來(lái)的工作就是將根精確化，即按某種方法將初始近似值逐步精確化，直到滿足所要求的精確度為止. 上節(jié)介紹的二分法是將根精確化的方法之一，但是它的迭代法的MATLAB程序1迭代法需要自行編制程序.下面提供的迭代法的MATLAB程序1使用時(shí)只

16、需輸入迭代初始值、迭代次數(shù)k、迭代公式xk+1=(xk)和一條命令，運(yùn)行后就可以輸出求迭代序列、迭代k次得到的迭代值和相鄰兩次迭代的偏差piancha =| （簡(jiǎn)稱偏差）和偏差的相對(duì)誤差xdpiancha=|的值，并且具有警報(bào)功能（若迭代序列發(fā)散，則提示用戶“請(qǐng)重新輸入新的迭代公式”；若迭代序列收斂，則屏幕會(huì)出現(xiàn)“祝賀您！此迭代序列收斂，且收斂速度較快”）.我們可以用這個(gè)程序來(lái)判斷迭代序列的斂散性，也可以用于比較由一個(gè)方程得到的幾個(gè)迭代公式的斂散性的優(yōu)劣.迭代法的MATLAB程序1 輸入的量：初始值、迭代次數(shù)k和迭代公式；運(yùn)行后輸出的量：迭代序列、迭代k次得到的迭代值、相鄰兩次迭代的偏差pia

17、ncha =|和它的偏差的相對(duì)誤差xdpiancha=的值.根據(jù)迭代公式（2.4）和已知條件，現(xiàn)提供名為diedai1.m的M文件如下function 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

18、(piancha >1)&(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,x

19、dpiancha,xk= diedai1(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.7

20、0828603859251 -6.88281250000000請(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.

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

22、624 790 355 40相差無(wú)幾，則迭代序列收斂，且收斂速度很快，此迭代法成功.迭代法的MATLAB程序2迭代法的MATLAB程序2 輸入的量：初初始值、精度tol和和迭代公式；運(yùn)行后輸出的量：迭代次數(shù)k，滿足精度tol的根的近似根的迭代序列，相鄰兩次迭代的偏差pianch=|和它的偏差的相對(duì)誤差xdpiancha =的值，yk=(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: ddmax x(i+1)=fun(

23、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;endend if i>ddmax disp('迭代次數(shù)超過(guò)給定的最大值ddmax') return;endP=(i-1),piancha,xdpiancha,xk,yk' 求x5-3x+1=0在

24、0.3附近的根，精確到4位小數(shù).解 構(gòu)造迭代公式.改寫原方程為等價(jià)方程.這時(shí)，初始值為x0=0.5, 迭代公式   .利用迭代法的MATLAB程序2計(jì)算精確到4位小數(shù)的根的近似值 y 和迭代次數(shù)k.在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk=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

25、.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程序 對(duì)于收斂的迭代過(guò)程，只要迭代足夠多次，就可以使結(jié)果達(dá)到任意的精度，但有時(shí)迭代過(guò)程收斂緩慢，從而使計(jì)算量變得很大，因此，迭代過(guò)程的加速是個(gè)重要的問(wèn)題.MATLAB程序加權(quán)迭代法的MATLAB主程序已知初始值、精度tol和迭代公式xk+1=(xk)，求滿足精度tol的根的近似根和它的函數(shù)值yk及迭代次數(shù)k.根據(jù)加權(quán)迭代的公式（2.10）和

26、已知條件，現(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-1;

27、xk=x1(i); return;endendif i>ddmax disp('迭代次數(shù)超過(guò)給定的最大值ddmax') return;endP=(i-1),xk,yk' 用（2.10）式求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.00000000000000 0.85260549975491 0.85

28、260550406236 3.00000000000000 0.85260550207699 0.85260550208255k =3； xk =0.852606； yk = 0.852606.MATLAB程序艾特肯加速方法的MATLAB程序已知初始值、精度tol和迭代公式xk+1=(xk)，求滿足精度tol的根的近似根和它的函數(shù)值yk及迭代次數(shù)k，p=k',x1',x2',x'.根據(jù)艾特肯加速方法的公式（2.11）和已知條件，現(xiàn)提供名為Aitken.m的M文件如下：function k,xk,yk,p= Aitken (x0,tol, ddmax)x(1)=x

29、0;for i=1: ddmax x1(i+1)=fun(x(i); x2(i+1)=fun(x1(i+1); x(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=x1(i); return;endendif i>ddmax disp('迭

30、代次數(shù)超過(guò)給定的最大值ddmax') 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.85000000000000 1.00000000000000 0.85482986389745 0.85071110652484 0.85260

31、683568607 2.00000000000000 0.85260436491811 0.85260647150826 0.85260550201407 3.00000000000000 0.85260550201343 0.85260550201398 0.852605502013732.6 牛頓(Newton)切線法及其MATLAB程序牛頓切線法的收斂性及其MATLAB程序判別牛頓切線法的局部收斂性的MATLAB主程序根據(jù)牛頓切線法局部收斂的條件 輸入的量：x是方程根或滿足的初始值，函數(shù)及其一階導(dǎo)數(shù)和二階導(dǎo)數(shù)；輸出的量：或的值.如果或時(shí)，運(yùn)行后輸出信息'恭喜您！此迭代序列收斂,(

32、x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下'；否則，運(yùn)行后輸出信息'請(qǐng)注意觀察下面顯示的 (x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值'.根據(jù)牛頓切線法局部收斂的條件和方程的函數(shù)及其一、二階導(dǎo)數(shù)，現(xiàn)提供程序：function y,f=newjushou(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ù)

33、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工作窗口輸入程序>> 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說(shuō)明：實(shí)際上，這個(gè)初始值所產(chǎn)生的迭代序列發(fā)散,見（2）輸入程序 >

34、;> y,f=newjushou(0)運(yùn)行后輸出結(jié)果請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值y = 8.0000f = 4說(shuō)明：實(shí)際上，這個(gè)初始值所產(chǎn)生的迭代序列收斂到方程的根2.925 32，而不是收斂到離它最近的根1.400 81,即不是局部收斂的，見（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.3566f = 1.7126說(shuō)明：實(shí)際上，這個(gè)初始值所產(chǎn)生的迭代序列收斂到離它最近的

35、根1.400 81,即是局部收斂的，見（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說(shuō)明：雖然y =1.259 3，不滿足牛頓切線法局部收斂的條件 .但是，實(shí)際上，這個(gè)初始值所產(chǎn)生的迭代序列收斂到離它最近的根1.400 81,即是局部收斂的.這說(shuō)明是充分條件，非必要條件（見 （5）輸入程序 >> y,f=newjushou(5.5)運(yùn)行后輸出結(jié)果請(qǐng)注意觀察下面顯示的(x)的導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和

36、方程f(x)=0的函數(shù)f(x)值y = 1.0447e+005f = 176.6400說(shuō)明：實(shí)際上，這個(gè)初始值所產(chǎn)生的迭代序列收斂到235.619,但不是此方程的根，比較初始值5.5與迭代值235.619，可見不是局部收斂的，見（6）輸入程序 >> y,f=newjushou(8)運(yùn)行后輸出結(jié)果恭喜您！此迭代序列收斂,(x)導(dǎo)數(shù)值的絕對(duì)值y=|d(x)/dx|和方程f(x)=0的函數(shù)f(x)值f如下y = 0.4038f = -2.9452e+003說(shuō)明：實(shí)際上，這個(gè)初始值所產(chǎn)生的迭代序列收斂到方程的根6.290 60，而不是收斂到離它最近的根9.424 46,即不是局部收斂的，見

37、牛頓切線法的MATLAB程序牛頓切線法求方程的近似根（精度為tol）需要自行編制程序.牛頓切線法的MATLAB主程序輸入的量：初始值、近似根的誤差限tol，近似根的函數(shù)值的誤差限ftol，迭代次數(shù)的最大值gxmax、函數(shù)fnq(x)及其導(dǎo)數(shù)dfnq(x).輸出的量：迭代序列、迭代k次得到的迭代值 (迭代次數(shù)超過(guò)gxmax時(shí)，運(yùn)行后輸出信息'請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax')、相鄰兩次迭代的偏差piancha =和它的偏差的相對(duì)誤差xdpiancha=的值.根據(jù)式（2.12）和已知條件，現(xiàn)提供名為newtonqx.m的M文件：function k,xk,yk,pian

38、cha,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+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);(

39、i-1) xk yk piancha xdpiancha return;endend if i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax。') return;end (i-1),xk,yk,piancha,xdpiancha' 用牛頓切線法求方程在的近似根，要求精度.然后判斷此方法分別在方程的根處的收斂速度.解 （1）先求方程的近似根.在MATLAB工作窗口輸入程序>> k,xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001, 0.001,100)運(yùn)行后輸出初始值的迭代結(jié)果如表 2-10所示

40、.表 2-10kxkyk1-0.516 7-0.076 70.116 70.225 82-0.500 4-0.001 60.016 30.032 63-0.500 0-0.000 00.000 40.000 7迭代次數(shù)k = 3,根的近似值xk =-0.500,它的函數(shù)值yk =-1.759e-013,偏差piancha =1.712e-007和相對(duì)誤差xdpiancha =3.423e-007.如果將步驟的迭代結(jié)果為k =7；piancha =7.290e-004；xdpiancha =7.295e-004；xk = 0.999；yk =1.590e-006. （2）再討論收斂速度比較初始值

41、分別是的兩組結(jié)果，在處迭代3次就得到單根的近似值xk =-0.500；而在處迭代7次才得到二重根的近似值xk單根處的迭代速度比二重根處的迭代速度快很多.這正如前面討論的結(jié)果一樣，即若是的單根，則牛頓切線法是平方收斂的；若是的二重根，則牛頓切線法是線性收斂的；我們也可以用階的定義判斷它們的收斂速度.留給讀者證明.求的方法及其MATLAB程序用的次根的迭代公式求的近似值（精度為tol）需要自行編制程序.求的次根（當(dāng)n是偶數(shù)時(shí)，）的MATLAB主程序輸入的量：初始值、根指數(shù)、被開方數(shù)、的誤差限tol，迭代次數(shù)的最大值gxmax.輸出的量：迭代序列、迭代k次得到的迭代值 (迭代次數(shù)達(dá)到最大值gxmax

42、時(shí)，運(yùn)行后輸出信息'請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax')、相鄰兩次迭代的偏差piancha =|和它的偏差的相對(duì)誤差xdpiancha=的值.根據(jù)求的次根的迭代公式(2.25)和已知條件，現(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); x(i+1)= x(i)-u(i); piancha=abs(x(i+1)-x(i); xdpiancha=pia

43、ncha/( 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); return;endendif i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax.')(i-1),xk,yk,piancha,xdpiancha return;endP=(i-1),xk,yk,piancha,xdpiancha' 求，要求精度為.解 本題介紹四種解法.方法1

44、用求的次根（當(dāng)n是偶數(shù)時(shí)，）的MATLAB程序計(jì)算.取初始值，根指數(shù)，被開方數(shù)，近似根的精度tol=，迭代的最大次數(shù)gxmax=100.在工作區(qū)間輸入程序>> k,xk,yk,piancha,xdpiancha=kainfang(10,113,2,1e-5,100)運(yùn)行后輸出結(jié)果可見,10.630 15，滿足精度.方法2 用牛頓迭代公式（2.12）計(jì)算.設(shè)，則，記，.由牛頓迭代公式得，即取初始值，計(jì)算結(jié)果列入表 2-12. 表 2-12迭代次數(shù)偏差根的近似值10.650 00010.650 00020.019 83610.630 16430.000 01910.630 14640.

45、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=newtonqx(10, 1e-5, 1e-5,100)運(yùn)行后，將輸出的結(jié)果列入下表 2-13. 迭代k=4次，得到精度為的結(jié)果 10.630 15.表 2-13kpianchaxdpianchaxkyk1

46、0.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ò)四舍五入后，得到精度為的結(jié)果10.630 15.牛頓切線法的加速及其兩種MATLAB程序當(dāng)是方程的重根時(shí)，由牛頓切線公式（2.12）產(chǎn)生的迭代序列收斂到的速度較慢，是線性收斂的.我們要提高收斂速度，可以有選擇地

47、采用已知根的重?cái)?shù)和未知根的重?cái)?shù)兩種求重根的修正牛頓切線公式.（一） 已知方程根的重?cái)?shù)已知方程根的重?cái)?shù),求重根的修正牛頓切線法的MATLAB主程序輸入的量：初始值，方程根的重?cái)?shù)，近似根的誤差限tol，近似根的函數(shù)值誤差限ftol，迭代次數(shù)的最大值gxmax、函數(shù)fnq(x)及其導(dǎo)數(shù)dfnq(x)；輸出的量：迭代序列、迭代k次得到的迭代值(迭代次數(shù)達(dá)到最大值gxmax時(shí)運(yùn)行后輸出信息'請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax')、相鄰兩次迭代的偏差piancha =|和它的偏差的相對(duì)誤差xdpiancha=的值. 根據(jù)式（2.26）和已知條件，現(xiàn)提供名為newtonxz.m的M文

48、件：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=x(i);yk=fnq(x(i); if (piancha<tol)|(xdpiancha< tol)&(abs(yk)<ftol) k=i-1; xk=x(i);(i-1)

49、piancha xdpiancha xk yk return;endend if i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過(guò)給定的最大值gxmax.') return;endP=(i-1),piancha,xdpiancha,xk,yk' 判斷是方程e的幾重根？在區(qū)間上，分別用牛頓切線法和求重根的修正牛頓切線公式求此根的近似值，使精確到.解 經(jīng)判斷知，是方程的三重根.在MATLAB工作窗口輸入程序>> k,piancha,xdpiancha,xk,yk= newtonqx (0.5,1e-4, 1e-4,100)運(yùn)行后整理結(jié)果列入表 2-20.表

50、 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.054 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

51、.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 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.0

52、00 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，精確到的根的近似值是xk =0.000 2，其函數(shù)值是yk =5.755 8e-011，是一階（線性）收斂速度.根據(jù)求重根的修正牛頓切線公式，在MATLAB工作窗口輸入程序 >> k,piancha,xdpiancha,xk,yk= newtonxz(3,0.5,1e-4, 1e-4,100)運(yùn)行后整理結(jié)果得表 2-21.表 2-21kpianchaxdpianchaxky

53、k10.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.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，精確到的根的近似值是xk =1.121 1e-006，piancha = 9.975 3e-005，xdpiancha=88.974 8，其函數(shù)值是yk = -8.881 8e-016，是二階收斂（平方收斂）速度.可見，求重根的修正牛頓切線公式比牛頓切線法收斂速度快得多.（二） 未知方程根的重?cái)?shù)未知方程根的重?cái)?shù)為，求重根的修正牛頓切線法的MATLAB主程序輸入的量：初始值，近似根的誤差限to