第二章 非线性方程（组）的数值解法的 MATLAB 程序 9 高等教育出版社 教育电子音像出版社 作者：任玉杰 本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等） ，求根的方法（二分法，迭代法，牛顿法，割线法，米勒（Mller）法和迭代法的加速等）及其 MATLAB 程序，求解非线性方程组的方法及其 MATLAB 程序. 2.1 2.1 2.1 2.1 方程方程方程方程（组组组组）的根及其的根及其的根及其的根及其 MATLABMATLABMATLABMATLAB 命令命令命令命令 2.1.2 2.1.2 2.1.2 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(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 2.1.3 2.1.3 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（按降幂排列） ，运行后输出 dfa 为多项式)(xf的导数)(xf的系数. dfx=poly2sym(dfa) 输入多项式)(xf的导数)(xf的系数 dfa（按降幂排列） ，运行后输出 dfx 为多项式)(xf的导数)(xf. 2.1.4 2.1.4 2.1.4 2.1.4 求解方程求解方程求解方程求解方程（组组组组）的的的的fsolve命令命令命令命令 如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots命令.如果非线性方程（组）是含有超越函数，则无法使用roots命令，需要调用MATLAB系统中提供的另一个程序fsolve来求解.当然，程序fsolve也可以用于多项式方程（组） ，但是它的计算量明显比roots命令的大. fsolve命令使用最小二乘法（least squares method）解非线性方程（组） (FX=)0 的数值解，其中 X 和 F（X）可以是向量或矩阵.此种方法需要尝试着输入解 X 的初始值（向量或矩阵） X0， 即使程序中的迭代序列收敛， 也不一定收敛到(FX=)0的根 （见例 2.1.8） . fsolve 的调用格式的调用格式的调用格式的调用格式： X=fsolve(F,X0) 第二章第二章第二章第二章 非线性方程非线性方程非线性方程非线性方程（组组组组）的数值解法的数值解法的数值解法的数值解法 第二章 非线性方程（组）的数值解法的 MATLAB 程序 10 高等教育出版社 教育电子音像出版社 作者：任玉杰 输入函数)(xF的M文件名和解X的初始值（向量或矩阵）X0，尝试着解方程（组）(FX=)0，运行后输出(FX=)0解的估计值（向量或矩阵）X. 要了解更多的调用格式和功能请输入:help fsolve,查看说明. 2.2 2.2 2.2 2.2 搜索根的方法及其搜索根的方法及其搜索根的方法及其搜索根的方法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）.搜索根的近似位置的常用方法有三种:作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理. 2.2.12.2.12.2.12.2.1 作图法及其作图法及其作图法及其作图法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等.下面介绍两种作图程序. 作函数作函数作函数作函数)(xfy = = = =在区间在区间在区间在区间 a,ba,ba,ba,b 的图形的的图形的的图形的的图形的 MATLABMATLABMATLABMATLAB 程序一程序一程序一程序一 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,ba,ba,ba,b 上的图形的上的图形的上的图形的上的图形的 MATLABMATLABMATLABMATLAB 程序二程序二程序二程序二 将)(xfy =化为)()(xgxh=，其中)()(xgxh和是两个相等的简单函数 x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2) grid,gtext( y1=h(x),y2=g(x) 说明： 此程序在 MATLAB 的工作区输入， 运行后即可出现函数)()(21xgyxhy=和的图形.两图形交点的横坐标即为所要求的根的近似值. 下面举例说明如何用计算机软件MATLAB的图形功能作图. 2.2.2 2.2.2 2.2.2 2.2.2 逐步搜索法逐步搜索法逐步搜索法逐步搜索法及其及其及其及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 逐步搜索法也称试算法.它是求方程0)(=xf根的近似值位置的一种常用的方法. 逐步搜索法依赖于寻找连续函数)(xf满足)(af与)(bf异号的区间,ba.一旦找到区间，无论区间多大，通过某种方法总会找到一个根. MATLAB 的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的收敛判定准则编写其主程序，命名为 zhubuss.m. 逐步搜索法的逐步搜索法的逐步搜索法的逐步搜索法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入区间端点 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); 第二章 非线性方程（组）的数值解法的 MATLAB 程序 11 高等教育出版社 教育电子音像出版社 作者：任玉杰 for k=2:n X(k)=a+k*h;Y(k)=funs(X(k); %程序中调用的funs.m为函数 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.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比较，将会发现逐步搜索法的MATLAB程序的优点.如果精度要求比较高，用这种逐步搜索法是不合算的. 2.3 2.3 2.3 2.3 二分法及其二分法及其二分法及其二分法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 2.3.1 2.3.1 2.3.1 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 2.3.2 2.3.2 2.3.2 二分法的二分法的二分法的二分法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 二分法需自行编制程序， 现提供用二分法求方程 f(x)=0 的根*x的近似值kx的步骤和式（2.3a）编写一个名为 erfen.m 的二分法的 MATLAB 主程序如下. 二分法的二分法的二分法的二分法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 求解方程0)(=xf在开区间（a,b）内的一个根的前提条件是)(xf在闭区间a,b上连续，且0)()(0, disp(注意：ya*yb0,请重新调整区间端点a和b.), return end max1=-1+ceil(log(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-7.从图像可以看出，此曲线有两个驻点33都在 x轴的上方，在(-2,-1)内曲线与 x 轴只有一个交点， 则该方程有唯一一个实根， 且在 (-2,-1)内. 方法方法方法方法 2 2 2 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( k,x,wuca,yx=erfen (-2,-1,0.001) 运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为 图 2-7 第二章 非线性方程（组）的数值解法的 MATLAB 程序 13 高等教育出版社 教育电子音像出版社 作者：任玉杰 k = 9，x=-1.7959， wuca = 9.7656e-004，yx = 0.0037 2.4 2.4 2.4 2.4 迭代法及其迭代法及其迭代法及其迭代法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 确定根的近似位置以后， 接下来的工作就是将根精确化，即按某种方法将初始近似值逐步精确化，直到满足所要求的精确度为止. 上节介绍的二分法是将根精确化的方法之一，但是它的收敛速度较慢，且不能求出偶重根.迭代法可以克服这种缺陷.迭代法是求解方程的根、线性和非线性方程组的解的基本而重要的方法. 2.4.2 2.4.2 2.4.2 2.4.2 迭代法的迭代法的迭代法的迭代法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 1 1 1 1 迭代法需要自行编制程序.下面提供的迭代法的MATLAB程序1使用时只需输入迭代初始值0x、迭代次数 k、迭代公式 xk+1=(xk)和一条命令，运行后就可以输出求迭代序列kx、迭代 k 次得到的迭代值kx和相邻两次迭代的偏差 piancha =|1kkxx| （简称偏差偏差偏差偏差）和偏差的相对误差 xdpiancha=|1kkxx|kx的值，并且具有警报功能（若迭代序列发散，则提示用户“请重新输入新的迭代公式”；若迭代序列收敛，则屏幕会出现“祝贺您！此迭代序列收敛，且收敛速度较快”）.我们可以用这个程序来判断迭代序列的敛散性，也可以用于比较由一个方程得到的几个迭代公式的敛散性的优劣. 迭代法的迭代法的迭代法的迭代法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 1 1 1 1 输入的量：初始值0x、迭代次数 k 和迭代公式), 2 , 1 , 0()(1L=+kxxkk； 运行后输出的量：迭代序列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)=fun1(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-1) piancha xdpiancha xk; 例例例例 2.4.12.4.12.4.12.4.1 求方程102)(2+=xxxf的一个正根. 解解解解 在 MATLAB 工作窗口输入程序 k,piancha,xdpiancha,xk= diedai1(2,5) 运行后输出用迭代公式2/ )10(21kkxx=+的结果 第二章 非线性方程（组）的数值解法的 MATLAB 程序 14 高等教育出版社 教育电子音像出版社 作者：任玉杰 k,piancha,xdpiancha,xk= 1.00000000000000 1.00000000000000 0.33333333333333 3.00000000000000 2.00000000000000 2.50000000000000 5.00000000000000 0.50000000000000 3.00000000000000 4.37500000000000 0.89743589743590 4.87500000000000 4.00000000000000 11.75781250000000 1.70828603859251 -6.88281250000000 5.00000000000000 11.80374145507813 0.63167031671297 -18.68655395507813 请用户注意：此迭代序列发散，请重新输入新的迭代公式 k=5，piancha = 11.80374145507813，xdpiancha = 0.63167031671297， xk = -18.68655395507813 由以上运行后输出的迭代序列与根*x=2.316 624 790 355 40 相差越来越大，即迭代序列kx发散，此迭代法就失败.这时偏差 piancha 逐渐增大且偏差的相对误差 xdpiancha 的值大于 0.5. 用迭代公式)2/(101+=+kkxx运行后输出的结果 k,piancha,xdpiancha,xk= 1.00000000000000 0.50000000000000 0.20000000000000 2.50000000000000 2.00000000000000 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,xdpiancha = 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(21+=+kkkkkxxxxx运行后输出的结果 k,piancha,xdpiancha,xk= 1.00000000000000 0.33333333333333 0.14285714285714 2.33333333333333 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; y = 2.31662479035540. 可见，偏差 piancha 和偏差的相对误差 xdpiancha 的值越来越小，且第 5 次的迭代值xk=2.316 624 790 355 40 与根*x=2.316 624 790 355 40 相差无几，则迭代序列kx收敛，且收敛速度很快，此迭代法成功. 2.4.5 2.4.5 2.4.5 2.4.5 迭代法的迭代法的迭代法的迭代法的 M M M MATLABATLABATLABATLAB 程序程序程序程序 2 2 2 2 迭代法的迭代法的迭代法的迭代法的 MATLAMATLAMATLAMATLAB B B B 程序程序程序程序 2 2 2 2 输入的量：初初始值0x、精度 tol 和和迭代公式), 2 , 1 , 0()(1L=+kxxkk； 运行后输出的量：迭代次数 k，满足精度 tol 的根的近似根kx的迭代序列kx，相邻 两 次 迭 代 的 偏 差pianch=|1kkxx| 和 它 的 偏 差 的 相 对 误 差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 xk yk if (pianchatol)|(xdpianchaddmax disp(迭代次数超过给定的最大值ddmax) return; end P=(i-1),piancha,xdpiancha,xk,yk; 例例例例2.4.5 2.4.5 2.4.5 2.4.5 求x5-3x+1=0在0.3附近的根，精确到4位小数. 解解解解 构造迭代公式构造迭代公式构造迭代公式构造迭代公式)(kkxx=. 改写原方程为等价方程3/) 1(5+=xx.这时3/ ) 1()(5+=xx，初始值为x0=0.5, 迭代公式 3/ ) 1(51+=+kkxx ), 2, 1, 0(L=k. 利用利用利用利用迭代法的迭代法的迭代法的迭代法的MATLABMATLABMATLABMATLAB程序程序程序程序2 2 2 2计算精确到计算精确到计算精确到计算精确到4 4 4 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.206089525390697e-005; xdpiancha =3.603129477781680e-005; xk =0.3347; yk =0.3347. 2.5 2.5 2.5 2.5 迭代过程的加迭代过程的加迭代过程的加迭代过程的加速方法及其速方法及其速方法及其速方法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的问题. 2.5.2 2.5.2 2.5.2 2.5.2 加权迭代法的加权迭代法的加权迭代法的加权迭代法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 加权迭代法的加权迭代法的加权迭代法的加权迭代法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 已知初始值0x、 精度 tol 和迭代公式 xk+1=(xk)， 求满足精度 tol 的根的近似根kx和它的函数值yk及迭代次数 k. 根据加权迭代的公式（2. . . .10）和已知条件，现提供名为 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 第二章 非线性方程（组）的数值解法的 MATLAB 程序 16 高等教育出版社 教育电子音像出版社 作者：任玉杰 disp(迭代次数超过给定的最大值ddmax) return; end P=(i-1),xk,yk; 例例例例 2.5.1 2.5.1 2.5.1 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.85260550406236 3.00000000000000 0.85260550207699 0.85260550208255 k =3； xk =0.852606； yk = 0.852606. 2.5.4 2.5.4 2.5.4 2.5.4 艾特肯艾特肯艾特肯艾特肯(Aitken)(Aitken)(Aitken)(Aitken)加速方法的加速方法的加速方法的加速方法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 艾特肯加速方法的艾特肯加速方法的艾特肯加速方法的艾特肯加速方法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 已知初始值0x、 精度 tol 和迭代公式 xk+1=(xk)， 求满足精度 tol 的根的近似根kx和它的函数值yk及迭代次数 k，p=k,x1,x2,x. 根据艾特肯加速方法的公式（2. . . .11）和已知条件，现提供名为 Aitken.m 的 M 文件如下： function k,xk,yk,p= Aitken (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; end m=0,1:i-1; p=m,x1,x2,x; 例例例例 2.5.3 2.5.3 2.5.3 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.85260436491811 0.85260647150826 0.85260550201407 3.00000000000000 0.85260550201343 0.85260550201398 0.85260550201373 2.6 2.6 2.6 2.6 牛顿牛顿牛顿牛顿(Newton)(Newton)(Newton)(Newton)切线切线切线切线法及其法及其法及其法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 第二章 非线性方程（组）的数值解法的 MATLAB 程序 17 高等教育出版社 教育电子音像出版社 作者：任玉杰 2.6.2 2.6.2 2.6.2 2.6.2 牛顿切线法的收敛性及其牛顿切线法的收敛性及其牛顿切线法的收敛性及其牛顿切线法的收敛性及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 判别牛顿切线法的局部收敛性的判别牛顿切线法的局部收敛性的判别牛顿切线法的局部收敛性的判别牛顿切线法的局部收敛性的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 根据牛顿切线法局部收敛的条件 . 1)()()()(2* =xfxfxfx 输入的量： x 是方程0)(=xf根*x或满足1*0 xx的初始值0x， 函数)(xf及其一阶导数)( xf和二阶导数)( xf； 输出的量：)(*x或)(0x的值.如果)(*x1或1)(0x时，运行后输出信息恭喜您！此迭代序列收敛,(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)*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 2 2 2）输入程序 y,f=newjushou(0) 运行后输出结果 请注意观察下面显示的(x)的导数值的绝对值 y=|d(x)/dx|和方程 f(x)=0 的函数 f(x)值 y = 8.0000 f = 4 说明说明说明说明：实际上，这个初始值所产生的迭代序列 kx收敛到方程的根 2.925 32，而不是收敛到离它最近的根 1.400 81,即不是局部收敛的，见例 2.6.6. （3 3 3 3）输入程序 y,f=newjushou(1) 第二章 非线性方程（组）的数值解法的 MATLAB 程序 18 高等教育出版社 教育电子音像出版社 作者：任玉杰 运行后输出结果 恭喜您！此迭代序列收敛,(x)导数值的绝对值 y=|d(x)/dx|和方程 f(x)=0 的函数 f(x)值 f 如下 y = 0.3566 f = 1.7126 说明说明说明说明： 实际上， 这个初始值所产生的迭代序列 kx收敛到离它最近的根1.400 81,即 kx是局部收敛的，见例 2.6.6. （4 4 4 4）输入程序： 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 y,f=newjushou(5.5) 运行后输出结果 请注意观察下面显示的(x)的导数值的绝对值 y=|d(x)/dx|和方程 f(x)=0 的函数 f(x)值 y = 1.0447e+005 f = 176.6400 说明说明说明说明： 实际上， 这个初始值所产生的迭代序列 kx收敛到 235.619,但不是此方程的根，比较初始值 5.5 与迭代值 235.619，可见 kx不是局部收敛的，见例 2.6.6. （6 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，而不是收敛到离它最近的根 9.424 46,即不是局部收敛的，见例 2.6.6. 2.6.3 2.6.3 2.6.3 2.6.3 牛顿切线牛顿切线牛顿切线牛顿切线法的法的法的法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 牛顿切线法求方程0)(=xf的近似根kx（精度为 tol）需要自行编制程序. 牛顿切线牛顿切线牛顿切线牛顿切线法的法的法的法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值0x、近似根kx的误差限tol，近似根kx的函数值)(kxf的误差限ftol，迭代次数的最大值gxmax、函数fnq(x)(xf=及其导数dfnq(x)( xf=. 输出的量：迭代序列kx、迭代 k 次得到的迭代值kx (迭代次数超过gxmax时，运行后输出信息请注意： 迭代次数超过给定的最大值 gxmax)、 相邻两次迭代的偏差piancha =1kkxx和它的偏差的相对误差xdpiancha=kkkxxx1的值. 第二章 非线性方程（组）的数值解法的 MATLAB 程序 19 高等教育出版社 教育电子音像出版社 作者：任玉杰 根据式（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); xdpiancha= piancha/( abs(x(i+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.32.6.32.6.32.6.3 用牛顿切线法求方程013223=+ xx在9 . 04 . 00和=x的近似根，要求精度310=.然后判断此方法分别在方程的根15 . 0*和=x处的收敛速度. 解解解解 （1 1 1 1）先求方程的近似根先求方程的近似根先求方程的近似根先求方程的近似根. . . . 在 MATLAB 工作窗口输入程序 k,xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001, 0.001,100) 运行后输出初始值4 . 00=x的迭代结果如表 2-10 所示. 表 2-10 k xk yk 1kkxx kkkxxx1 1 -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 和相对误差 xdpiancha =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)(=xf的单根，则牛顿切线法是平方收敛的；若*x是0)(=xf的二重根， 则牛顿切线法是线性收敛的； 我们也可以用阶的定义判断它们的收敛速度.留给读者证明. 2.6.4 2.6.4 2.6.4 2.6.4 求求求求nc 的方法及其的方法及其的方法及其的方法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 用c的n次根nc的迭代公式求nc的近似值kx（精度为tol）需要自行编制程序. 求求求求c的的的的n次根次根次根次根nc（当当当当 n n n n 是偶数时是偶数时是偶数时是偶数时，0c）的的的的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值0x、根指数n、被开方数c、nc的误差限tol，迭代次数的最大 第二章 非线性方程（组）的数值解法的 MATLAB 程序 20 高等教育出版社 教育电子音像出版社 作者：任玉杰 值gxmax. 输出的量：迭代序列kx、迭代 k 次得到的迭代值kx (迭代次数达到最大值gxmax时，运行后输出信息请注意：迭代次数超过给定的最大值 gxmax)、相邻两次迭代的偏差piancha =|1kkxx|和它的偏差的相对误差xdpiancha=kkkxxx1的值. 根据求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); 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.52.6.52.6.52.6.5 求113，要求精度为610. 解解解解 本题介绍四种解法. 方法方法方法方法 1 1 1 1 用求用求用求用求c的的的的n次根次根次根次根nc（当当当当 n n n n 是偶数时是偶数时是偶数时是偶数时，0c）的的的的 MATLABMATLABMATLABMATLAB 程序计算程序计算程序计算程序计算. . . . 取初始值100=x，根指数2=n，被开方数113=c，近似根的精度 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，满足精度510. 方法方法方法方法 2 2 2 2 用牛顿迭代公式用牛顿迭代公式用牛顿迭代公式用牛顿迭代公式（2.122.122.122.12）计算计算计算计算. . . . 设113=x， 则01132=x， 记113)(2= xxf，xxf2)(=.由牛顿迭代公式得，kkkkxxxx211321=+，即）（kkkxxx113211+=+)321 , 0(L，其中=k 取初始值100=x，计算结果列入表 2-12. 表 2-12 因 为 ， 迭 代 次 数k=4时 ， 偏 差34xx 810 k,xk,yk,piancha,xdpiancha=newtonqx(10, 1e-5, 1e-5,100) 运行后，将输出的结果列入下表 2-13. 迭代k=4次，得到精度为510的结果113 10.630 15. 表 2-13 k piancha xdpiancha xk yk 1 0.650 000 0.061 033 10.650 000 0.422 500 2 0.019 836 0.001 866 10.630 164 0.000 393 3 0.000 019 0.000 002 10.630 146 0.000 000 4 0.000 000 0.000 000 10.630 146 0.000 000 方法方法方法方法 4 4 4 4 在 MATLAB 工作空间输入程序 113(1/2) 运行后输出 ans = 10.63014581273465 经过四舍五入后，得到精度为510的结果11310.630 15. 2.6.6 2.6.6 2.6.6 2.6.6 牛顿切线牛顿切线牛顿切线牛顿切线法的加速及其两种法的加速及其两种法的加速及其两种法的加速及其两种 MATLABMATLABMATLABMATLAB 程序程序程序程序 当*x是方程0)(=xf的)2(mm重根时，由牛顿切线公式（2. . . .12）产生的迭代序列 kx收敛到*x的速度较慢，是线性收敛的.我们要提高收敛速度，可以有选择地采用已知根的重数和未知根的重数两种求重根的修正牛顿切线公式. （一一一一） 已知方程根的重数已知方程根的重数已知方程根的重数已知方程根的重数m 已知方程已知方程已知方程已知方程0)(= = = =xf根的重数根的重数根的重数根的重数m, , , ,求重根的修正牛顿切线法的求重根的修正牛顿切线法的求重根的修正牛顿切线法的求重根的修正牛顿切线法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值0x，方程0)(=xf根*x的重数m，近似根kx的误差限tol，近似根kx的函数值)(kxf误差限ftol， 迭代次数的最大值gxmax、 函数fnq(x)(xf=及其导数dfnq(x)( xf=； 输出的量：迭代序列kx、迭代 k 次得到的迭代值kx(迭代次数达到最大值gxmax时运行后输出信息请注意：迭代次数超过给定的最大值 gxmax)、相邻两次迭代的偏差piancha =|1kkxx|和它的偏差的相对误差xdpiancha=kkkxxx1的值. 根据式（2.26）和已知条件，现提供名为 newtonxz.m 的 M 文件： function k,piancha,xdpiancha,xk,yk=newtonxz(m,x0,tol,ftol,gxmax) x(1)=x0; for i=1: gxmax x(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 (pianchatol)|(xdpiancha tol)&(abs(yk)gxmax 第二章 非线性方程（组）的数值解法的 MATLAB 程序 22 高等教育出版社 教育电子音像出版社 作者：任玉杰 disp(请注意：迭代次数超过给定的最大值gxmax.) return; end P=(i-1),piancha,xdpiancha,xk,yk; 例例例例 2.6.7 2.6.7 2.6.7 2.6.7 判断0*=x是方程2)(=xfe026923=xxx的几重根？在区间 1 , 0上，分别用牛顿切线法和求重根的修正牛顿切线公式求此根的近似值kx，使精确到410=. 解解解解 经判断知，0*=x是方程0)(=xf的三重根.在MATLAB工作窗口输入程序 k,piancha,xdpiancha,xk,yk= newtonqx (0.5,1e-4, 1e-4,100) 运行后整理结果列入表 2-20. 表 2-20 k piancha xdpiancha xk yk 1 0.144 10 0.404 89 0.355 90 0.541 98 2 0.107 41 0.432 25 0.248 49 0.168 20 3 0.077 45 0.452 80 0.171 04 0.051 46 4 0.054 50 0.467 60 0.116 55 0.015 58 5 0.037 69 0.477 98 0.078 85 0.004 69 6 0.025 76 0.485 14 0.053 10 0.001 40 7 0.017 46 0.490 01 0.035 63 0.000 42 8 0.011 77 0.493 30 0.023 86 0.000 12 9 0.007 91 0.495 52 0.015 96 0.000 04 10 0.005 30 0.497 00 0.010 66 0.000 01 11 0.003 54 0.498 00 0.007 12 0.000 00 12 0.002 37 0.498 66 0.004 75 0.000 00 13 0.001 58 0.499 11 0.003 17 0.000 00 14 0.001 05 0.499 40 0.002 11 0.000 00 15 0.000 70 0.499 58 0.001 41 0.000 00 16 0.000 47 0.499 69 0.000 94 0.000 00 17 0.000 31 0.499 73 0.000 63 0.000 00 18 0.000 21 0.499 67 0.000 42 0.000 00 19 0.000 14 0.499 44 0.000 28 0.000 00 20 0.000 09 0.498 86 0.000 19 0.000 00 迭代次数 k=20，精确到410=的根的近似值是 xk =0.000 2，其函数值是 yk =5.755 8e-011，是一阶（线性）收敛速度. 根据求重根的修正牛顿切线公式，在 MATLAB 工作窗口输入程序 k,piancha,xdpiancha,xk,yk= newtonxz(3,0.5,1e-4, 1e-4,100) 运行后整理结果得表 2-21. 表 2-21 k piancha xdpiancha xk yk 1 0.432 30 6.385 79 0.067 70 0.002 94 2 0.066 54 57.310 77 0.001 16 0.000 00 3 0.001 16 2673.312 01 0.000 00 -0.000 00 4 0.000 13 0.996 70 0.000 13 0.000 00 5 0.000 13 151.535 85 0.000 00 -0.000 00 第二章 非线性方程（组）的数值解法的 MATLAB 程序 23 高等教育出版社 教育电子音像出版社 作者：任玉杰 6 0.000 10 0.991 44 0.000 10 0.000 00 7 0.000 10 88.974 84 0.000 00 -8.881 78e-016 迭代次数 k=7，精确到410=的根的近似值是 xk =1.121 1e-006，piancha = 9.975 3e-005，xdpiancha=88.974 8，其函数值是 yk = -8.881 8e-016，是二阶收敛（平方收敛）速度.可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多. （二二二二） 未知方程根的重数未知方程根的重数未知方程根的重数未知方程根的重数 未知方程未知方程未知方程未知方程0)(= = = =xf根的重数为根的重数为根的重数为根的重数为m，求重根的修正牛顿切线法的求重根的修正牛顿切线法的求重根的修正牛顿切线法的求重根的修正牛顿切线法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值0x，近似根kx的误差限tol，近似根kx的函数值)(kxf误差限ftol， 迭代次数的最大值gxmax、 函数fnq(x)(xf=及其一、 二阶导数dfnq(x)( xf=和ddfnq(x)( xf=. 运行后输出的量：迭代序列kx、迭代 k 次得到的迭代值kx (迭代次数达到最大值gxmax时，运行后输出信息请注意：迭代次数超过给定的最大值 gxmax)、相邻两次迭代的偏差piancha = |1kkxx|和它的偏差的相对误差xdpiancha=kkkxxx1的值. 根据式（2.27）和已知条件，现提供名为 newtonxz1.m 的 M 文件 function k,piancha,xdpiancha,xk,yk=newtonxz1(x0,tol,ftol,gxmax) x(1)=x0; for i=1: gxmax u(i)=fnq(x(i)/dfnq(x(i); du(i)=1-fnq(x(i)*ddfnq(x(i)/(dfnq(x(i)2+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)gxmax disp(请注意：迭代次数超过给定的最大值gxmax.) return; end P=(i-1),piancha,xdpiancha,xk,yk; 例例例例 2.6.8 2.6.8 2.6.8 2.6.8 用未知重数的求重根的修正牛顿切线法， 求方程 在区间 1 , 0上的根， 使精确到410=，并且与例 2.6.7 的结果比较. 解解解解 在MATLAB工作窗口输入程序 k,piancha,xdpiancha,xk,yk=newtonxz1(0.5,1e-4, 1e-4,100) 运行后整理结果得表 2-22. 表 2-22 k piancha xdpiancha xk yk 1 0.599 23 6.038 57 -0.099 23 -0.008 18 2 0.096 98 43.054 79 -0.002 25 -0.000 00 3 0.002 25 1590.394 68 -0.000 00 0.000 00 4 0.000 02 0.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.这些结果与例 2.6.7 的结果比较见下表 2-23. 第二章 非线性方程（组）的数值解法的 MATLAB 程序 24 高等教育出版社 教育电子音像出版社 作者：任玉杰 表 2-23 牛顿切线法 已知重数的修正牛顿切线法 未知重数的修正牛顿切线法 迭代次数k 20 7 4 根的近似值xk 1.858 1e-004 1.121 1e-006 -2.603 1e-005 偏差piancha 9.269 4e-005 9.975 3e-005 2.461 5e-005 相对误差xdpiancha 4.988 6e-001 8.897 5e+001 9.456 2e-001 xk的函数值yk 5.755 8e-011 -8.881 8e-016 -1.325 6e-013 收敛速度 线性收敛 平方收敛 平方收敛 可见，未知重数的修正牛顿切线法的迭代次数和偏差最小；已知重数的修正牛顿切线法的迭代次数比牛顿切线法少 13 次，比未知重数的修正牛顿切线法多 3 次，但是根的近似值最接近根 0；用牛顿切线法求得的近似根是三者中与根 0 距离最远的. 2.7 2.7 2.7 2.7 割线法及其割线法及其割线法及其割线法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 割线法又称弦截法弦截法弦截法弦截法、弦位法弦位法弦位法弦位法、弦割法弦割法弦割法弦割法和弦法弦法弦法弦法. 牛顿切线法的突出优点是收敛速度快，但它也有明显的缺点： 公式中含有导数，当)(xf较复杂时，使用不方便.为了避免切线法计算导数的麻烦，我们现在设法利用一些函数值)(kxf，)(1kxf，来回避导数值)( kxf的计算.下面具体介绍割线法. 2.7.2 2.7.2 2.7.2 2.7.2 割线法的割线法的割线法的割线法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 用割线法求方程0)(=xf的近似根kx（精度为tol）需要自行编制程序. 割线法的割线法的割线法的割线法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值x01,x02,要求的近似根kx的误差限tol，kx的函数值)(kxf的误差限 ftol，迭代次数的最大值gxmax及其函数fnq(x)； 输出的量： 迭代序列kx、 迭代 k 次得到的迭代值kx (迭代次数超过最大值gxmax时运行后输出信息请注意：迭代次数超过给定的最大值 gxmax)、相邻两次迭代的偏差piancha =1kkxx和它的偏差的相对误差xdpiancha=kkkxxx1的值. 使用初始值0x和1x，根据式（2.28） ，现提供名为 gexian.m 的 M 文件： function k,piancha,xdpiancha,xk,yk=gexian (x01,x02,tol,ftol,gxmax) x(1)=x01;x(2)=x02; for i=2: gxmax u(i)= fnq(x(i)*(x(i)-x(i-1); v(i)= fnq(x(i)-fnq(x(i-1); x(i+1)=x(i)- u(i)/( v(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-2) piancha xdpiancha xk yk if (abs(yk)ftol)&( piancha tol)|(xdpianchagxmax disp(请注意：迭代次数超过给定的最大值gxmax.) return; end P=(i-2),piancha,xdpiancha,xk,yk; 例例例例2.7.1 2.7.1 2.7.1 2.7.1 用割线法求方程=)(xfe0322= xx的一个实根的近似值kx，使精确到 第二章 非线性方程（组）的数值解法的 MATLAB 程序 25 高等教育出版社 教育电子音像出版社 作者：任玉杰 410=. 解解解解 先用作图法确定初始值 x01和 x02，然后在 MATLAB 工作窗口输入程序 k,piancha,xdpiancha,xk,yk= gexian (-0.4,-0.3,1e-4, 1e-4,100) 运行后输出结果经整理得表 2-24.由此可见，迭代 k =3 次，得到精度为310的根的近似值xk =-0.390 6，其函数值为 yk =3.860 7e-008，xk的相邻偏差为 piancha =3.327 1e-005，其相误差为 xdpiancha =8.516 8e-005. 表 2-24 k piancha xdpiancha xk yk 1 2 3 0.090 09 0.000 59 0.000 03 0.230 95 0.001 51 0.000 09 -0.390 09 -0.390 68 -0.390 65 0.001 81 -0.000 11 0.000 00 2.8 2.8 2.8 2.8 抛物线法及其抛物线法及其抛物线法及其抛物线法及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 2.8.2 2.8.2 2.8.2 2.8.2 抛物线法的抛物线法的抛物线法的抛物线法的 MATLABMATLABMATLABMATLAB 程序程序程序程序 用抛物线法求方程0)(=xf的近似根kx（精度为tol）需要自行编制程序. 抛物线法的抛物线法的抛物线法的抛物线法的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值 x1,x2,x3,要求的近似根kx的误差限tol， kx的函数值)(kxf的误差限 ftol，迭代次数的最大值gxmax及其函数fnq(x)； 输出的量：迭代序列kx、迭代 k 次数得到的迭代值kx、函数值)(kxf、相邻两次迭代的偏差piancha =max1kkxx，2kkxx，21kkxx和它的偏差的相对误差xdpiancha=piancha/kx的值(迭代次数超过最大值gxmax时运行后输出信息请注意：迭代次数超过给定的最大值 gxmax,请重新输入初始值 x0). 使用初始值0x，1x和2x， 根据式按公式 （2.35） ，（2.36） 和 （2.39） ， 现提供名为paowu.m的M文件 function k,piancha,xdpiancha,xk,yk=paowu(x1,x2, x3,tol,ftol,gxmax) X=x1, x2, x3; for i=1: gxmax h0= X(1)- X (3); h1= X (2)- X (3); u0= fnq (X (1)- fnq (X (3) ); u1= fnq (X (2)- fnq (X (3); c= fnq (X (3); fenmu= h02* h1- h12* h0; afenzi= u0* h1- u1* h0; bfenzi= u1*h02-u0*h12; a= afenzi/ fenmu; b= bfenzi/ fenmu; deta=b2-4*a*c; cha=abs(X(3)-X(1),abs(X(3)-X(2),abs(X(2)- X(1); piancha =min(cha); xdpiancha= piancha/( abs (X(3)+eps); if deta0 disp(提示：根的判别式值为负数.),detakf=sqrt(deta); elseif deta=0 disp(提示：根的判别式值为零.),detakf=0; else disp(提示：根的判别式值为正数.),detakf=sqrt(deta); end if b0 detakf=- detakf; 第二章 非线性方程（组）的数值解法的 MATLAB 程序 26 高等教育出版社 教育电子音像出版社 作者：任玉杰 end z=-2*c/(b+ detakf); xk= X(3)+ z;q1= xk - X (1); q2= xk - X (2); q3= xk - X (3); if abs(q2) abs(q1) X1=X(2), X(1), X(3); X= X1;Y= fnq(X); end if abs(q3) abs(q1) X2=X(1), X(3), X(2); X= X2;Y= fnq(X); end X(3)= xk; Y(3)=fnq(X(3); yk= Y(3); i piancha xdpiancha xk yk if (abs(yk)ftol)&( piancha tol)|(xdpianchagxmax disp(请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值x0) return; end P=i,piancha,xdpiancha,xk,yk; 例例例例2.8.12.8.12.8.12.8.1 用抛物线法求方程=)(xfe0322= xx的一个实根的近似值kx，使精确到410=. 解解解解 先用作图法确定初始值 x01，x02 和 x03，然后在 MATLAB 工作窗口输入程序 k,piancha,xdpiancha,xk,yk= paowu (-0.4,-0.3, -0.5,1e-4, 1e-4,100) 运行后输出结果 提示：根的判别式值为正数. ans = Columns 1 through 4 1.00000000000000 0.10000000000000 0.20000000000000 -0.39066350562239 Column 5 -0.00005581900299 提示：根的判别式值为正数. ans = Columns 1 through 4 2.00000000000000 0.00933649437761 0.02389906977038 -0.39064638365394 Column 5 -0.00000000923532 提示：根的判别式值为正数. ans = Columns 1 through 4 3.00000000000000 0.00001712196845 0.00004382983990 -0.39064638082059 Column 5 0.00000000000000 k = 3 piancha = 1.712196845282676e-005 xdpiancha = 4.382983989938760e-005 xk = -0.39064638082059 yk = 2.220446049250313e-016 即，迭代 k =3 次，得到精度为410的根的近似值 xk =-0.390 6，其函数值为 yk =2.220 4e-016，xk的相邻最小偏差为 piancha =1.712e-005，其相对误差为 xdpiancha =4.383 0e-005. 将抛物线法与割线法的迭代结果进行比较，见表 2-25.显然，抛物线法比割线法的迭代 第二章 非线性方程（组）的数值解法的 MATLAB 程序 27 高等教育出版社 教育电子音像出版社 作者：任玉杰 速度快. 表 2-25 迭代次数 抛物线法 割线法 k xk yk xk yk 1.0000 2.0000 3.0000 -0.390 66 -0.390 65 -0.390 65 -0.000 06 -0.000 00 0.000 00 -0.390 09 0.001 51 -0.390 65 0.001 81 -0.390 68 0.000 00 2.9 2.9 2.9 2.9 求解非线性方程组的牛顿法求解非线性方程组的牛顿法求解非线性方程组的牛顿法求解非线性方程组的牛顿法及其及其及其及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 求解非线性方程的牛顿法可以推广到二维非线性方程组的情形， 也可以推广到更高维的情形. 2.9.1 2.9.1 2.9.1 2.9.1 求解求解求解求解二元二元二元二元非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法及其及其及其及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 用牛顿法求二元非线性方程组（2.40）的近似根kkyx,（精度为tol）需要自行编制程序. 解解解解二元二元二元二元非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法的的的的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输入的量：初始值 X =x0,y0,要求的近似根),(kkyx的误差限tol， ),(kkyx的函数值2, 1),(=iyxfkki的误差限 ftol， 迭代次数的最大值ngmax和函数2, 1),(=iyxfi及其一阶偏导数； 输出的量：迭代 ci=k 次数得到的迭代向量值=X),(kkyx 、向量的函数值Y),(),(21kkkkyxfyxf=、Y的 2 范数hanfan、迭代方程组（2.44）的系数行列式D的值、相邻两次迭代向量的 2 范数 danfan=norm)(1kkXX+和它的 2 范数的相对误差xddf=danfan/norm)(1+kX的值.(当迭代次数超过最大值ngmax时运行后输出信息请注意：迭代次数超过给定的最大值 ngmax,请重新输入初始值 x0；当0=D时，运行后输出信息请注意！迭代方程组（2.44）的系数行列式的值等于零.) 使用两个初始值,0x0y，根据公式（2.45） ， （2.46） ， （2.47）和（2.49） ，现提供名为newtonzu2.m的M文件 function ci,D,danfan,xddf,hanfan,Xk,Yk=newtonzu2(X,tol,ftol,ngmax) Y=Z(X); for i=1:ngmax dY=JZ(X);D=det(dY); A1=Y(1),dY(1,2); Y(2),dY(2,2); A=det(A1); B1=dY(1,1), Y(1); dY(2,1), Y(2);B=det(B1);Xk=X-A,B/D; hanfan=norm(Y);danfan=norm(Xk-X); xddf=danfan/(norm(Xk)+eps); X=Xk;Y=Z(X);ci=i; if D=0 ci=i; Xk=X-A,B/D; Yk=Y; ci,D,danfan, xddf, hanfan , X, Y; else disp(请注意！迭代方程组（2.44）的系数行列式的值等于零.) end if (hanfan ftol)&( danfan tol)|( xddf gxmax 第二章 非线性方程（组）的数值解法的 MATLAB 程序 28 高等教育出版社 教育电子音像出版社 作者：任玉杰 disp(请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值x0.) return; end end 例例例例 2.9.12.9.12.9.12.9.1 用迭代公式（2.46）求解方程组 =+. 2,162222yxyx 解解解解 这里介绍两种求解已知方程组的方法. 方法方法方法方法 1 1 1 1 用解用解用解用解二元二元二元二元非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法的的的的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序求解求解求解求解. . . . （1）设16)(221+=yxxf， =)(2xf222 yx， 则 xxyxf2),(1=，yyyxf2),(1=， xxyxf2),(2=， yyyxf2),(2=. （2）作函数)(1xf和)(2xf的图形.输入程序 syms x y F1=x2-y2-2;F2=x2+y2-16; ezplot(F1,-4.5,4.5), hold on ezplot(F2,-4.5,4.5) ,hold off 运 行 后 屏 幕 显 示 图 2-29. 取 两 个 初 始 值, 20=x20=y. （3）建立并保存名为 Z.m 的 M 文件 function F=Z(X) x=X(1);y=X(2); F=zeros(1,2); F(1)= x2+y2-16; F(2)= x2-y2-2; （4）建立并保存名为 JZ.m 的 M 文件 function dF=JZ(X) x=X(1);y=X(2); dF=zeros(2,2); dF(1,1)=2*x; dF(1,2)=2*y; dF(2,1)=2*x; dF(2,2)=-2*y; （5）保存名为 newtonzu2.m 的 M 文件； 图 2-29 （6）在MATLAB工作窗口输入程序: X=2,2; k,D,danfan, xddf, hanfan , Xk, Yk= newtonzu2 (X,1e-4,1e-4,100) （7）运行后输出结果 k = 5 D = -63.49803146638493 danfan = 3.932201289951168e-011 xddf = 9.830503224877920e-012 hanfan = 3.336593498083849e-010 Xk = 2.99999999996068 2.64575131106449 Yk = 1.0e-015 * 0 -0.88817841970013 方法方法方法方法 2 2 2 2 用解矩阵方程的方法求解用解矩阵方程的方法求解用解矩阵方程的方法求解用解矩阵方程的方法求解. . . . 编写如下程序 第二章 非线性方程（组）的数值解法的 MATLAB 程序 29 高等教育出版社 教育电子音像出版社 作者：任玉杰 x0=2;y0=2; gxmax =10;tol=1e-6;x(1)=x0;y(1)=y0; i=1;u=1 1;k(1)=1; while(norm(u)tol*norm(x(i),y(i) A=2*x(i),y(i);x(i),-y(i); b=16-x(i)2-y(i)2,2-x(i)2+y(i)2; u=Ab;x(i+1)=x(i)+u(1); y(i+1)=y(i)+u(2); i=i+1;k(i)=i; if(i gxmax) error(请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值); end end k,x,y （2）运行后输出结果,取初始值=xT)2, 2(，迭代次数 n=6，精度=10-6，7512.645000,3.00066=yx与精确解的误差已不超过. 2.9.2 2.9.2 2.9.2 2.9.2 求解求解求解求解n元元元元非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法及其及其及其及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 解解解解n元元元元非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法的的的的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序 输 入 的量 ：初始值 X0,要 求 的近 似根)(kx的 误 差限tol， )(kx的 函数 值(if)(kxni, 2, 1),L=的 误 差 限ftol， 迭 代 次 数 的 最 大 值gxmax和 函 数(if)(kxni, 2, 1),L=及其一阶偏导数； 输出的量：迭代k次数得到的迭代向量值)(kx、向量的函数值(if)(kx) ni, 2, 1L=及其 2 范数hanfan、迭代方程组 （2.51）的雅可比矩阵（2.52） 的行列式D的值、相邻两次迭代向量的 2 范数 danfan=norm)()()1(kkxx+和它的 2 范数的相对误差xddf= danfan/norm)()1(+kx的值.(当迭代次数超过最大值gxmax时运行后输出信息请注意：迭代次数超过给定的最大值 gxmax,请重新输入初始值；当0=D时，运行后输出信息请注意！迭代方程组（2.51）的系数行列式的值等于零.) 使用两个初始值,0x0y， 根据公式 （2.50） ，（2.52） 和 （2.54） ， 现提供名为newtonzu2.m的M文件 function ci,D,danfan, xddf, hanfan , Xk, Yk= newtonzun (X, tol,ftol,gxmax) Y=Z(X); for i=1:gxmax dY=JZ(X);D=det(dY);Xk=X-(dYY); hanfan=norm(Y);danfan=norm(Xk-X); xddf=danfan/(norm(Xk)+eps);X=Xk;Y=Z(X);ci=i; if D=0 ci=i; Xk=X-(dYY);Yk=Y;ci,D,danfan,xddf,hanfan,X,Y; else disp(请注意！迭代方程组（2.51）的系数行列式的值等于零. ) end if (hanfan ftol)&( danfan tol)|( xddf gxmax disp(请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值.) return; end 第二章 非线性方程（组）的数值解法的 MATLAB 程序 30 高等教育出版社 教育电子音像出版社 作者：任玉杰 例例例例 2.9.2 2.9.2 2.9.2 2.9.2 用迭代公式（2.54）求解方程组 =+, 2,162222yxyx要求精度610=. 解解解解 这里介绍两种求解已知方程组的方法. 方法方法方法方法 1 1 1 1 用解用解用解用解n元元元元非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法非线性方程组的牛顿法的的的的 MATLABMATLABMATLABMATLAB 主程序主程序主程序主程序求解求解求解求解 在MATLAB工作窗口输入程序 X=2,2; k,D,danfan, xddf, hanfan , Xk, Yk= newtonzun (X,1e-6,1e-6,100) 运行后输出结果 k=5，D=-63.49803146638493，danfan=3.932201289951168e-011 xddf=9.830503224877920e-012，hanfan=3.336593498083849e-010 Xk = 3.00000000000000 2.64575131106459 Yk = 1.0e-015 * 0 -0.88817841970013 方法方法方法方法 2 2 2 2 用解矩阵方程的方法求解用解矩阵方程的方法求解用解矩阵方程的方法求解用解矩阵方程的方法求解 在 MATLAB 工作窗口输入程序 x0=2;y0=2; gxmax =10;tol=1e-6;x(1)=x0;y(1)=y0;i=1;u=1 1;k(1)=1; while(norm(u)tol*norm(x(i),y(i) A=2*x(i),y(i);x(i),-y(i); b=16-x(i)2-y(i)2,2-x(i)2+y(i)2; u=Ab; x(i+1)=x(i)+u(1); y(i+1)=y(i)+u(2); i=i+1;k(i)=i; if(i gxmax) error(请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值); end end k,x,y 运行后输出结果 ans = 1.00000000000000 2.00000000000000 2.00000000000000 2.00000000000000 3.25000000000000 2.75000000000000 3.00000000000000 3.00961538461538 2.64772727272727 4.00000000000000 3.00001536003932 2.64575204838080 5.00000000000000 3.00000000003932 2.64575131106469 6.00000000000000 3.00000000000000 2.64575131106459 即，取初始值=xT) 2, 2 (，迭代次数 n=6，精度=10-6，7512.645000,3.00066=yx与精确解的误差已不超过. 2.9.3 2.9.3 2.9.3 2.9.3 求解求解求解求解n元元元元非线性方程组的拟牛顿法非线性方程组的拟牛顿法非线性方程组的拟牛顿法非线性方程组的拟牛顿法及其及其及其及其 MATLABMATLABMATLABMATLAB 程序程序程序程序 用拟牛顿法求n元非线性方程组（2.50）的近似根)(kx （精度为tol） ，需要根据给定的方程组（2.50）和（2.58）分别自行编制名为NN.m 和NF.m 的 M 文件作为调用函数程序，然后根据（2.57）式，编写名为ninewton.m 的主程序 解解解解n元元元元非线性方程组的拟牛顿法非线性方程组的拟牛顿法非线性方程组的拟牛顿法非线性方程组的拟牛顿法的的的的 MATLABMATLABMATLABMATLAB 程序程序程序程序 输入的量：初始值0201, XX，要求的近似根)(kx的误差限tol， )(kx的函数值),()(kixfni, 2, 1L=的 误 差 限ftol， 迭 代 次 数 的 最 大 值gxmax和 函 数),()(kixfni, 2, 1L=及其一阶偏导数. 输出的量： 迭代k次数得到的迭代向量值)(kx、 向量的函数值),()(kixfni, 2, 1L=及其 2 范数hanfan、迭代方程组（2.51）的雅可比矩阵（2.52）的行列式D的值、相邻两 第二章 非线性方程（组）的数值解法的 MATLAB 程序 31 高等教育出版社 教育电子音像出版社 作者：任玉杰 次迭代向量的 2 范数 danfan=norm)()()1(kkxx+和它的 2 范数的相对误差xddf=danfan/(norm)1( +kx的值.(当矩阵（2.58）行列式的值等于零时，输出信息请注意！迭代方程组（2.59）中的矩阵（2.58）行列式的值等于零.；当迭代次数超过最大值gxmax时运行后输出信息请注意：迭代次数超过给定的最大值 gxmax,请重新输入初始值) . function k,hanfan,danfan,xddf,Xk,Yk=ninewton(X01,X02,tol, ftol,gxmax) X=X01;X02;n=length(X); for j=1:gxmax Y=NN(X); A=NF(X);D=det(A); Xk=X-(A+ones(n)*eps)Y); hanfan=norm(Y(2);danfan=norm(Xk(2)-X(2); xddf=danfan/(norm(Xk(2)+eps);j=j+1; X=Xk;Y=Z(X);Yk=Y;k=j;warning off MATLAB:singularMatrix if D=0 Xk=X-(AY);Yk=Y;j=j+1; j-1,D,danfan, xddf, hanfan ,Xk, Yk else disp(请注意！迭代方程组（2.59）中的矩阵（2.58）行列式的值等于零.) end if (hanfan ftol)|( danfan tol)|( xddf gxmax) error(请注意：迭代次数超过给定的最大值gxmax.); end