第二章 非线性方程（组）的数值解法的MATLAB程序
本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（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(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命令
如果f(x)为多项式，则可分别用如下命令求方程f(x)0的根，或求导数f'(x)（见表 2-1）.
表 2-1 求解多项式方程（组）的roots命令 命 令 xk =roots(fa) 功 能 输入多项式f(x)的系数fa（按降幂排列），运行后输出xk为f(x)0的全部根. 输入多项式f(x)的系数fa（按降幂排列），运行后输出dfa为多项式f(x)的导数f(x)的系数. 输入多项式f(x)的导数f(x)的系数dfa（按降幂排列），运行后输出dfx为多项式f(x)的导数f(x). 'dfa=polyder(fa) dfx=poly2sym(dfa) ''2.1.4 求解方程（组）的fsolve命令
如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots命令.如果非线性方程（组）是含有超越函数，则无法使用roots命令，需要调用MATLAB系统中提供的另一个程序fsolve来求解.当然，程序fsolve也可以用于多项式方程（组），但是它的计算量明显比roots命令的大.
fsolve命令使用最小二乘法（least squares method）解非线性方程（组）
F(X)0
的数值解，其中X和F（X）可以是向量或矩阵.此种方法需要尝试着输入解X的初始值（向量或矩阵）X0，即使程序中的迭代序列收敛，也不一定收敛到F(X)0的根（见例2.1.8）.
fsolve的调用格式： X=fsolve(F,X0)
9．
第二章 非线性方程（组）的数值解法的MATLAB程序 输入函数F(x)的M文件名和解X的初始值（向量或矩阵）X0，尝试着解方程（组）F(X)0，运行后输出F(X)0解的估计值（向量或矩阵）X.
要了解更多的调用格式和功能请输入:help fsolve,查看说明.
2.2 搜索根的方法及其MATLAB程序
求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根，需要搜索根所在的区间或确定根的初始近似值（简称初始值）.搜索根的近似位置的常用方法有三种:作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理.
2.2.1 作图法及其MATLAB程序
作函数的图形的方法很多，如用计算机软件的图形功能画图，或用高等数学中应用导数作图，或用初等数学的函数叠加法作图等.下面介绍两种作图程序.
作函数yf(x)在区间 [a,b] 的图形的MATLAB程序一
x=a:h:b; % h是步长 y=f(x); plot(x,y) grid, gtext('y=f(x)')
说明：⑴ 此程序在MATLAB的工作区输入，运行后即可出现函数yf(x)的图形.此图形与x轴交点的横坐标即为所要求的根的近似值.
⑵ 区间[a,b] 的两个端点的距离 b-a 和步长h的绝对值越小，图形越精确. 作函数yf(x)在区间 [a,b]上的图形的MATLAB程序二 将yf(x)化为h(x)g(x)，其中h(x)和g(x)是两个相等的简单函数
x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2)
grid,gtext(' y1=h(x),y2=g(x)') 说明：此程序在MATLAB的工作区输入，运行后即可出现函数y1h(x)和y2g(x)的图形.两图形交点的横坐标即为所要求的根的近似值.
下面举例说明如何用计算机软件MATLAB的图形功能作图.
2.2.2 逐步搜索法及其MATLAB程序
逐步搜索法也称试算法.它是求方程f(x)0根的近似值位置的一种常用的方法. 逐步搜索法依赖于寻找连续函数f(x)满足f(a)与f(b)异号的区间[a,b].一旦找到区间，无论区间多大，通过某种方法总会找到一个根.
MATLAB的库函数中没有逐步搜索法的程序，现提供根据逐步搜索法的计算步骤和它的收敛判定准则编写其主程序，命名为zhubuss.m. 逐步搜索法的MATLAB主程序 输入区间端点a和b的值，步长h和精度tol，运行后输出迭代次数 k=(b-a)/h+1,方程 f(x)0根的近似值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);
10．
第二章 非线性方程（组）的数值解法的MATLAB程序
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)) sin(cos2x3)0在区间[2,2]上的根的近似值，要求精度是0.000 1. 解 在MATLAB工作窗口输入如下程序 >> [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 二分法及其MATLAB程序 2.3.1 二分法 二分法也称逐次分半法.它的基本思想是：先确定方程f(x)0含根的区间（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 二分法的MATLAB程序 二分法需自行编制程序，现提供用二分法求方程f(x)=0的根x的近似值xk的步骤和式（2.3a）编写一个名为erfen.m的二分法的MATLAB主程序如下. 二分法的MATLAB主程序 求解方程f(x)0在开区间（a,b）内的一个根的前提条件是f(x)在闭区间[a,b]上连续，且f(a)f(b)0. 输入的量：a和b是闭区间[a,b]的左、右端点，abtol是预先给定的精度. 运行后输出的量：k是使用二分法的次数.x是方程 在（a,b）内的实根x*的近似值， 11． * 第二章 非线性方程（组）的数值解法的MATLAB程序 其精度是abtol.wuca=|bk-ak|/2是使用k次二分法所得到的小区间的长度的一半，即实根x*的近似值x的绝对精度限,满足wuca≤ abtol.yx=f(xk) ，即方程f(x)=0在实根x*的近似值x处的函数值. 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* yb>0, disp('注意：ya*yb>0,请重新调整区间端点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*yx>0 b=x;yb=yx; else a=x; ya=yx; end if b-a< abtol , return, end end k=max1; x; wuca; yx=fun(x); 例2.3.3 确定方程x3-x+4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001. 解 （1）先用两种方法确定方程x3-x+4=0 的实根的分布情况。 方法1 作图法. 在MATLAB工作窗口输入如下程序 >>x=-4:0.1:4; y=x.^3-x +4; plot(x,y) grid,gtext('y=x^3-x+4') 画出函数f(x)=x3-x+4的图像，如图2-7.从图像可以看出，此曲线有两个驻点3都在x3轴的上方，在(-2,-1)内曲线与x轴只有一个交点，则该方程有唯一一个实根，且在 (-2,-1)内. 方法2 试算法. 图2-7 在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)满足f(2)f(1)0，所以此方程在(-2,-1)内有一个实根. （2） 用二分法的主程序计算. 在MATLAB工作窗口输入程序 >> [k,x,wuca,yx]=erfen (-2,-1,0.001) 运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为 12． 第二章 非线性方程（组）的数值解法的MATLAB程序 k = 9，x=-1.7959， wuca = 9.7656e-004，yx = 0.0037 2.4 迭代法及其MATLAB程序 确定根的近似位置以后，接下来的工作就是将根精确化，即按某种方法将初始近似值逐步精确化，直到满足所要求的精确度为止. 上节介绍的二分法是将根精确化的方法之一，但是它的收敛速度较慢，且不能求出偶重根.迭代法可以克服这种缺陷.迭代法是求解方程的根、线性和非线性方程组的解的基本而重要的方法. 2.4.2 迭代法的MATLAB程序1 迭代法需要自行编制程序.下面提供的迭代法的MATLAB程序1使用时只需输入迭代初始值x0、迭代次数k、迭代公式xk+1=(xk)和一条命令，运行后就可以输出求迭代序列{xk}、迭代k次得到的迭代值xk和相邻两次迭代的偏差piancha =|xkxk1| （简称偏差）和偏差的相对误差xdpiancha=|xkxk1| xk的值，并且具有警报功能（若迭代序列发散， 则提示用户“请重新输入新的迭代公式”；若迭代序列收敛，则屏幕会出现“祝贺您！此迭代序列收敛，且收敛速度较快”）.我们可以用这个程序来判断迭代序列的敛散性，也可以用于比较由一个方程得到的几个迭代公式的敛散性的优劣. 迭代法的MATLAB程序1 输入的量：初始值x0、迭代次数k和迭代公式xk1(xk)(k0,1,2,)； 运行后输出的量：迭代序列{xk}、迭代k次得到的迭代值xk、相邻两次迭代的偏差piancha =|xkxk1|和它的偏差的相对误差xdpiancha=xkxk1xk的值. 根据迭代公式（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)&(xdpiancha>0.5)&(k>3) disp('请用户注意：此迭代序列发散，请重新输入新的迭代公式') return; end if (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝贺您！此迭代序列收敛，且收敛速度较快') return; end p=[(i-1) piancha xdpiancha xk]'; 例2.4.1 求方程f(x)x2x10的一个正根. 解 在MATLAB工作窗口输入程序 >> [k,piancha,xdpiancha,xk]= diedai1(2,5) 运行后输出用迭代公式xk1(10xk)/2的结果 13． 22 第二章 非线性方程（组）的数值解法的MATLAB程序 [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相差越来越大，即迭代序列{xk}发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5. 用迭代公式xk110/(xk2)运行后输出的结果 [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接近，则迭代序列{xk}收敛，但收敛速度较慢，此迭代法较为成功. 用迭代公式xk1xk(xk2xk10)/(2xk2)运行后输出的结果 [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. *2可见，偏差piancha和偏差的相对误差xdpiancha的值越来越小，且第5次的迭代值xk=2.316 624 790 355 40与根x=2.316 624 790 355 40相差无几，则迭代序列{xk}收敛，且收敛速度很快，此迭代法成功. * 2.4.5 迭代法的MATLAB程序2 迭代法的MATLAB程序2 输入的量：初初始值x0、精度tol和和迭代公式xk1(xk)(k0,1,2,)； 运行后输出的量：迭代次数k，满足精度tol的根的近似根xk的迭代序列{xk}，相邻两次迭代的偏差pianch=|xkxk1|和它的偏差的相对误差xdpiancha =xkxk1xk的值，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 (piancha 第二章 非线性方程（组）的数值解法的MATLAB程序 return; end end if i>ddmax disp('迭代次数超过给定的最大值ddmax') return; end P=[(i-1),piancha,xdpiancha,xk,yk]'; 例2.4.5 求x5-3x+1=0在0.3附近的根，精确到4位小数. 解 ⑴ 构造迭代公式xk(xk). 改写原方程为等价方程x(x51)/3.这时(x)(x51)/3，初始值为x0=0.5, 迭代公式 5xk1(xk1)/3 (k0,1,2,). ⑵利用迭代法的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.206089525390697e-005; xdpiancha =3.603129477781680e-005; xk =0.3347; yk =0.3347. 2.5 迭代过程的加速方法及其MATLAB程序 对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度，但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的问题. 2.5.2 加权迭代法的MATLAB程序 加权迭代法的MATLAB主程序 已知初始值x0、精度tol和迭代公式xk+1=(xk)，求满足精度tol的根的近似根xk和它的函数值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 (piancha if i>ddmax 15． 第二章 非线性方程（组）的数值解法的MATLAB程序 disp('迭代次数超过给定的最大值ddmax') return; end P=[(i-1),xk,yk]'; x例2.5.1 用（2.10）式求2ex0在0.85附近的一个近似根，要求精度10解 在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. 6. 2.5.4 艾特肯(Aitken)加速方法的MATLAB程序 艾特肯加速方法的MATLAB程序 已知初始值x0、精度tol和迭代公式xk+1=(xk)，求满足精度tol的根的近似根xk和它的函数值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 (piancha end if i>ddmax disp('迭代次数超过给定的最大值ddmax') return; end m=[0,1:i-1]; p=[m',x1',x2',x']; 例2.5.3 用艾特肯加速方法求2e xx0在0.85附近的一个近似根，要求精度 106. 解 在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 牛顿(Newton)切线法及其MATLAB程序 16． 第二章 非线性方程（组）的数值解法的MATLAB程序 2.6.2 牛顿切线法的收敛性及其MATLAB程序 判别牛顿切线法的局部收敛性的MATLAB主程序 f(x*)f(x*)1. 根据牛顿切线法局部收敛的条件 (x)*2[f(x)]*输入的量：x是方程f(x)0根x或满足x0x*1的初始值x0，函数f(x)及其一阶导数f'(x)和二阶导数f\"(x)； 输出的量：(x)或(x0)的值.如果(x)1或(x0)1时，运行后输出信息'恭喜您！此迭代序列收敛,φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下'；否则，运行后输出信息'请注意观察下面显示的 φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值'. 根据牛顿切线法局部收敛的条件和方程f(x)0的函数f(x)及其一、二阶导数，现提***供程序： 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)导数值的绝对值y=|dφ(x)/dx|和方程 f(x)=0的函数f(x)值f如下') else disp('请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程 f(x)=0的函数f(x)值') end P=[y,f]'; 例 2.6.2 用牛顿切线法的局部收敛性判别方程 esinx4的近似根时，由下列初始值x0产生的迭代序列是否收敛？ ⑴x0-1; ⑵x00; ⑶x01; ⑷x02; ⑸ x05.5；⑹x08. 解 在MATLAB工作窗口输入程序 >> [y,f]=newjushou(-1) 运行后输出结果 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 139.5644 f = 4.3096 x说明：实际上，这个初始值所产生的迭代序列xk发散,见例 2.6.6. （2）输入程序 >> [y,f]=newjushou(0) 运行后输出结果 说明：实际上，这个初始值所产生的迭代序列xk收敛到方程的根2.925 32，而不是收敛到离它最近的根1.400 81,即不是局部收敛的，见例 2.6.6. （3）输入程序 >> [y,f]=newjushou(1) 17． 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 8.0000 f = 4 第二章 非线性方程（组）的数值解法的MATLAB程序 运行后输出结果 恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下 y = 0.3566 f = 1.7126 说明：实际上，这个初始值所产生的迭代序列xk收敛到离它最近的根1.400 81,即xk是局部收敛的，见例 2.6.6. （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，不满足牛顿切线法局部收敛的条件(x)1*.但是，实 是局部收际上，这个初始值所产生的迭代序列xk收敛到离它最近的根1.400 81,即xk敛的.这说明(x0)1是充分条件，非必要条件（见例 2.6.6）. （5）输入程序 >> [y,f]=newjushou(5.5) 运行后输出结果 请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值 y = 1.0447e+005 f = 176.6400 比较初始值5.5与迭代值235.619，可见xk不是局部收敛的，见例 2.6.6. （6）输入程序 >> [y,f]=newjushou(8) 运行后输出结果 说明：实际上，这个初始值所产生的迭代序列xk收敛到235.619,但不是此方程的根， 说明：实际上，这个初始值所产生的迭代序列xk收敛到方程的根6.290 60，而不是收敛到离它最近的根9.424 46,即不是局部收敛的，见例 2.6.6. 恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值f如下 y = 0.4038 f = -2.9452e+003 2.6.3 牛顿切线法的MATLAB程序 牛顿切线法求方程f(x)0的近似根xk（精度为tol）需要自行编制程序. 牛顿切线法的MATLAB主程序 输入的量：初始值x0、近似根xk的误差限tol，近似根xk的函数值f(xk)的误差限ftol，迭代次数的最大值gxmax、函数fnq(x)f(x)及其导数dfnq(x)f'(x). 输出的量：迭代序列{xk}、迭代k次得到的迭代值xk (迭代次数超过gxmax时，运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =xkxk1和它的偏差的相对误差xdpiancha=xkxk1xk的值. 18． 第二章 非线性方程（组）的数值解法的MATLAB程序 根据式（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) disp('请注意：迭代次数超过给定的最大值gxmax。') return; end [(i-1),xk,yk,piancha,xdpiancha]'; 例 2.6.3 用牛顿切线法求方程2x3x10在x00.4和0.9的近似根，要求精度10.然后判断此方法分别在方程的根x0.5和1处的收敛速度. 解 （1）先求方程的近似根. 在MATLAB工作窗口输入程序 >> [k,xk,yk,piancha,xdpiancha]=newtonqx(-0.4,0.001, 0.001,100) 运行后输出初始值x00.4的迭代结果如表 2-10所示. 表 2-10 332*k xk yk xkxk1 xkxk1xk 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，则运行后输出初始值x00.4的迭代结果为 k =7；piancha =7.290e-004；xdpiancha =7.295e-004；xk = 0.999；yk =1.590e-006. （2）再讨论收敛速度 比较初始值分别是x00.4和0.9的两组结果，在x00.4处迭代3次就得到单根 x*0.5的近似值xk =-0.500；而在x00.9处迭代7次才得到二重根x*1的近似值 xk=0.999.可见，牛顿切线迭代法在单根处的迭代速度比二重根处的迭代速度快很多.这正如前面讨论的结果一样，即若x是f(x)0的单根，则牛顿切线法是平方收敛的；若x是 则牛顿切线法是线性收敛的；我们也可以用阶的定义判断它们的收敛速f(x)0的二重根， 度.留给读者证明. ** 2.6.4 求nc的方法及其MATLAB程序 用c的n次根nc的迭代公式求nc的近似值xk（精度为tol）需要自行编制程序. 求c的n次根nc（当n是偶数时，c0）的MATLAB主程序 输入的量：初始值x0、根指数n、被开方数c、nc的误差限tol，迭代次数的最大 19． 第二章 非线性方程（组）的数值解法的MATLAB程序 值gxmax. 输出的量：迭代序列{xk}、迭代k次得到的迭代值xk (迭代次数达到最大值gxmax时，运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =|xkxk1|和它的偏差的相对误差xdpiancha=xkxk1xk的值. 根据求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 (piancha if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax.') [(i-1),xk,yk,piancha,xdpiancha] return; end P=[(i-1),xk,yk,piancha,xdpiancha]'; 例2.6.5 求113，要求精度为10解 本题介绍四种解法. 6. 方法1 用求c的n次根nc（当n是偶数时，c0）的MATLAB程序计算. 5取初始值x010，根指数n2，被开方数c113，近似根的精度tol=10，迭代 的最大次数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，满足精度10. 方法2 用牛顿迭代公式（2.12）计算. 2设x113，则x1130，记f(x)x2113，f'(x)2x.由牛顿迭代公式得， 52xk1131113）(其中k0,1，2，3，) ，即xk1（xkxk1xk2xk2xk取初始值x010，计算结果列入表 2-12. 表 2-12 因为，迭代次数k=4时，偏差偏差xk1xk 迭代次数k 58x4x310，满足精度10，所以， 1 0.650 000 2 0.019 836 11310.630 15. 3 0.000 019 方法3 用牛顿切线法的MATLAB主程序计 4 0.000 000 算. 分别建立名为fnq.m和dfnq.m的M文件 function y=fnq(x) 根的近似值xk 10.650 000 10.630 164 10.630 146 10.630 146 20． 第二章 非线性方程（组）的数值解法的MATLAB程序 y=x^2-113; function y=dfnq(x) y=2*x; 在MATLAB工作窗口输入程序 >> [k,xk,yk,piancha,xdpiancha]=newtonqx(10, 1e-5, 1e-5,100) 5运行后，将输出的结果列入下表 2-13. 迭代k=4次，得到精度为10的结果113 10.630 15. 表 2-13 k 1 2 3 4 piancha 0.650 000 0.019 836 0.000 019 0.000 000 xdpiancha 0.061 033 0.001 866 0.000 002 0.000 000 xk 10.650 000 10.630 164 10.630 146 10.630 146 yk 0.422 500 0.000 393 0.000 000 0.000 000 方法4 在MATLAB工作空间输入程序 >> 113^(1/2) 运行后输出 ans = 10.63014581273465 5经过四舍五入后，得到精度为10的结果11310.630 15. 2.6.6 牛顿切线法的加速及其两种MATLAB程序 当x是方程f(x)0的m(m2)重根时，由牛顿切线公式（2.12）产生的迭代序列xk收敛到x的速度较慢，是线性收敛的.我们要提高收敛速度，可以有选择地采用已 **知根的重数和未知根的重数两种求重根的修正牛顿切线公式. （一） 已知方程根的重数m 已知方程f(x)0根的重数m,求重根的修正牛顿切线法的MATLAB主程序 输入的量：初始值x0，方程f(x)0根x的重数m，近似根xk的误差限tol，近似根xk的函数值f(xk)误差限ftol，迭代次数的最大值gxmax、函数fnq(x)f(x)及其导数dfnq(x)f'(x)； *输出的量：迭代序列{xk}、迭代k次得到的迭代值xk(迭代次数达到最大值gxmax时运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =|xkxk1|和它的偏差的相对误差xdpiancha=xkxk1xk的值. 根据式（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 ((piancha if i>gxmax 21． 第二章 非线性方程（组）的数值解法的MATLAB程序 disp('请注意：迭代次数超过给定的最大值gxmax.') return; end P=[(i-1),piancha,xdpiancha,xk,yk]'; 3x2例2.6.7 判断x0是方程f(x)2e9x6x20的几重根？在区间 *[0,1]上，分别用牛顿切线法和求重根的修正牛顿切线公式求此根的近似值xk，使精确到 104. 解 经判断知，x0是方程f(x)0的三重根.在MATLAB工作窗口输入程序 >> [k,piancha,xdpiancha,xk,yk]= newtonqx (0.5,1e-4, 1e-4,100) 运行后整理结果列入表 2-20. 表 2-20 *k 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 piancha 0.144 10 0.107 41 0.077 45 0.054 50 0.037 69 0.025 76 0.017 46 0.011 77 0.007 91 0.005 30 0.003 54 0.002 37 0.001 58 0.001 05 0.000 70 0.000 47 0.000 31 0.000 21 0.000 14 0.000 09 xdpiancha 0.404 89 0.432 25 0.452 80 0.467 60 0.477 98 0.485 14 0.490 01 0.493 30 0.495 52 0.497 00 0.498 00 0.498 66 0.499 11 0.499 40 0.499 58 0.499 69 0.499 73 0.499 67 0.499 44 0.498 86 xk 0.355 90 0.248 49 0.171 04 0.116 55 0.078 85 0.053 10 0.035 63 0.023 86 0.015 96 0.010 66 0.007 12 0.004 75 0.003 17 0.002 11 0.001 41 0.000 94 0.000 63 0.000 42 0.000 28 0.000 19 yk 0.541 98 0.168 20 0.051 46 0.015 58 0.004 69 0.001 40 0.000 42 0.000 12 0.000 04 0.000 01 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 0.000 00 迭代次数k=20，精确到10的根的近似值是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 4k 1 2 3 4 5 piancha xdpiancha xk 0.432 30 6.385 79 0.067 70 0.066 54 57.310 77 0.001 16 0.001 16 2673.312 01 0.000 00 0.000 13 0.996 70 0.000 13 0.000 13 151.535 85 0.000 00 yk 0.002 94 0.000 00 -0.000 00 0.000 00 -0.000 00 22． 第二章 非线性方程（组）的数值解法的MATLAB程序 6 0.000 10 7 0.000 10 迭代次数k=7，精确到10的根的近似值是xk =1.121 1e-006，piancha = 9.975 3e-005，xdpiancha=88.974 8，其函数值是yk = -8.881 8e-016，是二阶收敛（平方收敛）速度.可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多. （二） 未知方程根的重数 未知方程f(x)0根的重数为m，求重根的修正牛顿切线法的MATLAB主程序 输入的量：初始值x0，近似根xk的误差限tol，近似根xk的函数值f(xk)误差限ftol，迭代次数的最大值gxmax、函数fnq(x)f(x)及其一、二阶导数dfnq(x)f'(x)''和ddfnq(x)f(x). 0.991 44 0.000 10 0.000 00 88.974 84 0.000 00 -8.881 78e-016 4运行后输出的量：迭代序列{xk}、迭代k次得到的迭代值xk (迭代次数达到最大值gxmax时，运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha = |xkxk1|和它的偏差的相对误差xdpiancha=xkxk1xk的值. 根据式（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 ((piancha if i>gxmax disp('请注意：迭代次数超过给定的最大值gxmax.') return; end P=[(i-1),piancha,xdpiancha,xk,yk]'; 例2.6.8 用未知重数的求重根的修正牛顿切线法，求方程 在区间[0,1]上的根，使精确到10，并且与例2.6.7的结果比较. 解 在MATLAB工作窗口输入程序 >> [k,piancha,xdpiancha,xk,yk]=newtonxz1(0.5,1e-4, 1e-4,100) 运行后整理结果得表 2-22. 表 2-22 4xdpiancha 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. 23． k piancha 第二章 非线性方程（组）的数值解法的MATLAB程序 迭代次数k 根的近似值xk 偏差piancha 相对误差xdpiancha xk的函数值yk 收敛速度 牛顿切线法 20 1.858 1e-004 9.269 4e-005 4.988 6e-001 5.755 8e-011 线性收敛 表 2-23 已知重数的修正牛顿切线法 7 1.121 1e-006 9.975 3e-005 8.897 5e+001 -8.881 8e-016 平方收敛 未知重数的修正牛顿切线法 4 -2.603 1e-005 2.461 5e-005 9.456 2e-001 -1.325 6e-013 平方收敛 可见，未知重数的修正牛顿切线法的迭代次数和偏差最小；已知重数的修正牛顿切线法的迭代次数比牛顿切线法少13次，比未知重数的修正牛顿切线法多3次，但是根的近似值最接近根0；用牛顿切线法求得的近似根是三者中与根0距离最远的. 2.7 割线法及其MATLAB程序 割线法又称弦截法、弦位法、弦割法和弦法. 牛顿切线法的突出优点是收敛速度快，但它也有明显的缺点：公式中含有导数，当f(x)较复杂时，使用不方便.为了避免切线法计算导数的麻烦，我们现在设法利用一些函数值f(xk)，f(xk1)，…来回避导数值f'(xk)的计算.下面具体介绍割线法. 2.7.2 割线法的MATLAB程序 用割线法求方程f(x)0的近似根xk（精度为tol）需要自行编制程序. 割线法的MATLAB主程序 输入的量：初始值x01,x02,要求的近似根xk的误差限tol，xk的函数值f(xk)的误差限ftol，迭代次数的最大值gxmax及其函数fnq(x)； 输出的量：迭代序列{xk}、迭代k次得到的迭代值xk (迭代次数超过最大值gxmax时运行后输出信息'请注意：迭代次数超过给定的最大值gxmax')、相邻两次迭代的偏差piancha =xkxk1和它的偏差的相对误差xdpiancha=xkxk1xk的值. 使用初始值x0和x1，根据式（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) disp('请注意：迭代次数超过给定的最大值gxmax.') return; end P=[(i-2),piancha,xdpiancha,xk,yk]'; 例2.7.1 用割线法求方程f(x)e 2x3x20的一个实根的近似值xk，使精确到 24． 第二章 非线性方程（组）的数值解法的MATLAB程序 104. 解 先用作图法确定初始值x01和x02，然后在MATLAB工作窗口输入程序 >> [k,piancha,xdpiancha,xk,yk]= gexian (-0.4,-0.3,1e-4, 1e-4,100) 运行后输出结果经整理得表 2-24.由此可见，迭代k =3次，得到精度为10的根的近似值xk =-0.390 6，其函数值为yk =3.860 7e-008，xk的相邻偏差为piancha =3.327 1e-005，其相误差为xdpiancha =8.516 8e-005. 表 2-24 3k 1 2 3 piancha xdpiancha xk yk 0.090 09 0.230 95 -0.390 09 0.001 81 0.000 59 0.001 51 -0.390 68 -0.000 11 0.000 03 0.000 09 -0.390 65 0.000 00 2.8 抛物线法及其MATLAB程序 2.8.2 抛物线法的MATLAB程序 用抛物线法求方程f(x)0的近似根xk（精度为tol）需要自行编制程序. 抛物线法的MATLAB主程序 输入的量：初始值x1,x2,x3,要求的近似根xk的误差限tol， xk的函数值f(xk)的误差限ftol，迭代次数的最大值gxmax及其函数fnq(x)； 输出的量：迭代序列{xk}、迭代k次数得到的迭代值xk、函数值f(xk)、相邻两次迭代的偏差piancha =max｛xkxk1，xkxk2，xk1xk2｝和它的偏差的相对误差xdpiancha=piancha/xk的值(迭代次数超过最大值gxmax时运行后输出信息'请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值x0'). 使用初始值x0，x1和x2，根据式按公式（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= h0^2* h1- h1^2* h0; afenzi= u0* h1- u1* h0; bfenzi= u1*h0^2-u0*h1^2; a= afenzi/ fenmu; b= bfenzi/ fenmu; deta=b^2-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 deta<0 disp('提示：根的判别式值为负数.'),detakf=sqrt(deta); elseif deta<=0 disp('提示：根的判别式值为零.'),detakf=0; else disp('提示：根的判别式值为正数.'),detakf=sqrt(deta); end if b<0 detakf=- detakf; 25． 第二章 非线性方程（组）的数值解法的MATLAB程序 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) disp('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值x0') return; end P=[i,piancha,xdpiancha,xk,yk]'; 例2.8.1 用抛物线法求方程f(x)e 2x3x20的一个实根的近似值xk，使精确到 104. 解 先用作图法确定初始值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次，得到精度为10的根的近似值xk =-0.390 6，其函数值为yk =2.220 4e-016，xk的相邻最小偏差为piancha =1.712e-005，其相对误差为xdpiancha =4.383 0e-005. 将抛物线法与割线法的迭代结果进行比较，见表 2-25.显然，抛物线法比割线法的迭代 26． 4 第二章 非线性方程（组）的数值解法的MATLAB程序 速度快. 表 2-25 迭代次数 抛物线法 k xk yk 1.0000 -0.390 66 -0.000 06 2.0000 -0.390 65 -0.000 00 3.0000 -0.390 65 0.000 00 割线法 xk -0.390 09 0.001 51 -0.390 65 yk 0.001 81 -0.390 68 0.000 00 2.9 求解非线性方程组的牛顿法及其MATLAB程序 求解非线性方程的牛顿法可以推广到二维非线性方程组的情形，也可以推广到更高维的情形. 2.9.1 求解二元非线性方程组的牛顿法及其MATLAB程序 用牛顿法求二元非线性方程组（2.40）的近似根xk,yk（精度为tol）需要自行编制程 序. 解二元非线性方程组的牛顿法的MATLAB主程序 输入的量：初始值X =[x0,y0],要求的近似根(xk,yk)的误差限tol， (xk,yk)的函数值fi(xk,yk),i1,2的误差限ftol，迭代次数的最大值ngmax和函数fi(x,y),i1,2及其一阶偏导数； 输出的量：迭代ci=k次数得到的迭代向量值X(xk,yk) 、向量的函数值Yf1(xk,yk),f2(xk,yk)、Y的2范数hanfan、迭代方程组（2.44）的系数行列式D的值、相邻两次迭代向量的2范数danfan=norm(Xk1Xk)和它的2 范数的相对误差xddf=danfan/norm(Xk1)的值.(当迭代次数超过最大值ngmax时运行后输出信息’请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值x0’；当D0时，运行后输出信息’请注意！迭代方程组（2.44）的系数行列式的值等于零.’) 使用两个初始值x0,y0，根据公式（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 if i>gxmax 27． 第二章 非线性方程（组）的数值解法的MATLAB程序 disp('请注意：迭代次数超过给定的最大值ngmax,请重新输入初始值x0.') return; end end 22xy16,例2.9.1 用迭代公式（2.46）求解方程组 22xy2.解 这里介绍两种求解已知方程组的方法. 方法1 用解二元非线性方程组的牛顿法的MATLAB主程序求解. （1）设f1(x)xy16， 22f2(x)x2y22， 则 f1(x,y)xyx（2）作函数f1(x)和f2(x)的图形.输入程 2x， f1(x,y)2y， f2(x,y)2x， f2(x,y)2y. y序 >>syms x y F1=x^2-y^2-2;F2=x^2+y^2-16; ezplot(F1,[-4.5,4.5]), hold on ezplot(F2,[-4.5,4.5]) ,hold off 运行后屏幕显示图2-29.取两个初始值x02,y02. （3）建立并保存名为Z.m的M文件 function F=Z(X) x=X(1);y=X(2); F=zeros(1,2); F(1)= x^2+y^2-16; F(2)= x^2-y^2-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 用解矩阵方程的方法求解. 编写如下程序 28． 第二章 非线性方程（组）的数值解法的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=A\\b;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）运行后输出结果,取初始值x(2,2)，迭代次数n=6，精度ε=10-6， Tx63.000000,y62.645751与精确解的误差已不超过ε. 2.9.2 求解n元非线性方程组的牛顿法及其MATLAB程序 解n元非线性方程组的牛顿法的MATLAB主程序 输入的量：初始值X0,要求的近似根x(k)的误差限tol， x(k)的函数值fi(x(k)),i1,2,,n的误差限ftol，迭代次数的最大值gxmax和函数fi(x(k)),i1,2,,n及其一阶偏导数； 输出的量：迭代k次数得到的迭代向量值x(k)(k)、向量的函数值fi(x) i1,2,,n及其2范数hanfan、迭代方程组（2.51）的雅可比矩阵（2.52）的行列式D的值、相邻两次迭代向量的2 范数danfan=norm(x(k1)x(k))和它的2 范数的相对误差(k1))的值.(当迭代次数超过最大值gxmax时运行后输出信息’请注xddf= danfan/norm(x意：迭代次数超过给定的最大值gxmax,请重新输入初始值’；当D0时，运行后输出信息’请注意！迭代方程组（2.51）的系数行列式的值等于零.’) 使用两个初始值x0,y0，根据公式（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-(dY\\Y')'; 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-(dY\\Y')';Yk=Y;[ci,D,danfan,xddf,hanfan,X,Y]; else disp('请注意！迭代方程组（2.51）的系数行列式的值等于零. ') end if (hanfan disp('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值.') return; end 29． 第二章 非线性方程（组）的数值解法的MATLAB程序 22xy16,6例2.9.2 用迭代公式（2.54）求解方程组 要求精度10. 22xy2,解 这里介绍两种求解已知方程组的方法. 方法1 用解n元非线性方程组的牛顿法的MATLAB主程序求解 在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 用解矩阵方程的方法求解 在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=A\\b; 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 T即，取初始值x(2,2)，迭代次数n=6，精度ε=10-6，x63.000000,y62.645751与精确解的误差已不超过ε. 2.9.3 求解n元非线性方程组的拟牛顿法及其MATLAB程序 用拟牛顿法求n元非线性方程组（2.50）的近似根x （精度为tol），需要根据给定的方程组（2.50）和（2.58）分别自行编制名为NN.m和NF.m的M文件作为调用函数程序，然后根据（2.57）式，编写名为ninewton.m的主程序 解n元非线性方程组的拟牛顿法的MATLAB程序 输入的量：初始值X01,X02，要求的近似根x(k)(k)的误差限tol， x(k)的函数值fi(x(k)),i1,2,,n的误差限ftol，迭代次数的最大值gxmax和函数fi(x(k)),i1,2,,n及其一阶偏导数. 输出的量：迭代k次数得到的迭代向量值x(k)(k)、向量的函数值fi(x),i1,2,,n及其2范数hanfan、迭代方程组（2.51）的雅可比矩阵（2.52）的行列式D的值、相邻两 30． 第二章 非线性方程（组）的数值解法的MATLAB程序 次迭代向量的2范数danfan=norm(x(k1)x(k))和它的2范数的相对误差xddf=danfan/norm(x(k1))的值.(当矩阵（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-(A\\Y’)’;Yk=Y;j=j+1; [j-1,D,danfan, xddf, hanfan ,Xk, Yk]’ else disp(‘请注意！迭代方程组（2.59）中的矩阵（2.58）行列式的值等于零.’) end if (hanfan if(j> gxmax) error(‘请注意：迭代次数超过给定的最大值gxmax.’); end 31． 因篇幅问题不能全部显示，请点此查看更多更全内容