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

方程（组）的根及其命求解方程（组）的solve命x=solve('方程f(x)=q(x)','待求符号变量>>E1=sym('方程En=sym('方程[x1,x2,…,xn]=solve(E1,E2,…,En,求解方程（组）的fsolve命fsolve的调用格式搜索根的方法及其程作图法及其程作函y

f(x)在区间[a,b]的图形的程序x=a:h:b;h是步长y=f(x);grid,说明：⑴此程序在的工作区输入，运行后即可出现函数yf(x)的图形.此图x轴交点的横坐标即为所要求的根的近似值.a,b]b-ah作函y

f(x)在区间[a,b]上的图形的程序yf(xh(x)g(x)，其中h(x)g(xx=a:h:b;y1=h(x);y2=g(x);plot(x,y1,x,y2)grid,gtext('y1=h(x),y2=g(x)')说明：此程序 的工作区输入，运行后即可出现函数y1h(x)和y2逐步搜索法及其程逐步搜索法的主程function输入的量a和b是闭区间[a,b]%---h是步长%---tol是预先给定的精度运行后输出的量---k%---r[a,b]上的实根的近似值，其精度是X=a:h:b;Y=funs(X);n=(b-forX(k)=a+k*h;Y(k)=funs(X(k));程序中调用的funs.m为函数ifsk<=0,xielv=(Y(k+1)-Y(k))*(Y(k)-Y(k-if(abs(Y(k))<tol)&(xielv<=0)例7.2.1用逐步搜索法的程序分别求方程2x32x23x30

2x30在区间[2,2上的根的近似值，要求精度是0.000解将逐步搜索法的程序保存名为zhubuss.m的M文件.functiony=2.*x.^3+2.*x.^2-3.*x->>[k,r]=zhubuss(-kr=-1.2240-1.0000-1.0000- 即搜索点的个数为k=4001，其中有5个是方程2x32x23x30的近似根，即r-1.2240，-1.0000，-1.0000，-0.9990，1.2250，其精度为0.000在程序中将y=2.*x.^3+2.*x.^2-3.*x-3用y=sin(cos(2.*x.^3))

2x3)0在区间[2,2上的根的近似值如r=-1.9190-1.7640-1.5770-1.3300- 1.3310 二分法及其程二分法的程二分法的主程function[k,x,wuca,yx]=erfen(a,b,abtol)a(1)=a;b(1)=b;ya=fun(a(1));yb=fun(b(1));程序中调用的fun.mifya*disp('注意：ya*yb>0,请重新调整区间端点a和b.'),returnmax1=-1+ceil((log(b-a)log(abtol))/log(2));ceilfork=1:a;ya=fun(a);b;yb=fun(b);x=(a+b)/2;yx=fun(x);wuca=abs(b-a)/2;k=k-1;ifa=x;b=x;elseifa=x;ifb-a<abtol,return,endk=max1;x;wuca;7.3.1x3-x+4=02,-1)内的0.001.（1）先用两种方法确定方x3-x+4=0的实根的分布情况。1作图法.>>x=-y=x.^3-x+4;plot(x,y)grid,gtext('y=x^3-f(x)=x3-x+4的图像.从图像可以看出，此曲线有两个驻点

3x3在(-2,-1x2,-12试算>>x=-4:1:4,y=x.^3-xx----01234y---444f(xf(2f(1)0，所以此方程在(-2,-1)（2）用二分法的主程>>[k,x,wuca,yx]=erfen(-2,-2-3，k=9，x=-wuca=9.7656e-004，yx=表2-bk20-2.000-1.000-1.5000.500-2.0004.0002.1251-2.000-1.500-1.7500.250-2.0002.1250.3902-2.000-1.750-1.8750.125-2.0000.390-0.7163-1.875-1.750-1.8120.062-0.7160.390-0.1414-1.812-1.750-1.7810.031-0.1410.3900.1295-1.812-1.781-1.7960.015-0.1410.129-0.0046-1.796-1.781-1.7890.007-0.0040.1290.0627-1.796-1.789-1.7930.003-0.0040.0620.0298-1.796-1.793-1.7940.002-0.0040.0290.0129-1.796-1.794-1.7950.001-0.0040.0120.003迭代法及其程迭代法的程序迭代法的主程序function输入的量--x0是初始值,kforpiancha=abs(x(i+1)-x(i));xdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);[(i-1)pianchaxdpianchaxk]if(pianchaif(piancha<0.001)&(xdpiancha<0. p=[(i-1)pianchaxdpiancha7.4.1f(x)x22x10解在工作窗口输入程k>>[k,piancha,xdpiancha,xk]=k

k=5，piancha=11.80374145507813，xdpiancha=0.63167031671297，xk=-18.68655395507813由以上运行后输出的迭代序列与根x*=2.31662479035540相差越来越大，即迭代序列{xk}发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的0.5.

xk110xk2k=5,piancha=0.04230404765128,xdpiancha=0.01814486863115,xk=2.33146067415730可见，偏差piancha和偏差的相对误差xdpiancha的值逐渐变小，且第5的迭代值k=2.33146067415730x*=2.31662479035540{x收敛，但收k

xk1xk(xk22xk10)

0k=5;piancha=0;xdpiancha=0;y=kpianchaxdpiancha5次的迭代值xk=2.31662479035540x*=2.31662479035540相差无几，则迭代序列{x收敛，且k迭代法的程序迭代法的主程序function[k,piancha,xdpiancha,xk,yk]=diedai2(x0,tol,ddmax)fori=1:ddmaxxdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fun(x(i));[(i-1)pianchaxdpianchaxkyk]if(piancha<tol)|(xdpiancha<tol)k=i-1;xk=x(i);ifk=i-1;xk=x(i);yk=fun(x(i));[(i-1)pianchaxdpianchaxkyk];P=[(i-例7.4.2求x5-3x+1=0在0.3附近的根，精确到4位小数造kxk1(x51)/k

(k0,1,2,)⑵利用迭代法 程序2计算精确到4位小数的根的近似值y和迭代次数>>[k,piancha,xdpiancha,xk,yk]=diedai2(0.3,1e-[k,piancha,xdpiancha,xk,yk]00.034140.102180.3000010.034140.102180.3341420.000580.001730.3347230.000010.000040.334730.33473k=3;piancha=1.206089525390697e-005;xdpiancha=3.603129477781680e-005;xk=0.3347;yk=0.3347.迭代过程的加速方法及其程7.5.1迭代法的程迭代法的主程jasudd.mMfunction[k,xk,yk]=jasudd(x0,tol,L,ddmax)fori=1:x1(i+1)=x(i+1)+L*(x(i+1)-x1(i))/(1-piancha=abs(x1(i+1)-xdpiancha=piancha/(abs(x1(i+1))+eps);i=i+1;xk=x1(i);yk=fun(x1(i));[(i-1)xkyk]if(piancha<tol)|(xdpiancha<tol)k=i-1;xk=x1(i);ifk=i-1;xk=x1(i);[(i-1)xkyk];P=[(i-7.5.12exx00.85附近的一个近似根，要求精度106解在工作窗口输入程>>[k,xk,yk]=jasudd(0.85,1e-6,-[k,xk,yk]1.000000000000000.852603700280412.000000000000000.852605499754913.000000000000000.852605502076990.85260550208255k=3；xk=0.852606；yk=0.852606.7.5.2(Aitken)加速方法的程(Aitken)加速方法的主程Aitken.m的Mfunction[k,xk,yk,p]=Aitken(x0,tol,ddmax)fori=1:x1(i+1)=fun(x(i));x(i+1)=x2(i+1)-(x2(i+1)-x1(i+1))^2/(x2(i+1)-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'];ifk=i-1;xk=x(i);yk=fun(x(i));m=[0,1:i-1];p=[m',x1',x2',x'];m=[0,1:i-1];例 用加速方法求2exx0在0.85附近的一个近似根，要求精106解在工作窗口输入程>>[k,xk,yk,p]=Aitken(0.85,1e-6,p 7.6(Newton)切及其程切的收敛性及其程 的收敛性及其主程functionf=fnq(x);fz=fnq(x)*ddfnq(x)/((dfnq(x))^2+eps);if例7.6.1用切的局部收敛性判别方程ex始值x0产生的迭代序列是否收敛

x4⑴x0-

⑵x0

⑶x0

⑷x02;

x05.5；⑹x08 >>[y,f]=newjushou(-请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|f(x)=0的函数f(x)y=f>>请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|f(x)=0的函数f(x)的值yf=>>！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|f(x)=0的函数f(x)fyf=>>请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|f(x)=0的函数f(x)yf->>请注意观察下面显示的φ(x)的导数值的绝对值y=|dφ(x)/dx|f(x)=0的函数f(x)y=f=>>！此迭代序列收敛,φ(x)导数值的绝对值y=|dφ(x)/dx|f(x)=0的函数f(x)fyf- 的程 主程newtonqx.mMfori=1:xdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fnq(x(i));[(i-1)xkykpianchaxdpiancha]if(abs(yk)<ftol)&((piancha<tol)|(xdpiancha<tol))k=i-1;xk=x(i);[(i-1)xkykpianchaxdpiancha]ifdisp('请注意：迭代次数超过给定的最大值gxmax。 xk=x(i);[(i-1)xkykpianchaxdpiancha][(i-例7.6.2用 求方程2x33x210在x00.4和0.9的近似根，要求精度103.解>>[k,xk,yk,piancha,xdpiancha]=newtonqx(-0.4,0.001,x00.4求nc的方法及其程求nc的方法及其主程kai2fang.m的Mfunction[k,xk,yk,piancha,xdpiancha,P]=kainfang(x0,c,n,tol,fori=1:u(i)=(x(i)^n-c)/(n*x(i)^(n-1));x(i+1)=x(i)-u(i);xdpiancha=piancha/(abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fnq(x(i));[(i-if(piancha<tol)|(xdpiancha<tol)k=i-1;xk=x(i);yk=fnq(x(i));ifk=i-1;xk=x(i);yk=fnq(x(i));P=[(i-7.6.3求113，要求精度为106方法1用求c的n次根nc（当n是偶数时，c0）的程序计算x010n2，被开方数c113，tol=105，迭代的gxmax=100.在工作区间输入程序>>[k,xk,yk,piancha,xdpiancha]=kainfang(10,113,2,1e-k=4,piancha=1.610800381968147e-xdpiancha=1.515313534117706e-xk=10.63014581273465,yk可见 10.63015，满足精度105方法 迭代计算x 113x

x21130

f(x)x2113，f'(x)2x.由迭代得x2xk1xk xk1

（xk121

）(其中kx010偏差xk1偏差xk110.65010.65020.01910.63030.00010.63040.00010.630108，满足精度105，所以 10.630方法3用切的主程序fnq.mdfnq.mM件functionfunctiony=dfnq(x)在工作窗口输入程>>[k,xk,yk,piancha,xdpiancha]=newtonqx(10,1e-5,1e-运行后，将输出的结果列入下表2-13.迭代k=4次，得到精度为105的结 10.630k10.6500.06110.6500.42220.0190.00110.6300.00030.0000.00010.6300.00040.0000.00010.6300.000方法4在工作空间输入程>>ans=经过四舍五入后，得到精度为105的结 10.630割及其程 的程割 主程gexian.m的Mfunction[k,piancha,xdpiancha,xk,yk]=gexianfori=2:u(i)=fnq(x(i))*(x(i)-x(i-1));v(i)=fnq(x(i))-fnq(x(i-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)pianchaxdpianchaxkyk]if(abs(yk)<ftol)&((piancha<tol)|(xdpiancha<tol))k=i-2;xk=x(i);yk=fnq(x(i));[(i-2)pianchaxdpianchaxkyk]；ifk=i-2;xk=x(i);yk=fnq(x(i));

抛 及其程抛 的程抛物 主程function[k,piancha,xdpiancha,xk,yk]=paowu(x1,x2,X=[x1,x2,fori=1:h0=X(1)-X(3);h1=X(2)-Xu0=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;b