精通MATLAB科学计算（第3版）(王正林)10-3r.pdf

《精通MATLAB科学计算（第3版）(王正林)10-3r.pdf》由会员分享，可在线阅读，更多相关《精通MATLAB科学计算（第3版）(王正林)10-3r.pdf（37页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第 章非线性方程求解非线性方程是常见的一类方程，非线性方程(组)的理论远不如线性方程(组)成熟和有效，特别是非线性方程组解的存在唯一性还没有完全解决，判断其解的存在性和解的个数几乎没有可行的方法。利用MATLAB提供的函数，可以求解一些简单的非线性方程;通过编程，可以解决一些较为复杂的非线性方程。通过本章学习，读者能够熟练掌握MATLAB中的非线性方程求解相关的函数，而且能通过编程实现多种求解非线性方程的数值算法。10.1MATLAB中非线性方程求根函数MATLAB中求解非线性方程的函数主要是以下两个：fzero和 fsolveo其中fzero可用来求解非线性方程，fsolve可用来求解非线性

2、方程和非线性方程组。10.1.1 fzero 函数fzero函数用于求解单变量非线性方程1&的基本用法及功能如下：1 .x=fzero(/；x 0)求函数f 在 x0附近的根；2 JV=fzero(/；x0,options)使用优化选项,options可选的值很多，具体可参考MATLAB帮 助；3.x,fval=fzero(/；x0,options)返回值中的 fv al=/(x)；4.x,fval,exitflag=fzero(f,xO,o ptions)其中 exitflag 标志着求解状态。函数fzero的具体用法见下面的例子。【例 10-1 非线性方程求解函数fzero的应用实例。采用

3、fcero函数求方程x2+x-l=0在 x=0.5附近的根。解：在 MATLAB命令窗口中输入下列命令：x=fzero(1xA2+x-l1,0.5)输出计算结果为：X=0.6180精通MATLAB科学计算(第2忡-从结果可知，方程在x=0.5附近的根为x=0.6180 o函数自er。还可以求函数在某个区间上的一个根，例 如，求取方程sin(x)-0.5=0 在0,2区间上的根，继续在MATLAB命令窗口中输入下列命令：x=fzero(*sin(x)-0.5*,0 2)输出计算结果为：X=0.5236从结果可知，方程在0,2区间上的根为x=0.5236。需要注意的是，当用也er。求函数在某个区间

4、上的一个根时，此函数在区间两端点处的值符号要相反才行。函数fzero还有一些特殊的用法，例 如，求取函数fx,y)=#-3 当y=2 时在x=2 附近 的 根，在 MATLAB命令窗口中输入下列命令：u=inline(xAy-3,x,y);号建立非线性方程/(X J)=*-3 c=fzero(uz 2,2)输出计算结果为：C=1.7321从结果可知，所求的根为x=1.7321。优化选项options的使用方法如下，使用0Ptimset函 数，其中参数TolX选项用来设置求解精度，在 MATLAB命令窗口中输入下列命令：opt=optimset(TolX*,1.Oe-8);x=fzero(sin

5、(x)-0.5,0.5,opt)x=0.5236 opt=optimset(*TolX1,1.Oe-2);x=fzero(*sin(x)-0.5,0.5,opt)x=0.5283从运算结果可以看出，不同精度对根有影响。10.1.2 fsolve 函数fsolve函数用于求解多变量非线性方程组。它的基本用法及功能如下：1 .x=fsolve(fun,x0)用 fun定义向量函数，其定义方式可以参考例子；2.x=fsolve(fiin,x0,options)options 为优化选项；3.x,val=fsolve(.)val=F(x),即函数值向量;220 第 1 0 章非线性方程求解4.x,va

6、l,exitflag=fsolve(.)exitflag 标志着求解状态；5.x,val,exitflag,output=fsolve(.)output 包含着优化后的结果信息；6.x,val,exitflag,output,jacobian=fsolve(.)jacobian 为解 x 处的 Jacobian 阵。其余参数与前面参数相似。另 外，符号数学工具箱中也有一个可用于求解方程组的函数solve，它可求各种类型方程组的解析解。当找不到解析解时，solve会自动寻求一个近似解，并且精度很高。需要注意的是，函数fsolve不是MATLAB符号工具箱中的函数，它位于优化工具箱内，而 solv

7、e是一个符号函数。【例 10-2】非线性方程组求解函数fsolve应用实例1。采 用 fsolve函数求方程组x：_ y _ 2 =0 在 a,y)=(2,2)附近的根。y2-2 x-4 =0解：首先新建一个.m 文 件，命名为myfun.m,键入下列内容：function F=myfun(x)F=x(l)A2-x(2)-2;x(2)A2-2*x(l)-4;在 MATLAB命令窗口中输入下列命令：x0=2 2;opt=optimset(*Display,iter1);r,val=fsolve(myfun,xO,opt)输出计算结果为：Optimization terminated:first-

8、order optimality is less than options.TolFun.Norm ofFirst-orderTrust-regionIterationFunc-count f(x)stepoptimalityradius0316161160.69900414.621290.0001181560.1239650.03562.53124.3305e-0110.003045211.96e-0052.54156.63221e-0241.8626e-0066.78e-0122.52.2143 2.9032val=1.0e-011*0.20750.1525由计算结果可知，方程组的一个根为

9、(x,y)=(2.2143,2.9032)。【例 10-31 非线性方程组求解函数fsolve应用实例2 o 采用fsolve函数求矩阵方程 221精通MATLAB科学计算（第 2 回1=-的 一 个 根，初始解向量为24 3解：首先新建一个.m 文 件，命名为myfun.m,键入下列内容：function F=matrixfun(x)F=x*x-ll,-8;24z-3;在 MATLAB命令窗口中输入下列命令：xO=l,l;lzl;opt=optimset(Display,1 iter);r=fsolve(Gmatrixfun,xO)输出计算结果为：Optimization terminate

10、d:first-order optimality is less than options.TolFun.r=1.0000-2.00006.0000 3.0000可以验证矩阵一 2 的平方刚好等于10.2其他数值求根法-11-824-3下面讲述的几种数值方法，在不加说明的情况下，都是在某个区间内求得方程的一个根。其理论依据为：如果/伍）/3）0,令,转 至（1 ）;如 果/（岁）=0，则 誓 为 一 个 根。（2）如 果a 审 0)disp(1两端点函数值乘积大于0)；return;elseroot=FindRoots(f,a,b,eps);电调用求解子程序endfunction r=Find

11、Roots(f,a,b,eps)f_l=subs(sym(f),findsym(sym(f),a);f_2=subs(sym(f),findsym(sym(f),b);223精通MATLAB科学计算(第 2 忡-mf=subs(sym(f),findsym(sym(f),(a+b)/2);%中点函数值if(f_l*mf0)t=(a+b)/2;r=FindRoots(f,t,b,eps);%右递归elseif(f_l*mf=O)r=(a+b)/2;elseif(abs(b-a)r=HalfInterval(TxA3-3*x+l,0,1)输出计算结果为：r=0.34730.3473可 见，方程d-

12、3 x +1 =()在区间 0,1上的一个根为x=0.3473。10.2.2黄金分割法二分法是把求解区间的长度逐次减半，而黄金分割法是把求解区间逐次缩短为前次的0.618倍。它的求解步骤如下。(1)设八=a+(l-0.618)*(b-a),=0 +0.618*(6-。)，且力=/(),力=/出)；(2 )如果,-司(给定的最小区间长度)，则 输 出 方 程 的 根 为 号；否则转到(3)；(3 )如果力*力0,令a=f2,反之令b=,转 到(1卜在M ATLAB中编程实现的黄金分割法的函数为：hjo功 能：用黄金分割法求函数在某个区间上的一个零点。调用格式:root=hj(f,a,b,eps)

13、。224 第10章非线性方程求解其 中，/为 函 数 名；a 为区间左端点；/7为区间右端点；eps为根的精度；root为求出的函数零点。黄金分割法的MATLAB程序代码如下：function root=hj(f,a,b,eps)黄金分割法求函数f在区间 a,b上的一个零点5t函数名：f为区间左端点；a8区间右端点：b义根的精度：eps亳求出的函数零点：rootif(nargin=3)eps=l.Oe-4;endfl=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(fl=0)root=a;endif(f2=0)

14、root=b;endif(fl*f20)disp(,两端点函数值乘积大于0!，)；return;elsetl=a+(b-a)*0.382;t2=a+(b-a)*0.618;f_l=subs(sym(f),findsym(sym(f),tl);f_2=subs(sym(f),findsym(s ym(f),t2);tol=abs(tl-t2);while(toleps)号精度控制if(f_l*f_20)a=t2;扬同号左端点调整 225精 通MATLAB科 学 计 算(第2忡-elseb=tl;%异号右端点调整endendtl=a+(b-a)*0.382;t2=a+(b-a)*0.618;f_l

15、=subs(sym(f),findsym(sym(f),tl);f_2=subs(sym(f),findsym(sym(f),t2);tol=abs(t2-tl);endroot=(tl+t2)/2;W 输出根end【例 10-5】黄金分割法求解非线性方程应用实例。采用黄金分割法求方程1-3 x +1 =0 在区间 0,1 上的一个根。解：在 MATLAB命令窗口中输入：r=hj(xA3-3*x+l,0,l)输出计算结果为：r=0.3473可 见，方程在区间 0,1 上的一个根为x=0.3473,与二分法求出来的结果相同。10.2.3不动点迭代法求方程/(x)=0 的 根，可以先把方程改写成如

16、下的形式：x=/(x)+x于是得到不动点迭代法的其中一种迭代公式：X=/(x-i)+x_1此种不动点迭代法很有可能不收敛，因为它的本质是求函数y=/(x)+x 与直线y=x的交点，而它们不一定存在交点。即使收敛，其速度也十分慢，因此有了艾特肯加速迭代与史蒂芬森加速迭代。在 MATLAB中编程实现的不动点迭代法的函数为：StablePointo功 能：用不动点迭代法求函数的一个零点。调用格式：root,n =StablePoint(f,x0,eps)。其 中，/为 函 数 名；M)为初始迭代向量；eps为根的精度；226 第10章非线性方程求解root为求出的函数零点；为迭代步数。不动点迭代法的

17、MATLAB程序代码如下：function root,n=StablePoint(f,xO,eps)常用不动点迭代法求函数的一个零点3初始迭代向量：x0卷根的精度：eps%求出的函数零点：root%迭代步数：nif(nargin=2)eps=l.Oe-4;endtol=l;root=x0;n=0;while(toleps)n=n+l;rl=root;root=subs(sym(f),f indsym(sym(f),rl)+rl;超 迭 代的核心公式tol=abs(root-rl);end【例10-6】不动点迭代法求解非线性方程应用实例。采用不动点迭代法求方程-y=+x-2 =0 的一个根。解：

18、在MATLAB命令窗口中输入：rz n=StablePoint(11/sqrt(x)+x-2*,0.5)r=0.3820nV .J7 JQ Jo 2nn=4二从 计 算 结 果 可 以 看 出，经 过4步 迭 代，得 出 方 程J=+x-2 =0的一个根为x=0.3820 o1.艾特肯加速迭代法艾特肯加速迭代法是在计算出X,X+l,X+2后，对 作 以 下 修 正：*(X+1 x”)x+i=-x“+2 2x+i+xnY 227精通MATLAB科学计算(第2忡-然后用焉+i来逼近方程的根。艾特肯加速迭代法是先用不动点迭代法算出一系列 4 ,再对此系列作修正得到系列xn,用后者来逼近方程的根。在

19、MATLAB中编程实现的艾特肯加速迭代法的函数为：AtkenStablePointo功 能：用艾特肯加速迭代法求函数的一个零点。调用格式:root,/?=AtkenStablePoint(f,xO,eps)o其 中，/为 函 数 名；M)为初始迭代向量；eps为根的精度；root为求出的函数零点；为迭代步数。艾特肯加速迭代法的MATLAB程序代码如下：function root,n=AtkenStablePoint(f,xO,eps)Z用艾特肯加速迭代法求函数f的一个零点常 初始迭代向量：xOZ根的精度：eps年求出的函数零点：r。”金迭代步数：nif(nargin=2)eps=l.Oe-4;

20、endtol=l;root=xO;x(1:2)=0;n=0;m=0;a2=x0;while(toleps)n=n+l;al=a2;rl=root;root=subs(sym(f),findsym(sym(f),rl)+rl;%求出根x(n)=root;w 把所有的根都保存下来if(n2)m=m+1;a2=x(m)-(x(m+1)-x(m)2/(x(m+2)-2*x(m+1)+x(m);%对根进行艾特肯修正tol=abs(a2-al);endendroot=a2;228 -第1 0 章非线性方程求解【例 10-7】艾特肯加速迭代法求解非线性方程应用实例。采用艾特肯加速迭代法求方程+x-2 =0

21、的一个根，迭代初始值为0.5。耳解：在 MATLAB命令窗口中输入：r,n=AtkenStablePoint(11/sqrt(x)+x-2,0.5)r=0.3820n=4-4从计算结果可以看出，经过4 步迭代，得出方程9+x-2 =0 的一个根为x=0.3820。这个方程比较简单，所以通过简单的几步迭代就可得出根。从下面的结果可以看出艾特肯加速迭代法的优点。rz n=AtkenStablePoint(*1/sqrt(x)+x-2 1 0.999)r=1.0000Ji-n7 o2 n7 2An=4 r,n=StablePoint(1/sqrt(x)+x-2 1 0.999)r=0.38200.3

22、820n=21w上面的例子说明采用初始值0.999时，经过4 步艾特肯加速迭代法，可以得到方程+x-2 =0 的另一个根x=l,而普通的不动点迭代法却得不到。Y 229精通MATLAB科学计算(第2忡-2.史蒂芬森加速迭代法史蒂芬森加速与艾特肯加速不同的地方在于前者是在迭代的同时就进行修正，它的迭代公式为：yn=/(x“)+Xz”=/仇)+%(%-X)?X+l=X-z“-2Ml+xn在 MATLAB中编程实现的史蒂芬森加速迭代法的函数为：StevenStablePointo功 能：用史蒂芬森加速迭代法求函数的一个零点。调用格式：root,=StevenStablePoint(/；.rO,eps

23、)o其 中，/为函数名；xO为初始迭代向量；eps为根的精度；root为求出的函数零点；为迭代步数。史蒂芬森加速的MATLAB程序代码如下：function root,n=StevenStablePoint(f,xO,eps)徒用史蒂芬森加速迭代法求函数f的一个零点%初始迭代向量：xO%根的精度：eps学求出的函数零点：root为迭代步数：nif(nargin=2)eps=l.Oe-4;endtol=l;root=xO;n=0;while(toleps)n=n+l;rl=root;y=subs(sym(f),findsym(sym(f),rl)+rl;z=subs(sym(f),findsym

24、(sym(f),y)+y;root=rl-(y-rl)A2/(z-2*y+rl);%对每次算出的根立即进行修正tol=abs(root-rl);end【例 10-8】史蒂芬森加速迭代法求解非线性方程应用实例。采用史蒂芬森加速迭代230 第10章非线性方程求解法求方程9+x-2 =0 的一个根。解：在 MATLAB命令窗口中输入：r,n=StevenStablePoint(1/sqrt(x)+x-2,0.999)r=1.0000一i Vn nV oV oVn=22 r,n=StevenStablePoint(11/sqrt(x)+x-2z 0.5)r=0.3820n=44对比上例可以看出，史蒂芬

25、森加速迭代法不仅能求出方程1忑+x-2 =0 的一个根x=0.3820,而且它比艾特肯加速迭代法更快地求出另一个根x=l。10.2.4弦截法弦截法的算法过程如下：(1)过两点(a,/(a),(6 J 3)作一直线，它与X轴有一个交点，记 为 汨；(2 )如果fafx)0,过两点(aj(a),(x1,/(的)作一直线，它与X轴的交点记为,否则过两点(6,7(b),(X ,/但)作一直线，它与x 轴的交点记为打；(3)如此下去，直到,就可以认为X“为/(尤)=0 在区间a,b上的一/根。X k=a-(4)为,的递推公式为：.Xk=b-x i -aXk-i-b/()/(%*._I)-/(6)/3)/

26、W(x i)0且 X i=a-/(a)o在 MATLAB中编程实现的弦截法的函数为：Secanto功 能：用弦截法求函数在某个区间上的一个零点。调用格式:root=Secant(f,a,h,eps)o其 中，/为 函 数 名；a 为区间左端点；6 为区间右端点；eps为根的精度；Y 231精通MATLAB科学计算(第2贿 root为求出的函数零点。弦截法的MATLAB程序代码如下：function root=Secant(fz a,b,eps)%弦截法求函数f在 区 间 上 的 一 个 零 点函数名：f他区间左端点；a%区间右端点：b为根的精度：eps星求出的函数零点：rootif(nargi

27、n=3)eps=l.Oe-4;endfl=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(fl=0)root=a;endif(f2=0)root=b;endif(fl*f20)disp(两端点函数值乘积大于0 1;return;elsetol=l;fa=subs(sym(f),findsym(sym(f),a);fb=subs(sym(f),findsym(sym(f),b);root=a-(b-a)*fa/(fb-fa);while(toleps)rl=root;fx=subs(sym(f),findsym(

28、sym(f),rl);s=fx*fa;if(s=0)root=rl;elseif(s0)root=b-(rl-b)*fb/(fx-fb);elseroot=a-(rl-a)*fa/(fx-fa);endendtol=abs(root-rl);endend多迭代初始值232 -第1 0 章非线性方程求解【例 10-9】弦截法求解非线性方程应用实例。采用弦截法求方程lgx+&=2 在区间 1,4 上的一个根。解：在 MATLAB命令窗口中输入：r=Secant(*log(x)+sqrt(x)-21,lz 4)输出计算结果为：r=1.87731.8773由计算结果可知，lgx+&=2 在区间 1,4

29、 上的一个根为x=1.8773。10.2.5史蒂芬森弦截法史蒂芬森弦截法是弦截法的一种变形，它的递推公式为：Xk=Xk-f(XA-1)/(Xz+f g-)八且有x=a-f(a)o史蒂芬森弦截法比弦截法的精度要高，收敛速度要快。在 MATLAB中编程实现的史蒂芬森弦截法的函数为：StevenSecanto功 能：用史蒂芬森弦截法求函数在某个区间上的一个零点。调用格式:root=StevenSecant(f,a,b,eps)o其 中，/为函数名；a 为区间左端点；6 为区间右端点；eps为根的精度；root为求出的函数零点。史蒂芬森弦截法的MATLAB程序代码如下：function root=St

30、evenSecant(f,a,b,eps)亳史蒂芬森弦截法求函数f在区间 a,b上的一个零点会函数名：f 区间左端点；a彩区间右端点：b国根的精度：eps省求出的函数零点：rootY 233精通MATLAB科学计算(第2忡-if(nargin=3)eps=l.Oe-4;endfl=subs(sym(f),findsym(sym(f)z a);f2=subs(sym(f),findsym(sym(f),b);if(fl=O)root=a;endif(f2=0)root=b;endif(fl*f20)disp C两端点函数值乘积大于0!，)；return;elsetol=l;fa=subs(sym

31、(f),findsym(sym(f),a);fb=subs(sym(f),findsym(sym(f),b);faa=subs(s ym(f),findsym(sym(f),fa+a);root=a-fa*fa/(faa-fa);多迭代初始值while(toleps)rl=root;fx=subs(sym(f),findsym(sym(f),rl);v=fx+rl;fxx=subs(sym(f),findsym(sym(f),v);root=rl-fx*fx/(fxx-fx);S递推的核心公式tol=abs(root-rl);endend【例10-10史蒂芬森弦截法求解非线性方程应用实例。采用

32、史蒂芬森弦截法求方程lgx+=2在区间 1,4上的一个根。解：在MATLAB命令窗口中输 入：r=StevenSecant(*log(x)+sqrt(x)-2,1,4)输出计算结果为：r=1一 工结果出现错误，改正方法可以是把求解区间缩短，缩短为 1.5,2：r=StevenSecant(*log(x)+sqrt(x)-2,1.5,2)r=1.8773-3这样就得出了正确的根1.8773。从上面的例子可以看出求解区间的选择也是很重要的，一般来说，求解区间越短越好。234 -第；0章 非线性方程求解1 0.2.6抛物线法弦截法其实是用不断缩短的直线来近似函数/（X）,而抛物线法则采用三个点来近似

33、函数/（X），并且用抛物线与横轴的交点来逼近函数./a）的根。它的算法过程如下：（1 ）选定初始值工0,两,切,并计算/（%0）,/（两）,/（x2）和以下差分:即-Xo/1电 的,知=/8 闻一/氏项X2-XQ一般取X O=Q,x=b,a X2 b o注意不要使三点共线。（2）用牛顿插值法对三点（工0,/（工0）,（两,/（两）,（工2 J（X2）进行插值得到一条抛物线,它有两个根：%3=+-gVg2-4AC2C其中A=f(x2),C=fX2,X,Xo,B=fx2,X1+fx2,X|,Xo(x2-XI)。两个根中只取靠近X2的那个根，即土号取与8同 号，即2AX3=X2-f fB+sgn(B

34、)ylB2-4AC（3）用百户2,看代替X。,即,X2,重复以上步骤，并有以下递推公式：2AnB,i+sgn(z/)J B,r-4AnCn其中A”=,/(x),C”=fxn,Xn-,Xn-2,Bn=fxn,Xn-+fxn Xn-2(x-Xn-)（4）进行精度控制。在M A T L A B中编程实现的抛物线法的函数为：P arabolao功 能：用抛物线法求函数在某个区间上的一个零点。调用格式：root=P arabola（f,a,h,x,eps）o词 0)disp(两端点函数值乘积大于0!，)；return;elsetol=l;fa=subs(sym(f),findsym(sym(f),a);

35、fb=subs(sym(f),findsym(sym(f),b);fx=subs(sym(f),findsym(sym(f),x);dl=(fb-fa)/(b-a);d2=(fx-fb)/(x-b);d3=(d2-dl)/(x-a);B=d2+d3*(x-b);root=x-2*fx/(B+sign(B)*sqrt(BA2-4*fx*d3)t=zeros(3);t(1)=a;236 第10章非线性方程求解t(2)=b;t(3)=x;while(toleps)t(l)=t(2);t(2)=t(3);t(3)=root;fl=subs(sym(f),findsym(sym(f),t(1);f2=s

36、ubs(sym(f),findsym(sym(f)z t(2);f 3=subs(sym(f),f indsym(sym(f),t)；dl=(f2-fl)/(t(2)-t(l);d2=(f3-f2)/(t(3)-t(2);d3=(d2-dl)/(t(3)-t(l);B=d2+d3*(t(3)-t(2);称保存三个点%计算三点的函数值告计算三个差分%计算算法中的Broot=t(3)-2*f3/(B+sign(B)*sqrt(BA2-4*f3*d3);tol=abs(root-t(3);endend【例 10-11】抛 物 线 法 求 解 非 线 性 方 程 应 用 实 例 I。采用抛物线法求方程

37、lgx+4 =2 在区间 1,4 上的一个根，迭代初始值为2。解：在 MATLAB命令窗口中输入：“Parabola(sqrt(x)+log(x)-2,1,4,2)r=1.87731.8773【例 10-12】抛物线法求解非线性方程应用实例2o分别从两个区间 0.5,1.5 和 0.1,1上采用抛物线法求方程1忑+x-2 =0 的根。解：在 MATLAB命令窗口中输入：r=Parabola(1/sqrt(x)+x-2 0.5,1.5,0.8)告迭代区间为0 5,1.5 ,初始值为 0.8r=1工 r=Parabola(*1/sqrt(x)+x-2*z0.1zl,0.5)带迭代区间为0.1,1

38、,初始值为 0.5r=0.38201忑由此可以看出抛物线法也能求出方程+x-2 =0 的两个根。1 0.2.7牛顿法弦截法本质上是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为：Y 237精通MATLAB科 学 计 算(第2回f M初始值可以取/(幻和/(b)的较大者，这样可以加快收敛速度。和牛顿法有关的还有简化牛顿法和牛顿下山法。在 MATLAB中编程实现的牛顿法的函数为：NewtonRooto功 能：用牛顿法求函数在某个区间上的一个零点。调用格式：root=NewtonRoot(f,a,b,eps)。其 中，/为函数名；。

39、为区间左端点；6 为区间右端点；eps为根的精度；root为求出的函数零点。牛顿法的MATLAB程序代码如下：function root=NewtonRoot(f,a,b,eps)W 牛顿法求函数f 在区间 a,b 上的一 零点专函数名：f得区间左端点；a专区间右端点：b传根的精度：eps务求出的函数零点：rootif(nargin=3)eps=l.Oe-4;endfl=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(fl=0)root=a;endif(f2=0)root=b;endif(fl*f20)disp

40、(1两端点函数值乘积大于0!，)；return;elsetol=l;fun=diff(sym(f);fa=subs(sym(f),findsym(sym(f),a);%求导数238 第10章非线性方程求解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(dfadfb)root=a-fa/dfa;elseroot=b-fb/dfb;endwhile(toleps)rl=root;fx=subs(sym(f),findsym(

41、sym(f),rl);dfx=subs(sym(fun),findsym(sym(fun),rl);考初始值取两端点导数较大者root=rl-fx/dfx;tol=abs(root-rl);endend为求该点的导数值告迭代的核心公式【例10-13】牛顿法求解非线性方程应用实例。采用牛顿法求方程4-丫3+2=0在区间0.5,2上的一个根。解：在M A T L A B命令窗口中输入：r=NewtonRoot(*sqrt(x)-xA3+2*,0.5,2)输出计算结果为：r=1.4759.(759由计算结果可知，/+2=0的一个根为x=.4759。需要注意的是，初始值的选择不要使得其导数为0。1.简

42、化牛顿法将牛顿法的递推公式改为：X+l=x-A f(xn)为保证收敛，系数兀只须满足0彳一 一。如果7取常数，则称为简化牛顿/(X)./(Xo)法。在M A T L A B中编程实现的简化牛顿法的函数为：SimpleNewtono功 能：用简化牛顿法求函数在某个区间上的一个零点。调用格式：root=SimpleNewton(/；a,b,eps)o其 中，/为函数名；“为区间左端点；6为区间右端点；Y 239精 通MATLAB科学计算(第2忡-eps为根的精度；root为求出的函数零点。简化牛顿法的MATLAB代码如下：function root=SimpleNewton(f,a,b,eps)专

43、简化牛顿法求函数f在区间 a,b上的一个零点得函数名：f阳区间左端点；a专区间右端点：b传根的精度：eps国求出的函数零点：rootif(nargin=3)eps=l.Oe-4;endfl=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(fl=0)root=a;endif(f2=0)root=b;endif(fl*f20)disp(两端点函数值乘积大于0 1;return;elsetol=l;fun=diff(sym(f);%求导数fa=subs(sym(f),findsym(sym(f),a);fb=subs

44、(sym(f),findsym(sym(f),b);dfa=subs(sym(fun),findsym(sym(fun),a);dfb=subs(sym(fun),findsym(sym(fun),b);if(dfadfb)%初始值取两端点导数较大者df0=l/dfa;root=a-df0*fa;elsedfO=l/dfb;root=b-df0*fb;endwhile(toleps)rl=root;fx=subs(sym(f),findsym(sym(f),rl);root=rl-dfO*fx;当迭代的核心公式，其中dfO为常数240 第 1 0 章非线性方程求解tol=abs(root-rl

45、);endend【例 10-14】简化牛顿法求解非线性方程应用实例。采用简化牛顿法求方程yx-x3+2=0 在区间 1.2,2 上的一个根。解：在 MATLAB命令窗口中输入：r=SimpleNewton(*sqrt(x)-x人 3+2,1.2,2)输出计算结果为：r=1.47591.4759由计算结果可知，五-/+2=。的一个根为X=I.4759。2.牛顿下山法牛顿下山法收敛速度很快，很容易就可得出方程的根。牛顿下山法的递推公式为：x“+i=x-a (Ova,1)f (x)a 称为下山因子，它的确定方法如下：(1)先取a =l,求出;(2 )判断条件是否满足；不满足的话令c =0.5a;(4

46、 )不断验证直到|/(瑞+|)卜|/。)|。在 MATLAB中编程实现的牛顿下山法的函数为：NewtonDowno功 能：用牛顿下山法求函数在某个区间上的一个零点。调用格式:root=NewtonDown(f,a,b,eps)o其 中，/为函数名；。为区间左端点；6 为区间右端点；eps为根的精度；root为求出的函数零点。牛顿下山法的MATLAB程序代码如下：function root=NewtonDown(f,a,b,eps)W牛顿下山法求函数f 在区间 a,b 上的一个零点区函数名：f与区间左端点；aB区间右端点：b 241精通MATLAB科学计算(第2忡-与根的精度：eps%求出的函数

47、零点：rootif(nargin=3)eps=l.Oe-4;endfl=subs(sym(f),findsym(sym(f),a);f2=subs(sym(f),findsym(sym(f),b);if(fl=0)root=a;endif(f2=0)root=b;endif(fl*f20)disp(1两端点函数值乘积大于0!);return;elsetol=l;fun=diff(s ym(f);fa=subs(sym(f),findsym(sym(f),a);fb=subs(s ym(f),findsym(sym(f),b);dfa=subs(sym(fun),findsym(sym(fun)

48、,a);dfb=subs(sym(fun),findsym(sym(fun),b);if(dfadfb)金迭代初始值的选取root=a;elseroot=b;endwhile(toleps)rl=root;fx=subs(sym(f),findsym(sym(f),rl);dfx=subs(sym(fun),findsym(sym(fun),rl);toldf=l;alpha=2;while toldf0%此循环为下山因子的选取过程alpha=alpha/2;root=rl-alpha*fx/dfx;会迭代的核心公式fv=subs(sym(f),findsym(s ym(f),root);to

49、ldf=abs(fv)-abs(fx);endtol=abs(root-rl);endend【例 10-15】牛顿法下山求解非线性方程应用实例。采用牛顿下山法求方程V x-x3+2=0 在区间 1.2,2 上的一个根。242 第10章非线性方程求解解：在 MATLAB命令窗口中输入：r=NewtonDown(sqrt(x)-xA3+2 1 12,2)输出计算结果为：r=1.47591.4759由计算结果可知，6-/+2 =0 的一个根为x=1.4759。10.2.8两步迭代法两步迭代法的一般公式为：y=g(x)x“+i=h(yn)它比较常用的有以下两种形式：(1)(2)/(%,)(此一/)在

50、MATLAB中编程实现的两步迭代法的函数为：TwoStep。功 能：用两步迭代法求函数在某个区间上的一个零点。调用格式：root=TwoStep(f,a,b,type,eps)0其 中，/为 函 数 名；“为区间左端点；6 为区间右端点；type为两步迭代法的形式；eps为根的精度；root为求出的函数零点。两步迭代法的MATLAB代码如下：function root=TwoStep(f,a,b,type,eps)S两步迭代法求函数f在区间 a,b 上的一个零点覆函数名：f务区间左端点；a 243精通MATLAB科学计算(第2版i-南区间右端点：b专两步迭代法的形式：type(可 取 1 ,2