代碼資源精通matlab科學計算10-3rd

精通MATLAB科學計算（第2版）精通MATLAB科學計算（第2版）第10第10章非線性方程求解PAGE254PAGE255第10章非線性方程求解非線性方程是常見的一類方程，非線性方程（組）的理論遠不如線性方程（組）成熟和有效，特別是非線性方程組解的存在唯一性還沒有完全解決，判斷其解的存在性和解的個數(shù)幾乎沒有可行的方法。利用MATLAB提供的函數(shù)，可以求解一些簡單的非線性方程；通過編程，可以解決一些較為復雜的非線性方程。通過本章學習，讀者能夠熟練掌握MATLAB中的非線性方程求解相關的函數(shù)，而且能通過編程實現(xiàn)多種求解非線性方程的數(shù)值算法。10.1MATLAB中非線性方程求根函數(shù)MATLAB中求解非線性方程的函數(shù)主要是以下兩個：fzero和fsolve。其中fzero可用來求解非線性方程，fsolve可用來求解非線性方程和非線性方程組。10.1.1fzero函數(shù)fzero函數(shù)用于求解單變量非線性方程。它的，其基本用法及功能如下：1．x=fzero(f,x0)求函數(shù)f在x0附近的根；2．x=fzero(f,x0,options)使用優(yōu)化選項，options可選的值很多，具體可參考MATLAB幫助；3．[x,fval]=fzero(f,x0,options)返回值中的fval=；4．[x,fval,exitflag]=fzero(f,x0,options)其中exitflag標志著求解狀態(tài)。函數(shù)fzero的具體用法見下面的例子?！纠?0-1】非線性方程求解函數(shù)fzero的應用實例。采用fzero函數(shù)求方程x2+x1=0在附近的根。解：在MATLAB命令窗口中輸入下列命令：>>x=fzero('x^2+x-1',0.5)輸出計算結果為：x=0.6180從結果可知，方程在附近的根為。函數(shù)fzero還可以求函數(shù)在某個區(qū)間上的一個根，例如，求取方程在[0,2]區(qū)間上的根，繼續(xù)在MATLAB命令窗口中輸入下列命令：>>x=fzero('sin(x)-0.5',[02])輸出計算結果為：x=0.5236從結果可知，方程在[0,2]區(qū)間上的根為。需要注意的是，當用fzero求函數(shù)在某個區(qū)間上的一個根時，此函數(shù)在區(qū)間兩端點處的值符號要相反才行。函數(shù)fzero還有一些特殊的用法，例如，求取函數(shù)當時在附近的根，在MATLAB命令窗口中輸入下列命令：>>u=inline('x^y-3','x','y');%建立非線性方程>>c=fzero(u,2,[],2)輸出計算結果為：c=1.7321從結果可知，所求的根為。優(yōu)化選項options的使用方法如下，使用optimset函數(shù)，其中參數(shù)TolX選項用來設置求解精度，在MATLAB命令窗口中輸入下列命令：>>opt=optimset('TolX',1.0e-8);>>x=fzero('sin(x)-0.5',0.5,opt)x=0.5236>>opt=optimset('TolX',1.0e-2);>>x=fzero('sin(x)-0.5',0.5,opt)x=0.5283從運算結果可以看出，不同精度對根有影響。10.1.2fsolve函數(shù)fsolve函數(shù)用于求解多變量非線性方程組。它的基本用法及功能如下：1．x=fsolve(fun,x0)用fun定義向量函數(shù)，其定義方式可以參考例子；2．x=fsolve(fun,x0,options)options為優(yōu)化選項；3．[x,val]=fsolve(…)val=F(x)，即函數(shù)值向量；4．[x,val,exitflag]=fsolve(…)exitflag標志著求解狀態(tài)；5．[x,val,exitflag,output]=fsolve(…)output包含著優(yōu)化后的結果信息；6．[x,val,exitflag,output,jacobian]=fsolve(…)jacobian為解x處的Jacobian陣。其余參數(shù)與前面參數(shù)相似。另外，符號數(shù)學工具箱中也有一個可用于求解方程組的函數(shù)solve，它可求各種類型方程組的解析解。當找不到解析解時，solve會自動尋求一個近似解，并且精度很高。需要注意的是，函數(shù)fsolve不是MATLAB符號工具箱中的函數(shù)，它位于優(yōu)化工具箱內，而solve是一個符號函數(shù)?！纠?0-2】非線性方程組求解函數(shù)fsolve應用實例1。采用fsolve函數(shù)求方程組在附近的根。解：首先新建一個.m文件，命名為myfun.m，鍵入下列內容：functionF=myfun(x)F=[x(1)^2-x(2)-2;x(2)^2-2*x(1)-4];在MATLAB命令窗口中輸入下列命令：>>x0=[22];>>opt=optimset('Display','iter');>>[r,val]=fsolve(@myfun,x0,opt)輸出計算結果為：NormofFirst-orderTrust-regionIterationFunc-countf(x)stepoptimalityradius0316161160.69900414.621290.0001181560.1239650.03562.53124.3305e-0110.003045211.96e-0052.54156.63221e-0241.8626e-0066.78e-0122.5Optimizationterminated:first-orderoptimalityislessthanoptions.TolFun.r=2.21432.9032val=1.0e-011*0.20750.1525由計算結果可知，方程組的一個根為?！纠?0-3】非線性方程組求解函數(shù)fsolve應用實例2。采用fsolve函數(shù)求矩陣方程的一個根，初始解向量為x0=[1,1;1,1]。解：首先新建一個.m文件，命名為myfun.m，鍵入下列內容：functionF=matrixfun(x)F=x*x-[-11,-8;24,-3];在MATLAB命令窗口中輸入下列命令：>>x0=[1,1;1,1];>>opt=optimset('Display','iter');>>r=fsolve(@matrixfun,x0)輸出計算結果為：Optimizationterminated:first-orderoptimalityislessthanoptions.TolFun.r=1.0000-2.00006.00003.0000可以驗證矩陣的平方剛好等于。10.2其他數(shù)值求根法下面講述的幾種數(shù)值方法，在不加說明的情況下，都是在某個區(qū)間內求得方程的一個根。其理論依據(jù)為：如果，則在區(qū)間上方程至少存在一個實根。大體上可以把求解非線性方程的方法分為夾逼法（二分法和黃金分割法）和迭代法（其他方法）兩類，其中迭代法的內容十分豐富，包括一些加速和優(yōu)化的技巧。下面進行講述。10.2.1二分法二分法的具體求解步驟如下。（1）計算函數(shù)在區(qū)間中點的函數(shù)值，并作進行下面的判斷：如果，轉到（2）；如果，令，轉到（1）；如果，則為一個根。（2）如果（為預先給定的精度），則為一個根，否則令，轉到（1）。在MATLAB中編程實現(xiàn)的二分法的函數(shù)為：HalfInterval。功能：用二分法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=HalfInterval(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。二分法的MATLAB程序代碼如下：functionroot=HalfInterval(f,a,b,eps)%二分法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a); %兩端點的函數(shù)值f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elseroot=FindRoots(f,a,b,eps); %調用求解子程序endfunctionr=FindRoots(f,a,b,eps)f_1=subs(sym(f),findsym(sym(f)),a);f_2=subs(sym(f),findsym(sym(f)),b);mf=subs(sym(f),findsym(sym(f)),(a+b)/2); %中點函數(shù)值if(f_1*mf>0)t=(a+b)/2;r=FindRoots(f,t,b,eps); %右遞歸elseif(f_1*mf==0)r=(a+b)/2;elseif(abs(b-a)<=eps)r=(b+3*a)/4; %輸出根elses=(a+b)/2;r=FindRoots(f,a,s,eps); %左遞歸endendend【例10-4】二分法求解非線性方程應用實例。采用二分法求方程在區(qū)間[0,1]上的一個根。解：在MATLAB命令窗口中輸入：>>r=HalfInterval('x^3-3*x+1',0,1)輸出計算結果為：r=0.34730.3473可見，方程在區(qū)間[0,1]上的一個根為。10.2.2黃金分割法二分法是把求解區(qū)間的長度逐次減半，而黃金分割法是把求解區(qū)間逐次縮短為前次的0.618倍。它的求解步驟如下。（1）設，，且；（2）如果(給定的最小區(qū)間長度)，則輸出方程的根為；否則轉到（3）；（3）如果，則令，轉（1），否則如果，令，反之令，轉到（1）。在MATLAB中編程實現(xiàn)的黃金分割法的函數(shù)為：hj。功能：用黃金分割法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=hj(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。黃金分割法的MATLAB程序代碼如下：functionroot=hj(f,a,b,eps)%黃金分割法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elset1=a+(b-a)*0.382;t2=a+(b-a)*0.618;f_1=subs(sym(f),findsym(sym(f)),t1);f_2=subs(sym(f),findsym(sym(f)),t2);tol=abs(t1-t2);while(tol>eps)%精度控制if(f_1*f_2<0)a=t1;%區(qū)間兩端點都調整b=t2;elsefa=subs(sym(f),findsym(sym(f)),a);if(f_1*fa>0)a=t2;%同號左端點調整elseb=t1;%異號右端點調整endendt1=a+(b-a)*0.382;t2=a+(b-a)*0.618;f_1=subs(sym(f),findsym(sym(f)),t1);f_2=subs(sym(f),findsym(sym(f)),t2);tol=abs(t2-t1);endroot=(t1+t2)/2;%輸出根end【例10-5】黃金分割法求解非線性方程應用實例。采用黃金分割法求方程在區(qū)間[0,1]上的一個根。解：在MATLAB命令窗口中輸入：>>r=hj('x^3-3*x+1',0,1)輸出計算結果為：r=0.34730.3473可見，方程在區(qū)間[0,1]上的一個根為，與二分法求出來的結果相同。10.2.3不動點迭代法求方程的根，可以先把方程改寫成如下的形式：于是得到不動點迭代法的其中一種迭代公式：此種不動點迭代法很有可能不收斂，因為它的本質是求函數(shù)與直線的交點，而它們不一定存在交點。即使收斂，其速度也十分慢，因此有了艾特肯加速迭代與史蒂芬森加速迭代。在MATLAB中編程實現(xiàn)的不動點迭代法的函數(shù)為：StablePoint。功能：用不動點迭代法求函數(shù)的一個零點。調用格式：[root,n]=StablePoint(f,x0,eps)。其中，f為函數(shù)名；x0為初始迭代向量；eps為根的精度；root為求出的函數(shù)零點；n為迭代步數(shù)。不動點迭代法的MATLAB程序代碼如下：function[root,n]=StablePoint(f,x0,eps)%用不動點迭代法求函數(shù)的一個零點%初始迭代向量：x0%根的精度：eps%求出的函數(shù)零點：root%迭代步數(shù)：nif(nargin==2)eps=1.0e-4;endtol=1;root=x0;n=0;while(tol>eps)n=n+1;r1=root;root=subs(sym(f),findsym(sym(f)),r1)+r1;%迭代的核心公式tol=abs(root-r1);end【例10-6】不動點迭代法求解非線性方程應用實例。采用不動點迭代法求方程的一個根。解：在MATLAB命令窗口中輸入：>>[r,n]=StablePoint('1/sqrt(x)+x-2',0.5)r=0.38200.3820n=44從計算結果可以看出，經(jīng)過4步迭代，得出方程的一個根為。1．艾特肯加速迭代法艾特肯加速迭代法是在計算出后，對作以下修正：然后用來逼近方程的根。艾特肯加速迭代法是先用不動點迭代法算出一系列，再對此系列作修正得到系列，用后者來逼近方程的根。在MATLAB中編程實現(xiàn)的艾特肯加速迭代法的函數(shù)為：AtkenStablePoint。功能：用艾特肯加速迭代法求函數(shù)的一個零點。調用格式：[root,n]=AtkenStablePoint(f,x0,eps)。其中，f為函數(shù)名；x0為初始迭代向量；eps為根的精度；root為求出的函數(shù)零點；n為迭代步數(shù)。艾特肯加速迭代法的MATLAB程序代碼如下：function[root,n]=AtkenStablePoint(f,x0,eps)%用艾特肯加速迭代法求函數(shù)f的一個零點%初始迭代向量：x0%根的精度：eps%求出的函數(shù)零點：root%迭代步數(shù)：nif(nargin==2)eps=1.0e-4;endtol=1;root=x0;x(1:2)=0;n=0;m=0;a2=x0;while(tol>eps)n=n+1;a1=a2;r1=root;root=subs(sym(f),findsym(sym(f)),r1)+r1; %求出根x(n)=root; %把所有的根都保存下來if(n>2)m=m+1;a2=x(m)-(x(m+1)-x(m))^2/(x(m+2)-2*x(m+1)+x(m)); %對根進行艾特肯修正tol=abs(a2-a1);endendroot=a2;【例10-7】艾特肯加速迭代法求解非線性方程應用實例。采用艾特肯加速迭代法求方程的一個根，迭代初始值為0.5。解：在MATLAB命令窗口中輸入：>>[r,n]=AtkenStablePoint('1/sqrt(x)+x-2',0.5)r=0.38200.3820n=44從計算結果可以看出，經(jīng)過4步迭代，得出方程的一個根為。這個方程比較簡單，所以通過簡單的幾步迭代就可得出根。從下面的結果可以看出艾特肯加速迭代法的優(yōu)點。>>[r,n]=AtkenStablePoint('1/sqrt(x)+x-2',0.999)r=1.00001.0000n=44>>[r,n]=StablePoint('1/sqrt(x)+x-2',0.999)r=0.38200.3820n=2121上面的例子說明采用初始值0.999時，經(jīng)過4步艾特肯加速迭代法，可以得到方程的另一個根，而普通的不動點迭代法卻得不到。2．史蒂芬森加速迭代法史蒂芬森加速與艾特肯加速不同的地方在于前者是在迭代的同時就進行修正，它的迭代公式為：在MATLAB中編程實現(xiàn)的史蒂芬森加速迭代法的函數(shù)為：StevenStablePoint。功能：用史蒂芬森加速迭代法求函數(shù)的一個零點。調用格式：[root,n]=StevenStablePoint(f,x0,eps)。其中，f為函數(shù)名；x0為初始迭代向量；eps為根的精度；root為求出的函數(shù)零點；n為迭代步數(shù)。史蒂芬森加速的MATLAB程序代碼如下：function[root,n]=StevenStablePoint(f,x0,eps)%用史蒂芬森加速迭代法求函數(shù)f的一個零點%初始迭代向量：x0%根的精度：eps%求出的函數(shù)零點：root%迭代步數(shù)：nif(nargin==2)eps=1.0e-4;endtol=1;root=x0;n=0;while(tol>eps)n=n+1;r1=root;y=subs(sym(f),findsym(sym(f)),r1)+r1;z=subs(sym(f),findsym(sym(f)),y)+y;root=r1-(y-r1)^2/(z-2*y+r1); %對每次算出的根立即進行修正tol=abs(root-r1);end【例10-8】史蒂芬森加速迭代法求解非線性方程應用實例。采用史蒂芬森加速迭代法求方程的一個根。解：在MATLAB命令窗口中輸入：>>[r,n]=StevenStablePoint('1/sqrt(x)+x-2',0.999)r=1.00001.0000n=22>>[r,n]=StevenStablePoint('1/sqrt(x)+x-2',0.5)r=0.38200.3820n=44對比上例可以看出，史蒂芬森加速迭代法不僅能求出方程的一個根，而且它比艾特肯加速迭代法更快地求出另一個根。10.2.4弦截法弦截法的算法過程如下：（1）過兩點，作一直線，它與軸有一個交點，記為；（2）如果，過兩點，作一直線，它與軸的交點記為，否則過兩點，作一直線，它與軸的交點記為；（3）如此下去，直到，就可以認為為在區(qū)間上的一個根。（4）的遞推公式為：且。在MATLAB中編程實現(xiàn)的弦截法的函數(shù)為：Secant。功能：用弦截法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=Secant(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。弦截法的MATLAB程序代碼如下：functionroot=Secant(f,a,b,eps)%弦截法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);root=a-(b-a)*fa/(fb-fa);%迭代初始值while(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);s=fx*fa;if(s==0)root=r1;elseif(s>0)root=b-(r1-b)*fb/(fx-fb);elseroot=a-(r1-a)*fa/(fx-fa);endendtol=abs(root-r1);endend【例10-9】弦截法求解非線性方程應用實例。采用弦截法求方程在區(qū)間[1,4]上的一個根。解：在MATLAB命令窗口中輸入：>>r=Secant('log(x)+sqrt(x)-2',1,4)輸出計算結果為：r=1.87731.8773由計算結果可知，在區(qū)間[1,4]上的一個根為。10.2.5史蒂芬森弦截法史蒂芬森弦截法是弦截法的一種變形，它的遞推公式為：且有。史蒂芬森弦截法比弦截法的精度要高，收斂速度要快。在MATLAB中編程實現(xiàn)的史蒂芬森弦截法的函數(shù)為：StevenSecant。功能：用史蒂芬森弦截法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=StevenSecant(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。史蒂芬森弦截法的MATLAB程序代碼如下：functionroot=StevenSecant(f,a,b,eps)%史蒂芬森弦截法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);faa=subs(sym(f),findsym(sym(f)),fa+a);root=a-fa*fa/(faa-fa);%迭代初始值while(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);v=fx+r1;fxx=subs(sym(f),findsym(sym(f)),v);root=r1-fx*fx/(fxx-fx);%遞推的核心公式tol=abs(root-r1);endend【例10-10】史蒂芬森弦截法求解非線性方程應用實例。采用史蒂芬森弦截法求方程在區(qū)間[1,4]上的一個根。解：在MATLAB命令窗口中輸入：>>r=StevenSecant('log(x)+sqrt(x)-2',1,4)輸出計算結果為：r=11結果出現(xiàn)錯誤，改正方法可以是把求解區(qū)間縮短，縮短為[1.5,2]：>>r=StevenSecant('log(x)+sqrt(x)-2',1.5,2)r=1.87731.8773這樣就得出了正確的根1.8773。從上面的例子可以看出求解區(qū)間的選擇也是很重要的，一般來說，求解區(qū)間越短越好。10.2.6拋物線法弦截法其實是用不斷縮短的直線來近似函數(shù)，而拋物線法則采用三個點來近似函數(shù)，并且用拋物線與橫軸的交點來逼近函數(shù)的根。它的算法過程如下：（1）選定初始值，并計算和以下差分：一般取，，。注意不要使三點共線。（2）用牛頓插值法對三點進行插值得到一條拋物線，它有兩個根：其中兩個根中只取靠近的那個根，即號取與同號，即（3）用代替，重復以上步驟，并有以下遞推公式：其中（4）進行精度控制。在MATLAB中編程實現(xiàn)的拋物線法的函數(shù)為：Parabola。功能：用拋物線法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=Parabola(f,a,b,x,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；x為初始迭代點；eps為根的精度；root為求出的函數(shù)零點。拋物線法的MATLAB程序代碼如下：functionroot=Parabola(f,a,b,x,eps)%拋物線法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%初始迭代點：x%根的精度：eps%求出的函數(shù)零點：rootif(nargin==4)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);fx=subs(sym(f),findsym(sym(f)),x);d1=(fb-fa)/(b-a);d2=(fx-fb)/(x-b);d3=(d2-d1)/(x-a);B=d2+d3*(x-b);root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3));t=zeros(3);t(1)=a;t(2)=b;t(3)=x;while(tol>eps)t(1)=t(2); %保存三個點t(2)=t(3);t(3)=root;f1=subs(sym(f),findsym(sym(f)),t(1)); %計算三點的函數(shù)值f2=subs(sym(f),findsym(sym(f)),t(2));f3=subs(sym(f),findsym(sym(f)),t(3));d1=(f2-f1)/(t(2)-t(1)); %計算三個差分d2=(f3-f2)/(t(3)-t(2));d3=(d2-d1)/(t(3)-t(1));B=d2+d3*(t(3)-t(2)); %計算算法中的Broot=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3));tol=abs(root-t(3));endend【例10-11】拋物線法求解非線性方程應用實例1。采用拋物線法求方程在區(qū)間[1,4]上的一個根，迭代初始值為2。解：在MATLAB命令窗口中輸入：>>r=Parabola('sqrt(x)+log(x)-2',1,4,2)r=1.87731.8773【例10-12】拋物線法求解非線性方程應用實例2。分別從兩個區(qū)間[0.5,1.5]和[0.1,1]上采用拋物線法求方程的根。解：在MATLAB命令窗口中輸入：>r=11>>r=Parabola('1/sqrt(x)+x-2',0.1,1,0.5) %迭代區(qū)間為[0.1,1]，初始值為0.5r=0.38200.3820由此可以看出拋物線法也能求出方程的兩個根。10.2.7牛頓法弦截法本質上是一種割線法，它從兩端向中間逐漸逼近方程的根；牛頓法本質上是一種切線法，它從一端向一個方向逼近方程的根，其遞推公式為：初始值可以取和的較大者，這樣可以加快收斂速度。和牛頓法有關的還有簡化牛頓法和牛頓下山法。在MATLAB中編程實現(xiàn)的牛頓法的函數(shù)為：NewtonRoot。功能：用牛頓法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=NewtonRoot(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。牛頓法的MATLAB程序代碼如下：functionroot=NewtonRoot(f,a,b,eps)%牛頓法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fun=diff(sym(f)); %求導數(shù)fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb) %初始值取兩端點導數(shù)較大者root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1); %求該點的導數(shù)值root=r1-fx/dfx; %迭代的核心公式tol=abs(root-r1);endend【例10-13】牛頓法求解非線性方程應用實例。采用牛頓法求方程在區(qū)間[0.5,2]上的一個根。解：在MATLAB命令窗口中輸入：>>r=NewtonRoot('sqrt(x)-x^3+2',0.5,2)輸出計算結果為：r=1.47591.4759由計算結果可知，的一個根為。需要注意的是，初始值的選擇不要使得其導數(shù)為0。1．簡化牛頓法將牛頓法的遞推公式改為：為保證收斂，系數(shù)只須滿足。如果取常數(shù)，則稱為簡化牛頓法。在MATLAB中編程實現(xiàn)的簡化牛頓法的函數(shù)為：SimpleNewton。功能：用簡化牛頓法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=SimpleNewton(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。簡化牛頓法的MATLAB代碼如下：functionroot=SimpleNewton(f,a,b,eps)%簡化牛頓法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fun=diff(sym(f));%求導數(shù)fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb)%初始值取兩端點導數(shù)較大者df0=1/dfa;root=a-df0*fa;elsedf0=1/dfb;root=b-df0*fb;endwhile(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);root=r1-df0*fx;%迭代的核心公式，其中df0為常數(shù)tol=abs(root-r1);endend【例10-14】簡化牛頓法求解非線性方程應用實例。采用簡化牛頓法求方程在區(qū)間[1.2,2]上的一個根。解：在MATLAB命令窗口中輸入：>>r=SimpleNewton('sqrt(x)-x^3+2',1.2,2)輸出計算結果為：r=1.47591.4759由計算結果可知，的一個根為。2．牛頓下山法牛頓下山法收斂速度很快，很容易就可得出方程的根。牛頓下山法的遞推公式為：稱為下山因子，它的確定方法如下：（1）先取，求出；（2）判斷條件是否滿足；不滿足的話令；（4）不斷驗證直到。在MATLAB中編程實現(xiàn)的牛頓下山法的函數(shù)為：NewtonDown。功能：用牛頓下山法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=NewtonDown(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)零點。牛頓下山法的MATLAB程序代碼如下：functionroot=NewtonDown(f,a,b,eps)%牛頓下山法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fun=diff(sym(f));fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb) %迭代初始值的選取root=a;elseroot=b;endwhile(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1);toldf=1;alpha=2;whiletoldf>0 %此循環(huán)為下山因子的選取過程alpha=alpha/2;root=r1-alpha*fx/dfx; %迭代的核心公式fv=subs(sym(f),findsym(sym(f)),root);toldf=abs(fv)-abs(fx);endtol=abs(root-r1);endend【例10-15】牛頓法下山求解非線性方程應用實例。采用牛頓下山法求方程在區(qū)間[1.2,2]上的一個根。解：在MATLAB命令窗口中輸入：>>r=NewtonDown('sqrt(x)-x^3+2',1.2,2)輸出計算結果為：r=1.47591.4759由計算結果可知，的一個根為。10.2.8兩步迭代法兩步迭代法的一般公式為：它比較常用的有以下兩種形式：（1）（2）在MATLAB中編程實現(xiàn)的兩步迭代法的函數(shù)為：TwoStep。功能：用兩步迭代法求函數(shù)在某個區(qū)間上的一個零點。調用格式：root=TwoStep(f,a,b,type,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；type為兩步迭代法的形式；eps為根的精度；root為求出的函數(shù)零點。兩步迭代法的MATLAB代碼如下：functionroot=TwoStep(f,a,b,type,eps)%兩步迭代法求函數(shù)f在區(qū)間[a,b]上的一個零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%兩步迭代法的形式：type（可取1，2）%根的精度：eps%求出的函數(shù)零點：rootif(nargin==4)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fun=diff(sym(f));fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb)root=a;elseroot=b;endwhile(tol>eps)if(type==1) %第一種形式的兩步迭代法r1=root;fx1=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1);r2=r1-fx1/dfx; %第一步迭代fx2=subs(sym(f),findsym(sym(f)),r2);root=r2-fx2/dfx; %第二步迭代tol=abs(root-r1);else %第二種形式的兩步迭代法r1=root;fx1=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1);r2=r1-fx1/dfx; %第一步迭代fx2=subs(sym(f),findsym(sym(f)),r2);root=r2-fx2*(r2-r1)/(2*fx2-fx1); %第二步迭代tol=abs(root-r1);endendend【例10-16】兩步迭代法求解非線性方程應用實例。利用兩步迭代法的兩種形式求方程在區(qū)間[0.1,3]上的一個根。解：在MATLAB命令窗口中輸入：>>r=TwoStep('(log(x)-sin(x)+1)',0.1,3,1)%采用第一種形式的兩步迭代法r=0.70130.7013>>r=TwoStep('(log(x)-sin(x)+1)',0.1,3,2)%采用第二種形式的兩步迭代法r=0.70130.7013從上例可以看出，兩種形式的兩步迭代法得出的結果是一樣的。10.2.9重根迭代法如果方程存在重根，那么前面的方法都可能不收斂或者收斂速度很慢，因此針對二重根的情況有以下的迭代公式：但它要求是二階導數(shù)，如果方程比較復雜，計算量則比較大。在MATLAB中編程實現(xiàn)的重根的迭代法的函數(shù)為：MultiRoot。功能：用迭代法求函數(shù)在某個區(qū)間上的一個多重零點。調用格式：root=MultiRoot(f,a,b,eps)。其中，f為函數(shù)名；a為區(qū)間左端點；b為區(qū)間右端點；eps為根的精度；root為求出的函數(shù)多重零點。重根的迭代法的MATLAB代碼如下：functionroot=MultiRoot(f,a,b,eps)%求函數(shù)f在區(qū)間[a,b]上的一個多重零點%函數(shù)名：f%區(qū)間左端點；a%區(qū)間右端點：b%根的精度：eps%求出的函數(shù)多重零點：rootif(nargin==3)eps=1.0e-4;endf1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if(f1==0)root=a;endif(f2==0)root=b;endif(f1*f2>0)disp('兩端點函數(shù)值乘積大于0!');return;elsetol=1;fun=diff(sym(f)); %一階導數(shù)ddf=diff(sym(fun)); %二階導數(shù)fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb)root=a;elseroot=b;endwhile(tol>eps)r1=root;fx=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1); %一階導數(shù)值ddfx=subs(sym(ddf),findsym(sym(ddf)),r1); %二階導數(shù)值root=r1-fx*dfx/(dfx*dfx-fx*ddfx); %迭代公式tol=abs(root-r1);endend【例10-17】求解非線性方程綜合應用實例。用重根迭代法、兩步迭代法和弦截法求方程在區(qū)間[-2,3]上的一個根。解：在MATLAB命令窗口中輸入：>>r=MultiRoot('(sin(x)-x+2)*x*x',-2,3,1.0e-8)%重根迭代法r=00>>r=TwoStep('(sin(x)-x+2)*x*x',-2,3,1.0e-8)%兩步迭代法r=-1.0147e-005-1.0147e-005>>r=Secant('(sin(x)-x+2)*x*x',-2,3,1.0e-8)%弦截法r=2.55422.5542方程有一個二重根0，從計算結果可以看出，用其他方法得出的精度都不好，而用重根迭代法可以得到精確的結果。10.3非線性方程組的數(shù)值解法非線性方程組的數(shù)值解法大多來源于非線性方程的解法，只不過數(shù)的計算換成了矩陣的計算。10.3.1不動點迭代法非線性方程組的不動點迭代法與非線性方程的不動點迭代法基本一樣。此處公式中的變量皆為向量。在MATLAB中編程實現(xiàn)的非線性方程組的不動點迭代法的函數(shù)為：mulStablePoint。功能：用不動點迭代法求非線性方程組的一個解。調用格式：[r,n]=mulStablePoint(x0,eps)。其中，x0為初始迭代向量；eps為迭代精度；r為求出的解向量；n為迭代步數(shù)。不動點迭代法的MATLAB代碼如下：function[r,n]=mulStablePoint(x0,eps)%不動點迭代法求非線性方程組的一個解%初始迭代向量：x0%迭代精度：eps%解向量：r%迭代步數(shù)：nifnargin==1eps=1.0e-4;endr=myf(x0);n=1;tol=1;whiletol>epsx0=r;r=myf(x0)+x0; %迭代公式tol=norm(r-x0); %n=n+1;if(n>100000) %迭代步數(shù)控制disp('迭代步數(shù)太多，可能不收斂！');return;endend【例10-18】不動點迭代法求解非線性方程組應用實例。采用不動點迭代法求下面方程組的一個根。初始迭代值?。?,0）。解：首先建立myf.m函數(shù)文件，輸入以下內容：functionf=myf(x)f(1)=0.5*sin(x(1))+0.1*cos(x(2)*x(1))-x(1);f(2)=0.5*cos(x(1))-0.1*sin(x(2))-x(2);f=[f(1)f(2)];在MATLAB命令窗口中輸入：>>[r,n]=mulStablePoint([0,0])r=0.19790.4470n=11由計算結果可知，初始迭代值?。?，0）時，用11步迭代，得到原方程組的一組解（0.1979,0.4470）。10.3.2牛頓法設非線性方程組為，其中，牛頓法的迭代公式為：求解步驟為：（1）給出初始值；（2）對計算和；（3）求出，并進行精度控制。更一般的牛頓法迭代公式為：當時，就得到簡化牛頓法。在MATLAB中編程實現(xiàn)的非線性方程組的牛頓迭代法的函數(shù)為：mulNewton。功能：用牛頓迭代法求非線性方程組的一個解。調用格式：[r,n]=mulNewton(x0,eps)。其中，x0為初始迭代向量；eps為迭代精度；r為求出的解向量；n為迭代步數(shù)。牛頓法的MATLAB代碼如下：function[r,n]=mulNewton(x0,eps)%牛頓迭代法求非線性方程組的一個解%初始迭代向量：x0%迭代精度：eps%解向量：r%迭代步數(shù)：nifnargin==1eps=1.0e-4;endr=x0-myf(x0)/dmyf(x0);n=1;tol=1;whiletol>epsx0=r;r=x0-myf(x0)/dmyf(x0); %核心迭代公式tol=norm(r-x0);n=n+1;if(n>100000) %迭代步數(shù)控制disp('迭代步數(shù)太多，可能不收斂！');return;endend在MATLAB中編程實現(xiàn)的非線性方程組的簡化牛頓迭代法的函數(shù)為：mulSimNewton。功能：用簡化牛頓迭代法求非線性方程組的一個解。調用格式：[r,n]=mulSimNewton(x0,eps)。其中，x0為初始迭代向量；eps為迭代精度；r為求出的解向量；n為迭代步數(shù)。簡化牛頓法的MATLAB代碼如下：function[r,n]=mulSimNewton(x0,eps)%簡化牛頓迭代法求非線性方程組的一個解%初始迭代向量：x0%迭代精度：eps%解向量：r%迭代步數(shù)：nifnargin==1eps=1.0e-4;endr=x0-myf(x0)/dmyf(x0); %核心迭代公式c=dmyf(x0);n=1;tol=1;whiletol>epsx0=r;r=x0-myf(x0)/c;tol=norm(r-x0);n=n+1;if(n>100000) %迭代步數(shù)控制disp('迭代步數(shù)太多，可能不收斂！');return;endend【例10-19】牛頓法求解非線性方程組應用實例。采用牛頓法和簡化牛頓法求下面方程組的一個根。初始迭代值?。?,0）。解：首先建立myf.m函數(shù)文件，輸入以下內容：functionf=myf(x)f(1)=0.5*sin(x(1))+0.1*cos(x(2)*x(1))-x(1);f(2)=0.5*cos(x(1))-0.1*sin(x(2))-x(2);f=[f(1)f(2)];再建立dmyf.m導數(shù)的雅克比矩陣，輸入以下內容：functiondf=dmyf(x)df=[0.5*cos(x(1))-0.1*x(2)*sin(x(2)*x(1))-1,-0.1*x(1)*sin(x(2)*x(1));-0.5*sin(x(1)),-0.1*cos(x(2))-1];然后在MATLAB命令窗口中輸入：>>[r,n]=mulNewton(