第二章 非线性方程(组)的数值解法

7.1 方程（组）的根及其MATLAB 命令7.1.1 求解方程（组）的solve 命令求方程f (x )=q (x )的根可以用MATLAB 命令:>> x=solve('方程f(x)=q(x)','待求符号变量x ')求方程组f i (x 1,…,x n )=q i (x 1,…,x n ) (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)7.1.2 求解方程（组）的fsolve 命令fsolve 的调用格式： X=fsolve(F,X0)7.2 搜索根的方法及其MATLAB 程序7.2.1 作图法及其MATLAB 程序作函数)(x f y =在区间 [a,b ] 的图形的MATLAB 程序一x=a:h:b; % h 是步长 y=f(x); plot(x,y)grid, gtext('y=f(x)')说明：⑴ 此程序在MATLAB 的工作区输入，运行后即可出现函数)(x f y =的图形.此图形与x 轴交点的横坐标即为所要求的根的近似值.⑵ 区间[a,b ] 的两个端点的距离 b-a 和步长h 的绝对值越小，图形越精确. 作函数)(x f y =在区间 [a,b ]上的图形的MATLAB 程序二 将)(x f y =化为)()(x g x h =，其中)()(x g x h 和是两个相等的简单函数.x=a:h:b; y1=h(x); y2=g(x); plot(x, y1, x, y2)grid,gtext(' y1=h(x),y2=g(x)') 说明：此程序在MATLAB 的工作区输入，运行后即可出现函数)()(21x g y x h y ==和的图形.两图形交点的横坐标即为所要求的根的近似值.7.2.2 逐步搜索法及其MATLAB 程序逐步搜索法的MATLAB 主程序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); for k=2:nX(k)=a+k*h;Y(k)=funs(X(k)); %程序中调用的funs.m 为函数if sk<=0,m=m+1;r(m)=X(k); endxielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-1)); if (abs(Y(k))<tol)&( xielv<=0) m=m+1;r(m)=X(k); end end例7.2.1 用逐步搜索法的MATLAB 程序分别求方程0332223=--+x x x 和0)2sin(cos 3=x 在区间]2,2[-上的根的近似值，要求精度是0.000 1.解 将逐步搜索法的MATLAB 程序保存名为zhubuss.m 的M 文件. 建立M 文件funs.mfunction y=funs(x)y=2.*x.^3+2.*x.^2-3.*x-3 在MATLAB 工作窗口输入如下程序>> [k,r]=zhubuss(-2,2,0.001,0.0001)运行后输出的结果k =4001r = -1.2240 -1.0000 -1.0000 -0.9990 1.2250即搜索点的个数为k =4 001，其中有5个是方程0332223=--+x x x 的近似根，即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)) 代替，可得到方程0)2sin(cos 3=x 在区间]2,2[-上的根的近似值如下r = -1.9190 -1.7640 -1.5770 -1.3300 -0.9220 0.92301.3310 1.5780 1.7650 1.92007.3 二分法及其MATLAB 程序7.3.1二分法的MATLAB 程序 二分法的MATLAB 主程序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 endmax1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil 是向∞+ 方向取整for k=1: max1+1a;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==0a=x; b=x; elseif yb*yx>0b=x;yb=yx; elsea=x; ya=yx; endif b-a< abtol , return , end endk=max1; x; wuca; yx=fun(x);例7.3.1 确定方程x 3-x +4=0的实根的分布情况，并用二分法求在开区间 (-2,-1)内的实根的近似值，要求精度为0.001.解 （1）先用两种方法确定方程x 3-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 )=x 3-x+4的图像.从图像可以看出，此曲线有两个驻点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(<-⋅-f f ，所以此方程在(-2,-1)内有一个实根.（2） 用二分法的主程序计算. 在MATLAB 工作窗口输入程序>> [k,x,wuca,yx]=erfen (-2,-1,0.001)运行后屏幕显示用二分法计算过程被列入表 2-3，其余结果为k = 9，x=-1.7959，wuca = 9.7656e-004，yx = 0.00377.4 迭代法及其MATLAB 程序7.4.1迭代法的MATLAB 程序1迭代法的MATLAB 主程序1function [k,piancha,xdpiancha,xk]=diedai1(x0,k) % 输入的量--x0是初始值,k 是迭代次数 x(1)=x0; for i=1:kx(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]endif (piancha >1)&(xdpiancha>0.5)&(k>3)disp('请用户注意：此迭代序列发散，请重新输入新的迭代公式') return ; endif (piancha < 0.001)&(xdpiancha< 0.0000005)&(k>3) disp('祝贺您！此迭代序列收敛，且收敛速度较快') return ; endp=[(i-1) piancha xdpiancha xk]';例7.4.1 求方程102)(2-+=x x x f 的一个正根. 解 在MATLAB 工作窗口输入程序>> [k,piancha,xdpiancha,xk]= diedai1(2,5)运行后输出用迭代公式2/)10(21k k x x -=+的结果[k,piancha,xdpiancha,xk]=1.00000000000000 1.00000000000000 0.33333333333333 3.000000000000002.00000000000000 2.50000000000000 5.00000000000000 0.500000000000003.000000000000004.37500000000000 0.89743589743590 4.875000000000004.00000000000000 11.75781250000000 1.70828603859251 -6.882812500000005.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相差越来越大，即迭代序列}{k x 发散，此迭代法就失败.这时偏差piancha 逐渐增大且偏差的相对误差xdpiancha 的值大于0.5.用迭代公式)2/(101+=+k k x x 运行后输出的结果 [k,piancha,xdpiancha,xk]=1.00000000000000 0.50000000000000 0.200000000000002.500000000000002.00000000000000 0.27777777777778 0.12500000000000 2.222222222222223.00000000000000 0.14619883040936 0.06172839506173 2.368421052631584.00000000000000 0.07926442612555 0.03462603878116 2.289156626506025.00000000000000 0.04230404765128 0.01814486863115 2.33146067415730k =5,piancha =0.04230404765128,xdpiancha = 0.01814486863115, xk = 2.33146067415730可见，偏差piancha 和偏差的相对误差xdpiancha 的值逐渐变小，且第5次的迭代值x k =2.331 460 674 157 30与根*x =2.316 624 790 355 40接近，则迭代序列}{k x 收敛，但收敛速度较慢，此迭代法较为成功.用迭代公式)22/()102(21+-+-=+k k k k k x x x x x 运行后输出的结果[k,piancha,xdpiancha,xk ]=1.00000000000000 0.33333333333333 0.142857142857142.333333333333332.00000000000000 0.01666666666667 0.007194244604322.316666666666673.00000000000000 0.00004187604690 0.000018076318222.316624790619774.00000000000000 0.00000000026437 0.000000000114122.316624790355405.00000000000000 0 0 2.31662479035540 祝贺您！此迭代序列收敛，且收敛速度较快k = 5; piancha = 0; xdpiancha = 0; y = 2.31662479035540.可见，偏差piancha 和偏差的相对误差xdpiancha 的值越来越小，且第5次的迭代值x k =2.316 624 790 355 40与根*x =2.316 624 790 355 40相差无几，则迭代序列}{k x 收敛，且收敛速度很快，此迭代法成功.7.4.2迭代法的MATLAB 程序2迭代法的MATLAB 主程序2function [k,piancha,xdpiancha,xk,yk]=diedai2(x0,tol,ddmax) x(1)=x0;for i=1: ddmaxx(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<tol)|(xdpiancha< tol) k=i-1; xk=x(i); return ; end endif i>ddmaxdisp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x(i);yk=fun(x(i));[(i-1) piancha xdpiancha xk yk]; return ; endP=[(i-1),piancha,xdpiancha,xk,yk]';例7.4.2 求x 5-3x +1=0在0.3附近的根，精确到4位小数. 解 ⑴ 构造迭代公式3/)1(51+=+k k x x ),2,1,0( =k .⑵利用迭代法的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.33473k =3;piancha =1.206089525390697e-005; xdpiancha =3.603129477781680e-005;xk =0.3347; yk =0.3347.7.5 迭代过程的加速方法及其MATLAB 程序7.5.1加权迭代法的MATLAB 程序加权迭代法的MATLAB 主程序现提供名为jasudd.m 的M 文件如下：function [k,xk,yk]=jasudd (x0,tol,L,ddmax) x1(1)=x0;for i=1: ddmaxx(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<tol)|(xdpiancha< tol) k=i-1; xk=x1(i); return ;end endif i>ddmaxdisp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x1(i); [(i-1) xk yk]; return ; endP=[(i-1),xk,yk]';例7.5.1 求2e 0=--x x在0.85附近的一个近似根，要求精度610-=ε. 解 在MATLAB 工作窗口输入程序>> [k,xk,yk]=jasudd (0.85,1e-6,-0.855,100)运行后输出结果[k,xk,yk] =1.00000000000000 0.85260370028041 0.852607038305612.00000000000000 0.85260549975491 0.852605504062363.00000000000000 0.85260550207699 0.85260550208255 k =3； xk =0.852606； yk = 0.852606.7.5.2 艾特肯(Aitken)加速方法的MATLAB 程序艾特肯(Aitken)加速方法的MATLAB 主程序 现提供名为Aitken.m 的M 文件如下：function [k,xk,yk,p]= Aitken (x0,tol, ddmax) x(1)=x0;for i=1: ddmaxx1(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<tol)|(xdpiancha<tol) k=i-1; xk=x(i); yk=fun(x(i));m=[0,1:i-1]; p=[m',x1',x2',x'];return ;endendif i>ddmaxdisp('迭代次数超过给定的最大值ddmax') k=i-1; xk=x(i); yk=fun(x(i)); m=[0,1:i-1]; p=[m',x1',x2',x']; return ; endm=[0,1:i-1]; p=[m',x1',x2',x'];例7.5.3 用艾特肯加速方法求2e0=--x x在0.85附近的一个近似根，要求精度610-=ε.解 在MATLAB 工作窗口输入程序>> [k,xk,yk,p]= Aitken(0.85,1e-6, 100)运行后输出结果k=3，xk=0.85260550201343，yk=0.85260550201373p = 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.852605502013737.6 牛顿(Newton)切线法及其MATLAB 程序7.6.1 牛顿切线法的收敛性及其MATLAB 程序牛顿切线法的收敛性及其MATLAB 主程序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 如下')elsedisp('请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|和方程f(x)=0的函数f(x)值')endP=[y,f]';例 7.6.1 用牛顿切线法的局部收敛性判别方程 e 4sin =x x 的近似根时，由下列初始值0x 产生的迭代序列是否收敛？⑴-1;0=x ⑵0;0=x ⑶1;0=x ⑷2;0=x ⑸ 5.50=x ；⑹80=x . 解 在MATLAB 工作窗口输入程序>> [y,f]=newjushou(-1)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值y =139.5644 f =4.3096（2）输入程序>> [y,f]=newjushou(0)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)的值fy =8.0000 f = 4（3）输入程序>> [y,f]=newjushou(1)运行后输出结果恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值f 如下y =0.3566 f =1.7126（4）输入程序：>> [y,f]=newjushou(2)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值y =1.2593 f =-2.7188（5）输入程序>> [y,f]=newjushou(5.5)运行后输出结果请注意观察下面显示的φ(x)的导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值y =1.0447e+005 f =176.6400（6）输入程序>> [y,f]=newjushou(8)运行后输出结果恭喜您！此迭代序列收敛,φ(x)导数值的绝对值y=|d φ(x)/dx|和方程f(x)=0的函数f(x)值f 如下y =0.4038 f =-2.9452e+0037.6.2 牛顿切线法的MATLAB 程序牛顿切线法的MATLAB 主程序现提供名为newtonqx.m 的M 文件： function[k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax)x(1)=x0;for i=1: gxmaxx(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)&((piancha<tol)|(xdpiancha< tol)) k=i-1; xk=x(i);[(i-1) xk yk piancha xdpiancha] return ; end endif i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax 。