实验指导书非线性方程(组)的数值解法

1、仅供个人参考实验指导书非线性方程(组)的数值解法隔根区间的求法(确定根的初始近似值)：作图法，逐步搜索法等求根的方法：二分法，迭代法，牛顿法，割线法，抛物线法，迭代法的加速等Forpersonaluseonlyinstudyandresearch;notforcommercialuse实验平台：MATLAB件说明：与方程求根有关的MATLAB函数另见MATLAB学习资料Forpersonaluseonlyinstudyandresearch;notforcommercialuse一、搜索根的方法及其MATLAB序求解非线性方程根的近似值时，首先需要判断方程有没有根？如果有根，有几个根？如果有根

2、，需要搜索根所在的区间或确定根的初始近似值(简称初始值).搜索根的近似位置白常用方法有三种：作图法、逐步搜索法和二分法等，使用这些方法的前提是高等数学中的零点定理.1 .作图法及其MATLAB序(1)作函数y=f(x)在区间a,b的图形的MATLAB程序一x=a:h:b;%h是步长y=f(x);plot(x,y)grid,gtext('y=f(x)')说明：此程序在MATLAB的工作区输入，运行后即可出现函数y=f(x)的图形.此图形与x轴交点的横坐标即为所要求的根的近似值.区间a,b的两个端点的距离b-a和步长h的绝对值越小，图形越精确.(2)作函数y=f(x)在区间a,b的

3、图形的MATLAB程序二将y=f(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的工作区输入，运行后即可出现函数y1=h(x)和y2=g(x)的图形.两图形交点的横坐标即为所要求的根的近似值.2 .逐步搜索法及其MATLAB序逐步搜索法也称试算法.它是求方程f(x)=0根的近似值位置的一种常用的方法.逐步搜索法依赖于寻找连续函数f(x)满足f(a)与f(b)异号的区间a,b.一旦找到区间，无论区

4、间多大，通过某种方法总会找到一个根.MATLAB的库函数中没有逐步搜索法的程序，现根据逐步搜索法的计算步骤和它的收敛判定准则编写其MATLAB程序如下，命名为zhubussm.输入输出说明：输入区间端点a和b的值，步长h和精度tol,运行后输出迭代次数k=(b-a)/h+1,方程f(x)=0根的近似值r.functionk,r=zhubuss(a,b,h,tol)%输入的量：%琲Db是闭区间a,b的左、右端点；%h是步长；%tol是预先给定的精度.%运行后输出的量：%k是搜索点的个数；%r是方程在a,b上的实根的近似值，其精度是tol;X=a:h:b;Y=funs(X);%f(x)的M文件，对

5、应名称为funs.m的函数文件n=fix(b-a)/h)+1;m=0;X(n+1)=X(n);Y(n+1)=Y(n);fork=2:nX(k)=a+k*h;Y(k)=funs(X(k);%程序中调用的funs.m为函数sk=Y(k)*Y(k-1);ifsk<=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);endend例1用逐步搜索法的MATLAB程序分别求方程2x3+2x2-3x-3=0和sin(cos(2x3)=0在区间-

6、2,2上的根的近似值，要求精度是0.0001.解将逐步搜索法的MATLAB程序保存名为zhubuss.m的M文件.建立M文件funs.mfunctiony=funs(x)y=2.*xA3+2.*xA2-3.*x-3在MATLAB工作窗口输入如下程序>>k,r=zhubuss(-2,2,0.001,0.0001)运行后输出的结果k=4001r=-1.2240-1.0000-1.0000-0.99901.2250即搜索点的个数为k=4001,其中有5个是方程2-+2丁-31-3=0的近似根，即r=-1.2240,-1.0000,-1.0000,-0.9990,1.2250,其精度为0.

7、0001.在程序中将y=2.*x.A3+2.*x.A2-3.*x-3用y=sin(cos(2.*x.A3)代替，可得到方程sin(cos(2x3)=0在区间卜2,2上的根的近似值如下r=-1.9190-1.7640-1.5770-1.3300-0.92200.92301.33101.57801.76501.9200二分法及其MATLAB程序1.二分法的MATLAB程序二分法的MATLAB主程序一：functionk,x,wuca,yx=erfen(a,b,abtol)a(1)=a;b(1)=b;ya=fun(a(1);yb=fun(b(1);%f(x)的M文件，对应名称为fun.m的函数文件i

8、fya*yb>0,disp('注意：ya*yb>0,请重新调整区间端点a和b.'),returnendmax1=-1+ceil(log(b-a)-log(abtol)/log(2);%ceil是向+00方向取整%max1符合误差要求的迭代次数fork=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,yxifyx=0a=x;b=x;elseifyb*yx>0b=x;yb=yx;elsea=x;ya=yx;endifb-a<

9、;abtol,return,endendk=max1;x;wuca;yx=fun(x);(-2,-1)内例3确定方程x3-x+4=0的实根的分布情况，并用二分法求在开区问的实根的近似值，要求精度为0.001.解建立M文件fun.mfunctiony=fun(x)y=x.A3-x+4(1)先用两种方法确定方程x3-x+4=0的实根的分布情况。方法1作图法.在MATLAB工作窗口输入如下程序>>x=-4:0.1:4;y=x.A3-x+4;plot(x,y)grid,gtext('y=xA3-x+4')画出函数f(x)=x3-x+4的图像.从图像可以看出，此曲线有两个驻点

10、士1都在x轴的上方，在(-2,-1)内曲线与x轴只有一个交点，则该方程有唯一一个实根，且在(-2,-1)内.方法2试算法.在MATLAB工作窗口输入程序>>x=-4:1:4,y=x.A3-x+4运行后输出结果x=-4-3-2-101234y=-56-20-2444102864由于连续函数f(x)满足,所以此方程在(-2,-1)内有一个实根.(2)用二分法的主程序计算.在MATLAB工作窗口输入程序>>k,x,wuca,yx=erfen(-2,-1,0.001)运行后屏幕显示用二分法计算的过程(现已被列入下表)，其余结果为k=9,x=-1.7959,wuca=9.7656

11、e-004yx=0.0037次数k左端点ak右端点bk中点xk函数值f(ak)函数值f(bk)函数值f(xk)0-2.0000-1.0000-1.50000.5000-2.00004.00002.12501-2.0000-1.5000-1.75000.2500-2.00002.12500.39062-2.0000-1.7500-1.87500.1250-2.00000.3906-0.71683-1.8750-1.7500-1.81250.0625-0.71680.3906-0.14184-1.8125-1.7500-1.78130.0313-0.14180.39060.12965-1.8125

12、-1.7813-1.79690.0156-0.14180.1296-0.00486-1.7969-1.7813-1.78910.0078-0.00480.12960.06277-1.7969-1.7891-1.79300.0039-0.00480.06270.02908-1.7969-1.7930-1.79490.0020-0.00480.02900.01219-1.7969-1.7949-1.79590.0010-0.00480.01210.0037二分法的MATLAB主程序二：functionc,err,yc,k=bisect(f,a,b,delta)ya=feval(f,a);yb=fe

13、val(f,b);ifya*yb>0,return,endmax1=1+round(log(b-a)-log(delta)/log(2);fork=1:max1c=(a+b)/2;yc=feval(f,c);ifyc=0a=c;b=c;elseifyb*yc>0b=c;yb=yc;elsea=c;ya=yc;endifb-a<delta,break,endendc=(a+b)/2;err=abs(b-a);yc=feval(f,c);三、迭代法及其MATLAB程序1.迭代法的MATLAB程序1functionk,piancha,xdpiancha,xk=diedai1(x0,

14、k)%输入的量-x0是初始值，k是迭代次数x(1)=x0;fori=1:kx(i+1)=fun1(x(i);%程序中调用的funl.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)pianchaxdpianchaxkendif(piancha>1)&(xdpiancha>0.5)&(k>3)disp(请用户注意：此迭代序列发散，请重新输入新的迭代公式')return;endif(piancha<0.001)&(x

15、dpiancha<0.0000005)&(k>3)disp(祝贺您！此迭代序列收敛，且收敛速度较松)return;endp=(i-1)pianchaxdpianchaxk'例4求方程，0)二#+2尸10的一个正根.解首先建立迭代函数的M文件(注意调用主程序时要相应修改函数名)functiony=fun1(x)y=(10-xA2)/2;functiony=fun2(x)y=10/(x+2);在MATLAB工作窗口输入程序>>k,piancha,xdpiancha,xk=diedai1(2,5)运行后输出用迭代公式仆+1=(1°一只”2的结果k,p

16、iancha,xdpiancha,xk=1.000000000000001.000000000000000.3333333.000000000000002.000000000000002.0000005.000000000000000.0000003.000000000000004.0000000.7435904.0000004.0000000000000011.0000001.859251-6.0000005.0000000000000011.5078130.671297-18.507813请用户注意：此迭代序列发散，请重新输入新的迭代公式k=5,piancha=11.507813xdpia

17、ncha=0.67129%xk=-18.507813由以上运行后输出的迭代序列与根.=2.31662479035540相差越来越大，即迭代序列发散，此迭代法就失败.这时偏差piancha逐渐增大且偏差的相对误差xdpiancha的值大于0.5.不得用于商业用途用迭代公式工同=11M敌+2)运行后输出的结果k,piancha,xdpiancha,xk=1.000000000000000.0000000.0000002.0000002.000000000000000.7777780.0000002.2222223.000000000000000.1460.0061732.2631584.00000

18、0000000000.0125550.0781162.2895.000000000000000.0651280.0182.415730k=5,piancha=0.065128,xdpiancha=0.018,xk=2.415730可见，偏差piancha和偏差的相对误差xdpiancha的值逐渐变小，且第5次的迭代值xk=2.33146067415730与根，=2.3166247903554哦近，则迭代序列演)收敛，但收敛速度较慢，此迭代法较为成功.用迭代公式海+1二小一(aJ+2/-a+2)运行后输出的结果k,piancha,xdpiancha,xk=1.000000000000000.33

19、33330.2857142.3333332.000000000000000.0666673.000000000000000.0000900.0072.6666670.0000222.0619774.000000000000000.000000000264370.000000000114122.0355405.0000000000000002.035540祝贺您！此迭代序列收敛，且收敛速度较快k=5;piancha=0;xdpiancha=0;y=2.035540.可见，偏差piancha和偏差的相对误差xdpiancha的值越来越小，且第5次的迭*代值xk=2.31662479035540与根

20、工=2.31662479035540相差无几,则迭代序列%)收敛，且收敛速度很快，此迭代法成功.2.迭代法的MATLAB程序2迭代法的MATLAB主程序2functionk,piancha,xdpiancha,xk,yk=diedai2(x0,tol,ddmax)x(1)=x0;fori=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)pianchaxdpianchaxkykif(piancha<tol)|(xd

21、piancha<tol)k=i-1;xk=x(i);return;endendifi>ddmaxdisp(迭代次数超过给定的最大值ddmax)k=i-1;xk=x(i);yk=fun(x(i);(i-1)pianchaxdpianchaxkyk;return;endP=(i-1),piancha,xdpiancha,xk,yk'例5求x5-3x+1=0在0.3附近的根，精确到4位小数.解构造迭代公式利用迭代法的MATLAB程序2计算精确到4位小数的根的近似值y和迭代次数k.在MATLAB工作窗口输入程序>>k,piancha,xdpiancha,xk,yk=di

22、edai2(0.3,1e-4,100)运行后输出的结果k,piancha,xdpiancha,xk,yk=0.30000 0.334140.33414 0.334720.33472 0.334730.33473 0.3347300.034140.1021810.034140.1021820.000580.0017330.000010.00004k=3;piancha=1.5390697e-005;xdpiancha=3.7781680e-005;xk=0.3347;yk=0.3347.四、迭代过程的加速方法及其MATLAB程序1.艾特肯(Aitken)斯蒂芬森加速方法的MATLAB程序Aitk

23、en-Steffensen加速方法的MATLAB主程序现提供名为Aitkstf.m的M文件如下：functionk,xk,yk,p=Aitkstf(x0,tol,ddmax)x(1)=x0;fori=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)A2/(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<

24、;tol)|(xdpiancha<tol)k=i-1;xk=x(i);yk=fun(x(i);m=0,1:i-1;p=m',x1',x2',x'return;endendifi>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用艾特肯加速方法求23-1二。在0.85附近的一个近似根，要求精度二：.解首先建立

25、迭代函数的M文件(注意函数名要与主程序中的函数名一致，或者在调用主程序时要相应修改函数名),然后在MATLAB工作窗口输入程序>>k,xk,yk,p=Aitkstf(0.85,1e-6,100)运行后输出结果p =01.000000000000002.000000000000003.00000000000000k=3,xk=0.201343,yk=0.2013730.000000000.3897450.6524840.5686070.4918110.1508260.2014070.2013430.2013980.201373五、牛顿(Newton)切线法及其MATLAB程序1 .牛

26、顿切线法的收敛性及其MATLAB程序functiony,f=newshou(x)%fnq.m为函数f(x)的M文件，dfnq.m为函数f(x)的一阶导数的M文件%ddfnq.m为函数f(x)的二阶导数的M文件f=fnq(x);fz=fnq(x)*ddfnq(x)/(dfnq(x)A2+eps);y=abs(fz);if(y<1)disp('恭喜如止匕迭代序列收敛，日期数值的绝对值y=|d小(x)/dx|!方程f(x)=0的函数f(x)值f如下')elsedisp('请注意观察下面显示的小(x的导数值的绝对值y=|d小(x)/dx1方程f(x)=0的函数f(x)值)

27、endP=y,f'例8用牛顿切线法的局部收敛性判别方程exsinx=4的近似根时，由下列初始值人产生的迭代序列是否收敛？(1)1口二/；/二Q；1口=L(4)题=2,(5)风=5.5；%=8.解(1)在MATLAB工作窗口输入程序>>y,f=newshou(-1)运行后输出结果请注意观察下面显示的小(x)勺导数值的绝对值y=|d小(x)/dx1口方程f(x)=0的函数f(x)值y=139.5644f=4.3096(2)输入程序>>y,f=newshou(0)运行后输出结果请注意观察下面显示的小(x)勺导数值的绝对值y=|d小(x)/dx1口方程f(x)=0的函数

28、f(x)的值fy=8.0000f=4(3)输入程序>>y,f=newshou(1)运行后输出结果恭喜您！此迭代序列收敛，小(x异数值的绝对值y=|d小(x)/d洲方程f(x)=0的函数f(x)值f如下y=0.3566f=仅供个人参考1.7126( 4) 输入程序：>>y,f=newshou(2)运行后输出结果请注意观察下面显示的小(x)勺导数值的绝对值y=d小(x)/dX|口方程f(x)=0的函数f(x)值y=1.2593f=-2.7188( 5) 输入程序>>y,f=newshou(5.5)运行后输出结果请注意观察下面显示的小(x)勺导数值的绝对值y=|d

29、小(x)/dx1口方程f(x)=0的函数f(x)值y=1.0447e+005f=176.6400( 6) 输入程序>>y,f=newshou(8)运行后输出结果恭喜您！此迭代序列收敛，小(x异数值的绝对值y=|d6(x)/d洲方程f(x)=0的函数f(x)值f如下y=0.4038f=-2.9452e+0032牛顿切线法的MATLAB程序牛顿切线法的MATLAB主程序一:现提供名为newtonqx.m的M文件：functionk,xk,yk,piancha,xdpiancha=newtonqx(x0,tol,ftol,gxmax)x(1)=x0;fori=1:gxmaxx(i+1)=

30、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)xkykpianchaxdpianchaif(abs(yk)<ftol)&(piancha<tol)|(xdpiancha<tol)k=i-1;xk=x(i);(i-1)xkykpianchaxdpianchareturn;endendifi>gxmaxdisp(请注意：迭代次数超过给定的最大值gxmax。)k=i-1;xk=

31、x(i);(i-1)xkykpianchaxdpianchareturn;end(i-1),xk,yk,piancha,xdpiancha'例9用牛顿切线法求方程2-3/+1=0在的二-0,4和0.9的近似根，要求精度二二1解在MATLAB工作窗口输入程序>>k,xk,yk,piancha,xdpiancha=newtonqx(-0.4,0.001,0.001,100)运行后输出初始值而二14的迭代结果.牛顿切线法的MATLAB主程序二:functionp0,err,k,y=newton(f,df,p0,delta,epsilon,max1)fork=1:max1p1=p0

32、-feval(f,p0)/feval(df,p0);err=abs(p1-p0);relerr=2*err/(abs(p1)+delta);p0=p1;y=feval(f,p0);if(err<delta)|(relerr<delta)|(abs(y)<epsilon),break,endend3.求加的方法及其MATLAB程序现提供名为kai2fang.m的M文件：functionk,xk,yk,piancha,xdpiancha,P=kainfang(x0,c,n,tol,gxmax)x(1)=x0;fori=1:gxmaxu(i)=(x(i)An-c)/(n*x(i)A

33、(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,xdpianchaif(piancha<tol)|(xdpiancha<tol)k=i-1;xk=x(i);yk=fnq(x(i);return;endendifi>gxmaxdisp(请注意：迭代次数超过给定的最大值gxmax.')k=i-1;xk=x(i);yk=fnq(x(i);(i-1),xk,yk,pian

34、cha,xdpianchareturn;endP=(i-1),xk,yk,piancha,xdpiancha'例10求叵,要求精度为10*.解本题介绍四种解法.方法1用求C的册次根唬（当n是偶数时，e>0）的MATLAB程序计算.取初始值即二10,根指数胃二2,被开方数1r=113,近似根的精度toi=10-',迭代的最大次数gxmax=100.在工作区间输入程序>>k,xk,yk,piancha,xdpiancha=kainfang（10,113,2,1e-5,100）运行后输出结果k=4,piancha=1.1968147e-011,xdpiancha=1

35、.5153e-012xk=10.630,yk=1.7109e+009可见，10.63015,满足精度11一'.方法2用牛顿迭代公式(2.12)计算.设Tx=H3,则/-113=0,记拉）=7-113,/任）*.由牛顿迭代公X-113V113、工上川二办-工=-I工止+-.式得，端，即2"“（其板二0,123）取初始值刖二10,计算结果列入下表迭代次数k偏差根的近似值现10.65000010.65000020.01983610.63016430.000019110.63014640.00000010.630146因为，迭代次数七=4时，偏差卜4-讨<10,满足精度104,

36、所以，行主10.63015.方法3用牛顿切线法的MATLAB主程序计算.分别建立名为fnq.m和dfnq.m的M文件functiony=fnq(x)y=xA2-113;functiony=dfnq(x)y=2*x;在MATLAB工作窗口输入程序>>k,xk,yk,piancha,xdpiancha=newtonqx(10,1e-5,1e-5,100)运行后，将输出的结果列入下表.迭代k=4次，得到精度为10,的结果10.63015.pianchaxdpianchaxkyk0.6500000.06103310.6500000.4225000.0198360.00186610.6301

37、640.0003930.0000190.00000210.6301460.0000000.0000000.00000010.6301460.000000方法4在MATLAB工作空间输入程序>>113A(1/2)运行后输出ans=10.630经过四舍五入后，得到精度为的结果711310.63015.4.牛顿切线法的加速及其两种MATLAB程序(选学)(一)已知方程根的重数加已知方程f(x)=0根的重数血，求重根的修正牛顿切线法的MATLAB主程序现提供名为newtonxz.m的M文件：functionk,piancha,xdpiancha,xk,yk=newtonxz(m,x0,to

38、l,ftol,gxmax)x(1)=x0;fori=1:gxmaxx(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);(i-1)pianchaxdpianchaxkyk;if(piancha<tol)|(xdpiancha<tol)&(abs(yk)<ftol)k=i-1;xk=x(i);yk=fnq(x(i);(i-1)pianchaxdpianchaxkyk;return

39、;endendifi>gxmaxdisp(请注意：迭代次数超过给定的最大值gxmax.')k=i-1;xk=x(i);yk=fnq(x(i);(i-1)pianchaxdpianchaxkyk;return;end例11判断二0是方程，(x)=2eA-9x2-6工-2=0的几重根？在区间0上，分别用牛顿切线法和求重根的修正牛顿切线公式求此根的近似值牛，使精确到，I解经判断知，二0是方程1)二0的三重根.在MATLAB工作窗口输入程>>k,piancha,xdpiancha,xk,yk=newtonqx(0.5,1e-4,1e-4,100)运行后整理结果列入下表kpia

40、nchaxdpianchaxkyk10.144100.404890.355900.5419820.107410.432250.248490.1682030.077450.452800.171040.0514640.054500.467600.116550.0155850.037690.477980.078850.0046960.025760.485140.053100.0014070.017460.490010.035630.0004280.011770.493300.023860.0001290.007910.495520.015960.00004010.005300.497000.0106

41、60.00001110.003540.498000.007120.00000210.002370.498670.004750.00000310.001580.499110.003170.00000410.001050.499410.002110.00000510.000700.499600.001410.00000610.000470.499740.000940.00000710.000310.499820.000630.00000810.000210.499880.000420.00000910.000140.499920.000280.00000020.000090.499950.0001

42、90.00000迭代次数k=20,精确到日=1。-4的根的近似值是xk=0.0002,其函数值是yk=5.7520e-011,收敛速度较慢.根据求重根的修正牛顿切线公式，在MATLAB工作窗口输入程序>>k,piancha,xdpiancha,xk,yk=newtonxz(3,0.5,1e-4,1e-4,100)运行后整理结果得下表.pianchaxdpianchaxkyk0.432306.385790.067700.002940.0665457.310850.001160.000000.001163443.447270.000000.000000.000431.00078-0.0

43、0043-0.000000.000439228.792130.00000-0.000000.011151.000000.011150.000010.01112356.869170.000030.000000.000031638.927030.000000.00000迭代次数k=8,精确到"1CT4的根的近似值是xk=1.9003e-008,其函数值是yk=4.4409e-016,piancha=3.1145e-005,xdpiancha=1638.9270,是二阶收敛(平方收敛).可见，求重根的修正牛顿切线公式比牛顿切线法收敛速度快得多.(二)未知方程根的重数(选学)未知方程f(x)

44、=0根的重数为用，求重根的修正牛顿切线法的MATLAB主程现提供名为newtonxz1.m的M文件functionk,piancha,xdpiancha,xk,yk=newtonxz1(x0,tol,ftol,gxmax)x(1)=x0;fori=1:gxmaxu(i尸fnq(x(i)/dfnq(x(i);du(i)=1-fnq(x(i)*ddfnq(x(i)/(dfnq(x(i)A2+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=f

45、nq(x(i);if(piancha<tol)|(xdpiancha<tol)&(abs(yk)<ftol)k=i-1;xk=x(i);(i-1)pianchaxdpianchaxkykreturn;endendifi>gxmaxdisp(请注意：迭代次数超过给定的最大值gxmax.')k=i-1;xk=x(i);yk=fnq(x(i);(i-1)pianchaxdpianchaxkyk;return;end例12用未知重数的求重根的修正牛顿切线法，求方程在区间0上的根,使精确到二,，L解在MATLAB工作窗口输入程序>>k,piancha,

46、xdpiancha,xk,yk=newtonxz1(0.5,1e-4,1e-4,100)运行后整理结果得下表.kpianchaxdpianchaxkyk110.5992316.03857r0.09923-0.00818120.09698：43.0548510.00225-0.0000030.002251778.32784-0.000000401。10.000000六、割线法及其MATLAB程序1.割线法的MATLAB程序割线法的MATLAB主程序现提供名为gexian.m的M文件：不得用于商业用途仅供个人参考functionk,piancha,xdpiancha,xk,yk=gexian(x0

47、1,x02,tol,ftol,gxmax)x(1)=x01;x(2)=x02;fori=2:gxmaxu(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)pianchaxdpianchaxkykif(abs(yk)<ftol)&(piancha<tol)|(xdpiancha<tol)k=i-2

48、;xk=x(i);yk=fnq(x(i);(i-2)pianchaxdpianchaxkyk；return;endendifi>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax.')k=i-2;xk=x(i);yk=fnq(x(i);return;end七、抛物线法及其MATLAB程序1抛物线法的MATLAB程序抛物线法的MATLAB主程序functionk,piancha,xdpiancha,xk,yk=paowu(x1,x2,x3,tol,ftol,gxmax)X=x1,x2,x3;fori=1:gxmaxh0=X(1)-X(3);h1=X(2)-X(

49、3);u0=fnq(X(1)-fnq(X(3);u1=fnq(X(2)-fnq(X(3);c=fnq(X(3);fenmu=h0A2*hi-h1A2*h0;afenzi=u0*hi-u1*h0;bfenzi=u1*h0A2-u0*h1A2;a=afenzi/fenmu;b=bfenzi/fenmu;deta=bA2-4*a*c;cha=abs(X(3)-X(i),abs(X(3)-X(2),abs(X(2)-X(i);piancha=min(cha);xdpiancha=piancha/(abs(X(3)+eps);ifdeta<0disp('提示：根的判别式值为负数)detak

50、f=sqrt(deta);elseifdeta<=0disp('提示：根的判别式值为零.'),detakf=0;elsedisp('提示：根的判别式值为正数)detakf=sqrt(deta);endifb<0detakf=-detakf;endz=-2*c/(b+detakf);xk=X(3)+z;q1=xk-X(1);q2=xk-X(2);q3=xk-X(3);ifabs(q2)<abs(q1)X1=X(2),X(1),X(3);X=X1;Y=fnq(X);endifabs(q3)<abs(q1)X2=X(1),X(3),X(2);X=X2;

51、Y=fnq(X);endX(3)=xk;Y(3)=fnq(X(3);yk=Y(3);ipianchaxdpianchaxkykif(abs(yk)<ftol)&(piancha<tol)|(xdpiancha<tol)k=i;X(3)=xk;Y(3)=fnq(X(3);yk=Y(3);return;endendifi>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax,请重新输入初始值x1,x2,x3')return;endP=i,piancha,xdpiancha,xk,yk'例13用抛物线法求方程/W=e镰-3y=0的一

52、个实根的近似值k使精确到；二一.解先用作图法确定初始值x01,x02和x03,然后在MATLAB工作窗口输入程>>k,piancha,xdpiancha,xk,yk=paowu(-0.4,-0.3,-0.5,1e-4,1e-4,100)运行后输出结果提示：根的判别式值为正数.ans=Columns1through41.000000000000000.0000000.000000-0.562239Column5-0.000099提示：根的判别式值为正数.ans=Columns1through42.000000000000000.0077610.077038-0.365394Colum

53、n5-0.00000000923532提示：根的判别式值为正数.ans=Columns1through43.000000000000000.0000450.000090-0.082059Column50.00000000000000k=3piancha=1.5282676e-005xdpiancha=4.9938760e-005xk=-0.082059yk=2.9250313e-016即，迭代k=3次，得到精度为的根的近似值xk=-0.3906,其函数值为yk=2.2204e-016,xk的相邻最小偏差为piancha=1.712e-005,其相对误差为xdpiancha=4.3830e-00

54、5.八、求解非线性方程组的牛顿法及其MATLAB程序(选学)1.求解二元非线性方程组的牛顿法及其MATLAB程序解二元非线性方程组的牛顿法的MATLAB主程序functionci,D,danfan,xddf,hanfan,Xk,Yk=newtonzu2(X,tol,ftol,ngmax)Y=Z(X);fori=1:ngmaxdY=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)+