2022年第二章非线性方程的数值解法的matlab程序 .pdf

《2022年第二章非线性方程的数值解法的matlab程序 .pdf》由会员分享，可在线阅读，更多相关《2022年第二章非线性方程的数值解法的matlab程序 .pdf（23页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等） ，求根的方法（二分法，迭代法，牛顿法，割线法，米勒（M ller）法和迭代法的加速等）及其MATLAB 程序，求解非线性方程组的方法及其MATLAB 程序 . 2.1 方程（组）的根及其MATLAB 命令2.1.2 求解方程（组）的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

2、(x1, ,xn)=q1(x1, ,xn) ); . En=sym( 方程 fn(x1, ,xn)=qn(x1, ,xn) ); x1,x2, ,xn=solve(E1,E2,En, x1, ,xn) 2.1.3 求解多项式方程（组）的roots 命令如果)(xf为多项式，则可分别用如下命令求方程0)(xf的根，或求导数)(xf（见表 2-1）. 表 2-1 求解多项式方程（组）的roots命令命 令功能xk =roots(fa) 输入多项式)(xf的系数fa（按降幂排列） ，运行后输出xk为0)(xf的全部根 . dfa=polyder(fa) 输入多项式)(xf的系数fa（按降幂排列） ，

3、运行后输出 dfa为多项式)(xf的导数)(xf的系数 .dfx=poly2sym(dfa) 输入多项式)(xf的导数)(xf的系数dfa（按降幂排列） ，运行后输出dfx为多项式)(xf的导数)(xf. 2.1.4 求解方程（组）的 fsolve 命令如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的 roots 命令 . 如果非线性方程（组）是含有超越函数，则无法使用roots 命令，需要调用 MATLAB 系统中提供的另一个程序fsolve 来求解 .当然，程序 fsolve 也可以用于多项式方程（组），但是它的计算量明显比roots 命令的大 . f

4、solve 命令使用最小二乘法（least squares method）解非线性方程（组）(FX )0的数值解，其中X 和 F（X）可以是向量或矩阵.此种方法需要尝试着输入解X 的初始值（向量或矩阵） X0， 即使程序中的迭代序列收敛，也不一定收敛到(FX )0的根 （见例 2.1.8 ） . fsolve的调用格式 ： X=fsolve(F,X0)输入函数)(xF的 M 文件名和解X的初始值（向量或矩阵）X0，尝试着解方程（组）第二章非线性方程（组）的数值解法精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 23 页(FX )0，运行

5、后输出(FX )0解的估计值（向量或矩阵）X.要了解更多的调用格式和功能请输入: help fsolve , 查看说明 . 2.2 搜索根的方法及其MATLAB 程序求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）. 搜索根的近似位置的常用方法有三种: 作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理 . 2.2.1 作图法及其 MATLAB 程序作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等. 下面介绍两种作图程序.

6、 作函数)(xfy在区间 a,b 的图形的 MATLAB 程序一x=a:h:b; % h 是步长y=f(x); plot(x,y) grid, gtext(y=f(x)说明：此程序在MATLAB 的工作区输入，运行后即可出现函数)(xfy的图形 . 此图形与x轴交点的横坐标即为所要求的根的近似值. 区间 a,b 的两个端点的距离b-a 和步长h的绝对值越小，图形越精确. 作函数)(xfy在区间 a,b 上的图形的MATLAB 程序二将)(xfy化为)()(xgxh，其中)()(xgxh和是两个相等的简单函数x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2

7、) grid,gtext( y1=h(x),y2=g(x)说明：此程序在MATLAB 的工作区输入， 运行后即可出现函数)()(21xgyxhy和的图形 . 两图形交点的横坐标即为所要求的根的近似值. 下面举例说明如何用计算机软件MATLAB 的图形功能作图. 2.2.2 逐步搜索法及其MATLAB 程序逐步搜索法也称试算法. 它是求方程0)(xf根的近似值位置的一种常用的方法. 逐步搜索法依赖于寻找连续函数)(xf满足)(af与)(bf异号的区间,ba. 一旦找到区间，无论区间多大，通过某种方法总会找到一个根. MATLAB的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的

8、收敛判定准则编写其主程序，命名为zhubuss.m.逐步搜索法的MATLAB 主程序输入区间端点a 和 b 的值，步长h和精度tol，运行后输出迭代次数k=( b-a)/ h+1,方程0)(xf根的近似值r.function k,r=zhubuss(a,b,h,tol) % 输入的量 - a 和 b是闭区间 a,b 的左、右端点；%-h 是步长；%-tol 是预先给定的精度. % 运行后输出的量-k 是搜索点的个数；% - 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);

9、for k=2:n X(k)=a+k*h;Y(k)=funs(X(k); %程序中调用的 funs.m 为函数精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 23 页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 k,r=zhubuss(-2,2,0.001,0.0001)运行后输出的结果k =4001 r = -1.2240 -1.0000 -1.0000 -0.9990 1.

10、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.* x.3+2.* x.2-3.* x-3用y=sin(cos(2.* x.3) 代替，可得到方程在区间2 ,2上的根的近似值如下r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.9230 1.3310 1.5780 1.7650 1.9200 如果读者分别将方程的结果与例2.2.3 比较， 方程的结果与例2.1.2 比较，将会发现逐步

11、搜索法的MATLAB 程序的优点 . 如果精度要求比较高，用这种逐步搜索法是不合算的.2.3 二分法及其 MATLAB 程序2.3.1 二分法二分法也称 逐次分半法 . 它的基本思想是：先确定方程0)(xf含根的区间（a,b） ，再把区间逐次二等分. 我们可以根据式 （2.3b ） 、 区间 a,b 和误差， 编写二分法求方程根的迭代次数的MATLAB命令 . 已知闭区间 a,b 和误差. 用二分法求方程误差不大于的根的迭代次数k的MATLAB 命令k=-1+ceil(log(b-a)- log(abtol)/ log(2) % ceil 是向 +方向取整， abtol 是误差. 2.3.2

12、二分法的 MATLAB 程序二分法需自行编制程序，现提供用二分法求方程f( x)=0 的根*x的近似值kx的步骤和式（2.3a ）编写一个名为erfen.m 的二分法的MATLAB主程序如下 . 二分法的MATLAB 主程序求解方程0)(xf在开区间（ a,b）内的一个根的前提条件是)(xf在闭区间 a,b上连续，且0)()(bfaf. 输入的量： a和 b是闭区间 a,b的左、右端点，abtol是预先给定的精度. 运行后输出的量：k 是使用二分法的次数.x 是方程在（ a,b）内的实根x*的近似值，其精度是abtol.wuca=|bk-ak|/2 是使用 k 次二分法所得到的小区间的长度的一

13、半，即实根 x*的近似值x 的绝对精度限 ,满足 wuca abtol.yx=f(xk) ，即方程 f(x)=0 在实根 x*的近似值x精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 23 页处的函数值 .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* yb0, disp( 注意： ya*yb0, 请重新调整区间端点a和b.), return end max1=-1+ceil(l

14、og(b-a)- log(abtol)/ log(2); % ceil 是向方向取整for k=1: max1+1 a;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,yx if yx=0 a=x; b=x;elseif yb*yx0 b=x;yb=yx; else a=x; ya=yx; end if b-ax=-4:0.1:4; y=x.3-x +4; plot(x,y) grid,gtext(y=x3-x+4) 画出函数f(x)=x3-x+ 4 的图像，如图2

15、-7. 从图像可以看出，此曲线有两个驻点33都在x轴的上方，在(-2,-1)内曲线与x 轴只有一个交点，则该方程有唯一一个实根，且在 (-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) 满足0)1()2(ff，所以此方程在(-2,-1)内有一个实根.（2）用二分法的主程序计算. 在 MATLAB 工作窗口输入程序 k,x,wuca,yx=erfen (-2,-1,0.001) 运行后

16、屏幕显示用二分法计算过程被列入表 2-3 ，其余结果为k = 9 ，x=-1.7959 ，wuca = 9.7656e-004，yx = 0.0037图 2-7精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 23 页2.4 迭代法及其 MATLAB 程序确定根的近似位置以后，接下来的工作就是将根精确化，即按某种方法将初始近似值逐步精确化， 直到满足所要求的精确度为止. 上节介绍的二分法是将根精确化的方法之一，但是它的收敛速度较慢，且不能求出偶重根. 迭代法可以克服这种缺陷. 迭代法是求解方程的根、线性和非线性方程组的解的基本而重要的方法

17、. 2.4.2 迭代法的 MATLAB 程序 1 迭代法需要自行编制程序. 下面提供的迭代法的MATLAB 程序 1 使用时只需输入迭代初始值0 x、迭代次数k、迭代公式xk+ 1=( xk) 和一条命令，运行后就可以输出求迭代序列kx、迭代k次得到的迭代值kx和相邻两次迭代的偏差piancha =|1kkxx| （简称 偏差 ）和偏差的相对误差xdpiancha=|1kkxx|kx的值，并且具有警报功能（若迭代序列发散，则提示用户“请重新输入新的迭代公式”；若迭代序列收敛，则屏幕会出现“祝贺您！此迭代序列收敛，且收敛速度较快”）. 我们可以用这个程序来判断迭代序列的敛散性，也可以用于比较由一

18、个方程得到的几个迭代公式的敛散性的优劣. 迭代法的MATLAB 程序 1 输入的量：初始值0 x、迭代次数k和迭代公式), 2 , 1 ,0()(1kxxkk；运行后输出的量：迭代序列kx、迭代k次得到的迭代值kx、相邻两次迭代的偏差piancha =|1kkxx| 和它的偏差的相对误差xdpiancha=kkkxxx1的值 .根据迭代公式（2.4 ）和已知条件，现提供名为diedai1.m 的M文件如下function k,piancha,xdpiancha,xk=diedai1(x0,k) % 输入的量 -x0 是初始值 ,k是迭代次数x(1)=x0; for i=1:k x(i+1)=f

19、un1(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 xk end if (piancha 1)&(xdpiancha0.5)&(k3) disp( 请用户注意：此迭代序列发散，请重新输入新的迭代公式) return ; end if (piancha 0.001)&(xdpiancha3) disp( 祝贺您！此迭代序列收敛，且收敛速度较快) return ; end p=(i

20、-1) piancha xdpiancha xk; 例 2.4.1求方程102)(2xxxf的一个正根 . 解在 MATLAB 工作窗口输入程序 k,piancha,xdpiancha,xk= diedai1(2,5)运行后输出用迭代公式2/ )10(21kkxx的结果k,piancha,xdpiancha,xk= 1.00000000000000 1.00000000000000 0.33333333333333 3.00000000000000 2.00000000000000 2.50000000000000 5.00000000000000 0.50000000000000 3.000

21、00000000000 4.37500000000000 0.89743589743590 4.87500000000000 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 5 页，共 23 页4.00000000000000 11.75781250000000 1.70828603859251 -6.88281250000000 5.00000000000000 11.80374145507813 0.63167031671297 -18.68655395507813 请用户注意：此迭代序列发散，请重新输入新的迭代公式k=5 ，piancha =

22、 11.80374145507813，xdpiancha = 0.63167031671297，xk = -18.68655395507813由以上运行后输出的迭代序列与根*x=2.316 624 790 355 40相差越来越大，即迭代序列kx发散，此迭代法就失败. 这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于 0.5. 用迭代公式)2/(101kkxx运行后输出的结果k,piancha,xdpiancha,xk= 1.00000000000000 0.50000000000000 0.20000000000000 2.50000000000000 2.00000

23、000000000 0.27777777777778 0.12500000000000 2.22222222222222 3.00000000000000 0.14619883040936 0.06172839506173 2.36842105263158 4.00000000000000 0.07926442612555 0.03462603878116 2.28915662650602 5.00000000000000 0.04230404765128 0.01814486863115 2.33146067415730 k =5,piancha =0.04230404765128,xdpi

24、ancha = 0.01814486863115, xk = 2.33146067415730可见， 偏差piancha和偏差的偏差的相对误差xdpiancha的值逐渐变小， 且第 5 次的迭代值 xk=2.331 460 674 157 30与根*x=2.316 624 790 355 40接近，则迭代序列kx收敛，但收敛速度较慢，此迭代法较为成功. 用迭代公式)22/()102(21kkkkkxxxxx运行后输出的结果k,piancha,xdpiancha,xk= 1.00000000000000 0.33333333333333 0.14285714285714 2.3333333333

25、3333 2.00000000000000 0.01666666666667 0.00719424460432 2.31666666666667 3.00000000000000 0.00004187604690 0.00001807631822 2.31662479061977 4.00000000000000 0.00000000026437 0.00000000011412 2.31662479035540 5.00000000000000 0 0 2.31662479035540 祝贺您！此迭代序列收敛，且收敛速度较快k = 5; piancha = 0; xdpiancha = 0;

26、 y = 2.31662479035540.可见，偏差piancha 和偏差的相对误差xdpiancha 的值越来越小，且第5 次的迭代值xk=2.316 624 790 355 40 与根*x=2.316 624 790 355 40 相差无几，则迭代序列kx收敛，且收敛速度很快，此迭代法成功. 2.4.5 迭代法的 MATLAB 程序 2 迭代法的MATLAB 程序 2 输入的量：初初始值0 x、精度tol和和迭代公式),2, 1 , 0()(1kxxkk；运行后输出的量：迭代次数k，满足精度tol的根的近似根kx的迭代序列kx，相邻 两 次 迭 代 的 偏 差pianch=|1kkxx|

27、 和 它 的 偏 差 的 相 对 误 差xdpiancha=kkkxxx1的值， yk= (xk).根据迭代公式（2.4 ）和已知条件，现提供名为diedai2.m 的 M 文件如下function 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

28、 xk yk if (pianchatol)|(xdpianchaddmax 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 23 页disp( 迭代次数超过给定的最大值ddmax ) return ; end P=(i-1),piancha,xdpiancha,xk,yk; 例2.4.5 求 x5-3x+1=0在 0.3附近的根，精确到4位小数 .解 构造迭代公式)(kkxx. 改写原方程为等价方程3/) 1(5xx.这时3/ ) 1()(5xx，初始值为 x0=0.5, 迭代公式3/ ) 1(51kkxx),2, 1,0(k. 利用

29、迭代法的MATLAB 程序 2计算精确到 4位小数的根的近似值y 和迭代次数 k. 在MATLAB 工作窗口输入程序 k,piancha,xdpiancha,xk,yk=diedai2(0.3,1e-4,100)运行后输出的结果k,piancha,xdpiancha,xk,yk = 0 0.03414 0.10218 0.30000 0.33414 1 0.03414 0.10218 0.33414 0.33472 2 0.00058 0.00173 0.33472 0.33473 3 0.00001 0.00004 0.33473 0.33473 k =3;piancha =1.206089

30、525390697e-005; xdpiancha =3.603129477781680e-005; xk =0.3347; yk =0.3347. 2.5 迭代过程的加速方法及其MATLAB 程序对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的问题.2.5.2 加权迭代法的 MATLAB 程序加权迭代法的MATLAB 主程序已知初始值0 x、 精度tol和迭代公式xk+1=(xk) ， 求满足精度tol的根的近似根kx和它的函数值yk及迭代次数k.根据加权迭代的公式（2. 10）和已知条件，现提供名

31、为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 (pianchatol)|(xdpianchaddmax disp( 迭代次数超过给定的最大值dd

32、max ) return ; end P=(i-1),xk,yk; 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 23 页例 2.5.1 用（2.10 ）式求 2e0 xx在 0.85 附近的一个近似根，要求精度610. 解在 MATLAB 工作窗口输入程序 k,xk,yk=jasudd (0.85,1e-6,-0.855,100)运行后输出结果k,xk,yk = 1.00000000000000 0.85260370028041 0.85260703830561 2.00000000000000 0.85260549975491 0

33、.85260550406236 3.00000000000000 0.85260550207699 0.85260550208255 k =3 ； xk =0.852606； yk = 0.852606.2.5.4 艾特肯 (Aitken)加速方法的 MATLAB 程序艾特肯加速方法的MATLAB 程序已知初始值0 x、 精度tol和迭代公式xk+1=(xk) ， 求满足精度tol的根的近似根kx和它的函数值yk及迭代次数k，p= k ,x 1,x2,x .根据艾特肯加速方法的公式（2.11）和已知条件，现提供名为Aitken.m 的 M文件如下：function k,xk,yk,p= Ait

34、ken (x0,tol, ddmax) x(1)=x0; 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 (pianchatol)|(xdpianchaddmax disp( 迭代次数超过给定的最大值ddmax ) return ;

35、end m=0,1:i-1; p=m,x1,x2,x; 例2.5.3 用艾特肯加速方法求2e0 xx在 0.85 附近的一个近似根，要求精度610. 解在 MATLAB 工作窗口输入程序 k,xk,yk,p=Aitken(0.85,1e-6, 100)运行后输出结果k=3 ，xk=0.85260550201343，yk=0.85260550201373 p = 0 0 0 0.85000000000000 1.00000000000000 0.85482986389745 0.85071110652484 0.85260683568607 2.00000000000000 0.85260436

36、491811 0.85260647150826 0.85260550201407 3.00000000000000 0.85260550201343 0.85260550201398 0.852605502013732.6 牛顿(Newton) 切线法及其 MATLAB 程序2.6.2 牛顿切线法的收敛性及其MATLAB 程序判别牛顿切线法的局部收敛性的MATLAB 主程序精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 23 页根据牛顿切线法局部收敛的条件. 1)()()()(2*xfxfxfx输入的量： x 是方程0)(xf根*x或满

37、足1*0 xx的初始值0 x， 函数)(xf及其一阶导数)( xf和二阶导数)( xf；输出的量：)(*x或)(0 x的值 . 如果)(*x1或1)(0 x时，运行后输出信息恭喜您！此迭代序列收敛, (x)的导数值的绝对值y=|d(x)/dx|和方程f(x)=0的函数 f(x)值 f 如下 ；否则，运行后输出信息 请注意观察下面显示的(x)的导数值的绝对值 y=|d(x)/dx|和方程 f(x)=0的函数 f(x)值.根据牛顿切线法局部收敛的条件和方程0)(xf的函数)(xf及其一、二阶导数，现提供程序：function y,f=newjushou(x) f=fnq(x); fz=fnq(x)

38、*ddfnq(x)/(dfnq(x)2+eps); y=abs(fz); if (y y,f=newjushou(-1)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d (x)/dx| 和方程 f(x)=0 的函数 f(x)值y = 139.5644 f = 4.3096 说明 ：实际上，这个初始值所产生的迭代序列kx发散 , 见例 2.6.6. （2）输入程序 y,f=newjushou(0)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d (x)/dx| 和方程 f(x)=0 的函数 f(x)值y = 8.0000 f = 4 说明 ：实际上，这个初始值所产

39、生的迭代序列kx收敛到方程的根2.925 32 ，而不是收敛到离它最近的根1.400 81,即不是局部收敛的，见例 2.6.6. （3）输入程序 y,f=newjushou(1)运行后输出结果恭喜您！此迭代序列收敛,(x)导数值的绝对值y=|d (x)/dx| 和方程 f(x)=0 的函数 f(x) 值 f 如下y = 0.3566 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 9 页，共 23 页f = 1.7126说明： 实际上，这个初始值所产生的迭代序列kx收敛到离它最近的根1.400 81, 即kx是局部收敛的，见例 2.6.6. （4）

40、输入程序： y,f=newjushou(2)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d (x)/dx| 和方程 f(x)=0 的函数 f(x)值y = 1.2593 f = -2.7188 说明： 虽然y =1.259 31，不满足牛顿切线法局部收敛的条件1)(*x. 但是，实际上，这个初始值所产生的迭代序列kx收敛到离它最近的根1.400 81,即kx是局部收敛的 . 这说明1)(0 x是充分条件，非必要条件（见例 2.6.6）. （5）输入程序 y,f=newjushou(5.5)运行后输出结果请注意观察下面显示的(x)的导数值的绝对值y=|d (x)/dx| 和方程

41、 f(x)=0 的函数 f(x)值y = 1.0447e+005 f = 176.6400说明：实际上， 这个初始值所产生的迭代序列kx收敛到 235.619, 但不是此方程的根，比较初始值5.5 与迭代值235.619 ，可见kx不是局部收敛的，见例 2.6.6. （6）输入程序 y,f=newjushou(8) 运行后输出结果恭喜您！此迭代序列收敛,(x)导数值的绝对值y=|d(x)/dx| 和方程 f(x)=0 的函数 f(x)值 f 如下y = 0.4038 f = -2.9452e+003说明： 实际上，这个初始值所产生的迭代序列kx收敛到方程的根6.290 60 ，而不是收敛到离它

42、最近的根9.424 46,即不是局部收敛的，见例 2.6.6. 2.6.3 牛顿切线法的 MATLAB 程序牛顿切线法求方程0)(xf的近似根kx（精度为tol）需要自行编制程序. 牛顿切线法的MATLAB 主程序输入的量： 初始值0 x、近似根kx的误差限tol，近似根kx的函数值)(kxf的误差限ftol，迭代次数的最大值gxmax 、函数 fnq(x)(xf及其导数dfnq(x)( xf. 输出的量：迭代序列kx、迭代k次得到的迭代值kx (迭代次数超过gxmax时，运行后输出信息 请注意：迭代次数超过给定的最大值gxmax) 、 相邻两次迭代的偏差piancha=1kkxx和它的偏差的

43、相对误差xdpiancha=kkkxxx1的值 .根据式（ 2. 12）和已知条件，现提供名为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); 精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 10 页，共 23 页xdpiancha= piancha/( abs(x(i+

44、1)+eps); i=i+1; xk=x(i);yk=fnq(x(i); (i-1) xk yk piancha xdpiancha if (abs(yk)ftol)&(pianchatol)|(xdpianchagxmax disp( 请注意：迭代次数超过给定的最大值gxmax 。) return ; end (i-1),xk,yk,piancha,xdpiancha;例 2.6.3用牛顿切线法求方程013223xx在9. 04 .00和x的近似根， 要求精度310.然后判断此方法分别在方程的根15.0*和x处的收敛速度 . 解（1）先求方程的近似根.在 MATLAB 工作窗口输入程序 k,

45、xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001, 0.001,100) 运行后输出初始值4.00 x的迭代结果如表 2-10 所示 . 表 2-10 k xkyk1kkxxkkkxxx11 -0.516 7 -0.076 7 0.116 7 0.225 8 2 -0.500 4 -0.001 6 0.016 3 0.032 6 3 -0.500 0 -0.000 0 0.000 4 0.000 7 迭代次数k = 3,根的近似值xk=-0.500, 它的函数值yk=-1.759e-013, 偏差piancha=1.712e-007 和相对误差xdpian

46、cha =3.423e-007. 如果将步骤命令中的-0.4 改作 0.9，则运行后输出初始值4 .00 x的迭代结果为k =7；piancha =7.290e-004；xdpiancha =7.295e-004；xk = 0.999；yk =1.590e-006. （2）再讨论收敛速度比较初始值分别是9 .04 .00和x的两组结果，在4. 00 x处迭代3 次就得到单根5. 0*x的近似值xk=-0.500；而在9 .00 x处迭代7 次才得到二重根1*x的近似值xk=0.999.可见， 牛顿切线迭代法在单根处的迭代速度比二重根处的迭代速度快很多.这正如前面讨论的结果一样，即若*x是0)(

47、xf的单根，则牛顿切线法是平方收敛的；若*x是0)(xf的二重根， 则牛顿切线法是线性收敛的；我们也可以用阶的定义判断它们的收敛速度.留给读者证明 . 2.6.4 求nc的方法及其 MATLAB 程序用c的n次根nc的迭代公式求nc的近似值kx（精度为tol）需要自行编制程序. 求c的n次根nc（当n是偶数时，0c）的 MATLAB 主程序输入的量：初始值0 x、根指数n、被开方数c、nc的误差限tol，迭代次数的最大值 gxmax . 输出的量：迭代序列kx、迭代k次得到的迭代值kx ( 迭代次数达到最大值gxmax时，运行后输出信息 请注意：迭代次数超过给定的最大值gxmax) 、相邻两次

48、迭代的偏差piancha =|1kkxx| 和它的偏差的相对误差xdpiancha=kkkxxx1的值 .精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 11 页，共 23 页根据求c的n次根nc的迭代公式 (2.25) 和已知条件， 现提供名为kai2fang.m 的 M文件：function k,xk,yk,piancha,xdpiancha,P=kainfang(x0,c,n,tol, gxmax) x(1)=x0; for i=1: gxmax u(i)= (x(i)n-c)/(n*x(i)(n-1); x(i+1)= x(i)-u(i);

49、 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 (pianchatol)|(xdpianchagxmax disp( 请注意：迭代次数超过给定的最大值gxmax. ) (i-1),xk,yk,piancha,xdpiancha return ; end P=(i-1),xk,yk,piancha,xdpiancha;例 2.6.5求113，要求精度为610. 解 本题介绍四种解法.方法 1 用求c

50、的n次根nc（当n是偶数时，0c）的 MATLAB 程序计算 . 取初始值100 x，根指数2n，被开方数113c，近似根的精度tol=510，迭代的最大次数gxmax=100. 在工作区间输入程序 k,xk,yk,piancha,xdpiancha=kainfang(10,113,2,1e-5,100)运行后输出结果k=4,piancha=1.610800381968147e-011, xdpiancha=1.515313534117706e-012 xk =10.63014581273465,yk =1.710910332060509e+009可见 ,11310.630 15 ，满足精度5