MATLAB数学建模与仿真（第2版·微课视频版）程序代码 第8章-Matlab在高等数学的应用

【例8-1】QUOTEQUOTE的几何解释。下面的示意图中共描了40个点，取如图8-1所示，N可以取2或大于2的正整数，即当n>2时，QUOTEQUOTEn+(-1)n-1nn+(-1)n-1n的值落在直线y=0.5,y=1.5之间；取时，如图8-2所示，N可以取20，即当n>20时，QUOTEQUOTEn+(-1)n-1nn+(-1)n-1n的值落在直线y=0.9，y=1.1之间。可见下面是作示意图8-1，8-2的源程序，可以修改点数大于40，观察数列趋于1的变化，修改的大小，观察N随的变化，加深对数列极限的理解。clfsubplot(1,2,1)holdongridonn=1:40;%描40个点m=1+(-1).^(n-1)./n;plot(n,m,'.')fplot(@(t)0.5*ones(size(t)),[0,40])%取0.5fplot(@(t)1.5*ones(size(t)),[0,40])axis([0,40,0,2])subplot(1,2,2)holdongridonplot(n,m,'.')fplot(@(t)0.9*ones(size(t)),[0,40])%取0.1fplot(@(t)1.1*ones(size(t)),[0,40])图8-1当时图8-2当时【例8-2】函数QUOTEQUOTEfx=x-1x<00x=0x+1x>0fx=x-1x<00x=0x+1x>0，当QUOTEQUOTEf(x)的极限不存在。观察图7-3，当时，QUOTEQUOTEf(x)的左极限为－1，右极限为1，QUOTEQUOTE，所QUOTE不存在,x=0是函数f(x)的跳跃间断点。图8-3当QUOTEQUOTEf(x)的极限作图的源程序如下：holdonfplot(@(x)x-1,[-5,0])fplot(@(x)x+1,[0,5])plot(0,0,'.')axis([-5,5,-6,6])gridonplot(0,1,'o')plot(0,-1,'o')（2）函数f(x)当QUOTEQUOTE时极限为A的几何解释：任意给定一正数，作平行于x轴的两条直线和，介于这两条直线之间是一横条区域。 【例8-3】QUOTEQUOTE的几何解释。图8-4当时图8-5当时取QUOTEQUOTE蔚=1蔚=1时如图8-4所示，QUOTEQUOTE或更小，即当1.3<x<1.5时，函数的图形落在两直线y=0.44,y=2.44之间；取QUOTEQUOTE时如图8-5所示，QUOTEQUOTE，即当1.39<x<1.41时，函数的图形落在两直线y=1.34，y=1.54之间。可见QUOTEQUOTE未未的取值依赖于QUOTEQUOTE取值，且不唯一。下面是作示意图8-4，8-5的源程序，修改QUOTEQUOTE的大小，观察QUOTEQUOTE未未随QUOTEQUOTE的变化，加深对自变量趋于有限数时函数极限的理解。clfsubplot(1,2,1)holdongridonfplot(@(x)-(x-3).^2+4,[1,2])fplot(@(t)0.44*ones(size(t)),[1,2])%取1时fplot(@(t)2.44*ones(size(t)),[1,2])axis([1,2,-1,3])plot(1.4,-1.6*1.6+4,'o')subplot(1,2,2)holdongridonfplot(@(x)-(x-3).^2+4,[1,2])fplot(@(t)1.34*ones(size(t)),[1,2])%取0.1时fplot(@(t)1.54*ones(size(t)),[1,2])axis([1.38,1.42,1.2,1.7])plot(1.4,-1.6*1.6+4,'o')（3）自变量趋于无穷大时函数的极限如果QUOTEQUOTE则直线y=c是函数y=f(x)的图形的水平渐近线。（4）函数f(x)当时的极限为A的几何解释：任意给定一正数QUOTE蔚蔚，作平行于x轴的两条直线和，则总有一个正数X存在，使得当x<-X或x>X时，函数y=f(x)的图形位于这两直线之间。【例8-4】QUOTE的几何解释。取QUOTEQUOTE如图8-6所示，X可以取50或更大，即当x>X时，函数的图形落在两直线y=-0.1,y=0.1之间；取QUOTEQUOTE如图8-7所示，X可以取2000或更大，即当x<-X时，函数的图形落在两直线y=-0.001,y=0.001之间。可见N的取值依赖于QUOTEQUOTE蔚蔚取值，且不唯一。图8-6当时图8-7当时下面是绘制图8-6，8-7的源程序，修改QUOTEQUOTE蔚蔚的大小，观察X随的变化，加深对自变量趋于无穷大时函数极限的理解。clfsubplot(1,2,1)gridonholdonfplot(@(x)1./x,[2,200])fplot(@(t)0.1*ones(size(t)),[2,200])%取0.1时fplot(@(t)-0.1*ones(size(t)),[2,200])subplot(1,2,2)gridonholdonfplot(@(x)1./x,[-10000,-100])fplot(@(t)-0.001*ones(size(t)),[-10000,-100])%取0.001时fplot(@(t)0.001*ones(size(t)),[-10000,-100])【例8-5】QUOTEQUOTE图8-8x趋于0时图8-9x=0时当x趋于0时,如图8-8所示,QUOTEQUOTEsin1xsin1x的值在-1与1之间来回波动,有界,但没有极限,x=0是函数的振荡间断点.如图8-9所示,的值不断振荡,但向0趋近。生成图8-8，8-9的程序如下：clfsubplot(1,2,1)fplot(@(x)sin(1./x),[-0.0001,0.0001])gridonsubplot(1,2,2)fplot(@(x)sin(1./x).*x,[-0.0001,0.0001])gridon【例8-6】求下列极限（1）QUOTEQUOTE在MATLAB的命令窗口输入：symsxlimit((cos(x)-exp(-x^2/2))/x^4,x,0)运行结果为ans=-1/12理论上用洛必达法则或泰勒公式计算该极限：方法1：方法2：（2）(自变量趋于无穷大,带参数t)在MATLAB的命令窗口输入：symsxtlimit((1+2*t/x)^(3*x),x,inf)运行结果为ans=exp(6*t)理论上用重要极限计算：（3）求右极限。在MATLAB的命令窗口输入：symsxlimit(1/x,x,0,’right’)运行结果为ans=inf【例8-7】在（1，1）处切线方程为，如图8-12所示。函数在处连续是函数在处可导的必要条件。生生成图8-12的程序如下:clfholdon fplot(@(x)x.^2,[0,2])fplot(@(x)2.*x-1,[0,2])plot(1,1,'*')图8-12的切线方程下面的程序可演示函数在(1,1)处的切线的生成过程。symsxfplot(@(x)x.*x,[0,2],'k')fori=1:9;x1=2-0.1*i;y1=x1*x1;pauseeplot(@(x)1+(y1-1)/(x1-1)*(x-1),[0.8,2])endfplot(@(x)2.*x-1,[0.8,2],'r')【例8-8】函数在x=0处连续,但在x=0处不可导，如图8-13所示。图8-13图形生成图8-13的程序为:clf holdon fplot(@(x)abs(x),[-1,1]) 【例8-9】已知,求。在MATLAB的命令窗口输入如下命令序列:symsxyzz=x^2*sin(2*y);diff(z,x)diff(z,x,2)diff(diff(z,x),y)输出结果：ans=2*x*sin(2*y)ans=2*sin(2*y)ans=4*x*cos(2*y)已知,求（复合函数求导偏导数）。在MATLAB的命令窗口输入如下命令序列:symsxyzuz=x^2+y^2;u=(x-y)^z;diff(u,x)diff(u,y,2)(x-y)^(x^2+y^2)*(2*log(x-y)-4*y/(x-y)-(x^2+y^2)/(x-y)^2)diff(diff(u,x),y)输出结果：ans=(x-y)^(x^2+y^2)*(2*x*log(x-y)+(x^2+y^2)/(x-y))ans=(x-y)^(x^2+y^2)*(2*y*log(x-y)-(x^2+y^2)/(x-y))^2+ans=(x-y)^(x^2+y^2)*(2*y*log(x-y)-(x^2+y^2)/(x-y))*(2*x*log(x-y)+(x^2+y^2)/(x-y))+(x-y)^(x^2+y^2)*(-2*x/(x-y)+2*y/(x-y)+(x^2+y^2)/(x-y)^2)【例8-10】已知函数,用MATLAB软件计算。求函数f(x)的一阶，二阶导数，并画出它们相应的曲线。观察函数的单调区间，凹凸区间，以及极值点和拐点。在MATLAB的命令窗口输入如下命令序列：symsxy=x^3-x^2-x+1d1=diff(y,x)%求一阶导数d2=diff(d1,x)%求二阶导数clfsubplot(1,1,1)holdongridonezplot(y,[-22])gtext('f(x)')ezplot(d1,[-2,2])gtext('f’(x)')ezplot(d2,[-2,2])gtext('f’’(x)')title('导数的应用')gtext('o')gtext('(x1,y1)')gtext('o')gtext('(x2,y2)')gtext('o')gtext('(x3,y3)')f1=char(d1)x1=fzero(f1,0)%求一阶导函数在x＝0附近的零点x2=fzero(f1,1)%求一阶导函数在x＝1附近的零点f2=char(d2)x3=fzero(f2,0)%求二阶导数在x＝0附近的零点图8-14曲线从图8-14中可以清楚地看到：和为极值点，对应的一阶导数为零，为拐点，对应的二阶导数为0，在单调上升区间函数的一阶导数大于零，在单调下降区间函数的一阶导数小于0，在极大值点处二阶导数小于0，在极小值点处二阶导数大于0，在凸区间函数的二阶导数小于0，在凹区间函数的二阶导数大于0。【例8-11】求函数的极值，并作图。在MATLAB的命令窗口输入如下命令序列：symsxf=2.*x.^3-6.*x.^2-18.*x+7;xmin=fminbnd('2.*x.^3-6.*x.^2-18.*x+7',-5,5)x=xmin;miny3=subs(f)a31='-2.*x.^3+6.*x.^2+18.*x-7';xmax=fminbnd(a31,-5,5)x=xmax;maxy3=subs(f)fplot(@(x)2.*x.^3-6.*x.^2-18.*x+7,[-55])gridon执行结果：xmin=3.0000%在x＝3处取极小值miny3=-47.0000%极小值为-47xmax=-1.0000%在x＝－1处取极大值maxy3=17.0000％极大值为17图8-15的极值图【例8-12】求函数在附近的极小值。具体操作步骤：打开MATLAB软件，单击菜单file，new，m-file，进入M文件编辑窗口。在编辑窗口输入：functiony=f1(x)y=(x(1)*x(1)-4*x(2))^2+120*(1-2*x(2))^2;单击菜单file→save命令存盘，命名为f1.m。回到MATLAB命令窗口，输入以下命令计算函数极小值点：d1=fminsearch('f1',[-2,2])%计算函数的极小值点。执行结果：d1=-1.41420.5000再输入f1(d1)计算极小值f1(d1)执行结果：ans=9.7459e-009【例8-13】计算不定积分。QUOTEQUOTEexcos2xdxexcos2xsymsx;int(exp(x)*cos(2*x),x)结果： ans=(exp(x)*(cos(2*x)+2*sin(2*x)))/5【例8-14】计算不定积分QUOTE。symsx;int(1/(x^4*sqrt(1+x^2)))结果：ans=(x^2+1)^(1/2)*(2/(3*x)-1/(3*x^3))【例8-15】求由抛物线和直线y=-x+4所围成的面积。首先画出图形，如图7-16。图8-16和y=-x+4图形x=0:0.1:9;plot(x,-x+4,'b',x,sqrt(2*x),'r',x,-sqrt(2*x),'r')%结果为8-16图求解方程组，得到两曲线交点symsxyeqns=[y^2-2*x==0,y+x-4==0];S=solve(eqns,[xy])结果：S=包含以下字段的struct:x:[2×1sym]y:[2×1sym]求x，y的值S.xS.y结果：ans=28ans=2-4以y为积分变量求面积int(-y+4-y^2/2,y,-4,2)结果：ans= 18【例8-16】用MATLAB软件，计算下列不定积分。（1）在MATLAB的命令窗口输入如下命令：symsxint(x^3*exp(-x^2),x)执行结果：ans=-(exp(-x^2)*(x^2+1))/2理论推导：(2)在MATLAB的命令窗口输入如下命令：symsxy=[sin(x),x^3;x*exp(x),tan(x)]；int(y)执行结果：ans=[-cos(x),x^4/4][exp(x)*(x-1),-log(cos(x))]（3）求不定积分，并取a=2，b=3绘制其函数图像并说明不定积分的几何意义。试探讨参数a和b对积分曲线的影响。在MATLAB输入程序：symsxaC%定义符号F=int((sin(a*x/2))^2+(x^3)/35)%计算不定式F=2/a*(-1/2*cos(1/2*a*x)*sin(1/2*a*x)+1/4*a*x)+1/140*x^4;%不定积分的符号解y=simplify(F)+Cy=1/140*(-70*sin(a*x)+70*a*x+x^4*a)/a+cx=-2*pi:0.01:2*pi;a=2;forC=-28:28y=1/140*(-70*sin(a*x)+70*a*x+x.^4*a)/a+C;plot(x,y)holdonendgridholdoffaxis([-2*pi,2*pi,-8,8])xlabel('x')ylabel('y')title（'函数y=sin(a*x/2)^+(x^3)/35的积分曲线'）legend（'函数y=sin(a*x/2)^+(x^3)/35的积分曲线'）得到图像7-17图8-17函数图像【例8-17】求由圆r=3cos和双纽线1+r=cos所围图形的面积。首先在极坐标系下画出两曲线的图形，如图7-18。Th=0:0.5:2*pi;th=0:0.05:2*pi;r1=3*cos(th);r2=1+cos(th);polar(th,r1,'b');holdon;polar(th,r2,'r');holdoff%由对称性,求另一个交点(x/3.3/2)或(pi/3,3/2)s1=int(1/2*r2.^2,th,0,pi/3);s2=int(1/2*r1.^2,th,pi/3,pi/2);S=2*(s1+s2)结果：S=5/4*pi图8-18r=3cos和1+r=cos图形示意图【例8-18】利用MATLAB软件找出满足定积分中值定理的点,使得。在MATLAB的命令窗口输入如下命令序列：symsxy=sqrt(1-x^2);zhi=int(y,0,1);%计算z=y-zhi;zf=char(z);fzero(zf,0.5)%求满足的运行结果：ans=0.6190【例8-19】用MATLAB软件求下列定积分。（1）（2）在MATLAB的命令窗口输入如下命令序列：symsx;y=log(x)*x^(-0.5);int(y,1,4)运行结果:ans=log(256)-4在MATLAB的命令窗口输入如下命令序列：symsx;y=(x*(1+x)^3)^(-0.5);int(y,0,+inf)运行结果：ans=2【例8-20】求QUOTEQUOTEddx0xsinttdtddxsymstx;y=sin(t)/t;diff(int(y,t,0,x),x)运行结果：ans=sin(x)/x【例8-21】求下列极限。symsxt;f=cos(t^2);int(f,t,sin(x),0);f1=diff(int(f,t,sin(x),0),x)f2=f1/1limit(f2)运行结果：f1=-cos(sin(x)^2)*cos(x)f2=-cos(sin(x)^2)*cos(x)ans= -1【例8-23】用多种方法计算反常积分的近似值。（1）方法1：调用MATLAB函数计算，在MATLAB的命令窗口输入：quad('sin(x)./x',0,1,0.001)运行结果：ans=0.946083070076534（2）方法2：数值积分法，等分区间数为10，100，500。在M文件编辑窗口输入上述函数，分别存盘为zjx.m,djx.m,rjx.m,txf.m,smp.m在MATLAB的命令窗口输入一下命令序列：formatlongdigits(20)symsxy=sin(x)/x;zjx(y,0,1,10)zjx(y,0,1,100)zjx(y,0,1,500)rjx(y,0,1,10)rjx(y,0,1,100)rjx(y,0,1,500)djx(y,0,1,10)djx(y,0,1,100)djx(y,0,1,500)smp(y,0,1,10)smp(y,0,1,100)smp(y,0,1,500)表8-1数值积分法结果n10100500左矩形法0.9537585226265100.9468732057016930.946241498992812中矩形法0.9462085788431450.9460843252388310.946083120561967右矩形法0.9379056211073000.9452879155497720.945924440962428梯形法0.945832071866910.946080560625730.94608296997762辛甫生公式法0.9460830765177320.9460830703677980.946083070367185运行结果如下表8-1。结果精度差距很大，与积分的准确值I=0.946083070076534比较，辛甫生公式收敛的速度最高，而且精度高。（3）方法3：近似函数法，分别用5阶，10阶，20阶泰勒展开式进行近似计算。在MATLAB的命令窗口输入：ty=taylor(sin(x)/x,x,'Order',5);%5阶泰勒展开式，默认x=0int(ty,0,1)ty=taylor(sin(x)/x,x,'Order',10);%10阶泰勒展开式，默认x=0int(ty,0,1)ty=taylor(sin(x)/x,x,'Order',20);%20阶泰勒展开式，默认x=0int(ty,0,1)表8-2近似函数法结果n51020近似函数法0.9461111111111110.9460830726323450.946083070367183精度相较之前增加一位精度相较之前增加一位从表8-2看出，10阶泰勒展开式的近似结果精度已经很高。【例8-23】求由抛物线QUOTEQUOTEy2y2与所围图形的面积A。在MATLAB的命令窗口输入一下命令序列，画出积分区域的图形：y=linspace(-1,1,60);x1=5*y.^2;x2=1+y.^2;plot(x1,y,x2,y)运行结果如图8-20所示。图8-20QUOTEQUOTEy2y2与函数图形观察曲线，再计算面积：symsyf=(1+y^2)-5*y^2;A=int(f,y,-0.5,0.5)运行结果：A=2/3即所求平面图形的面积为2/3.【例8-24】计算,其中D为直线围成区域。具体步骤如下。划定积分区域：y1=2*x;y2=x/2;y3=12-x;ezplot(y1,[-2,12])holdonezplot(y2,[-2,12])ezplot(y3,[-2,12])title('积分区域')结果如图8-21，三条直线相交所围区域即为积分区域。图8-21围成区域确定交点的横坐标：xa=fzero('2*x-x/2',0)xb=fzero('2*x-12+x',4)xc=fzero('12-x-x/2',8)结果为:xa=0xb=4xc=8化二重积分为累次积分。在MATLAB的命令窗口输入以下命令序列：symsxyzz=x^2/y^2;dx1=int(z,y,x/2,2*x);j1=int(dx1,0,4);dx2=int(z,y,x/2,12-x);j2=int(dx2,4,8);jf=j1+j2jf=运行结果：132-144*log(2)【例8-25】计算数值积分，可将此二重积分转化为累次积分QUOTEQUOTE-11-1+x21+x2(1+x+y)dy输入MATLAB代码为：clear;symsxy;iy=int(1+x+y,y,-sqrt(1-x^2),sqrt(1-x^2));int(iy,x,-1,1)运行结果：ans=pi【例8-26】设。首先输入以下代码：forn=1:25a=1;form=1:na=a*10/m;endplot(n,a,'*')holdonend得到的输出如图8-22所示。从散点图可见QUOTEQUOTEanan的变化趋势。图8-21QUOTEQUOTEanan的变化趋势输入以下命令：symsn;symsum(10^n/gamma(n+1),n,1,inf)输出：ans=exp(10)-1【例8-27】判别级数的收敛性。因为所以故级数收敛。（d）根值审敛法。设是正项级数，并且，则=1\*GB3①当时，级数收敛；=2\*GB3②当(或)时，级数发散；=3\*GB3③当时，级数可能收敛，也可能发散。【例8-29】求QUOTEQUOTE的收敛性与和函数。输入：clear;symsnxa1=4^(2*n)*(x-3)^n/(n+1);a2=subs(a1,n,n+1);p=limit(a2/a1,n,inf)输出为： p= 16*x-48注意,这里对a2和a1都没有加绝对值。因此上式的绝对值小于1时,幂级数收敛,大于1时发散。为了求出收敛区间的端点,输入：x1=solve(16*x-48==1)x2=solve(16*x-48==-1)输出为: x1=49/16x2=47/16由此可知QUOTEQUOTE4716<x<49164716<x<49为了判断端点的敛散性,输入：simplify(subs(a1,'x',49/16))得到x为右端点时幂级数的一般项为： ans= 1/(n+1)因此当x=49/16时发散。再输入： simplify(subs(a1,'x',47/16))输出结果为： ans=(-1)^n/(n+1)因此当x=47/16时,级数收敛。【例8-30】求函数在x=1处的3阶泰勒展式。symsxtaylor(exp(x),x,1,'Order',3)执行后得到:ans=exp(1)+exp(1)*(x-1)+(exp(1)*(x-1)^2)/2【例8-31】求函数的不同阶数的麦克劳林展式,并作图观察不同阶数展式对函数的近似程度,计算的近似值。symsxsin3sin5sin7sin9sin3=taylor(sin(x),'Order',3);sin5=taylor(sin(x),'Order',5);sin7=taylor(sin(x),'Order',7);sin9=taylor(sin(x),'Order',9);%sin9=x-1/6*x^3+1/120*x^5-1/5040*x^7clfholdonezplot(sin(x))gtext('sinx')ezplot(sin3)gtext('3阶')ezplot(sin5)gtext('5阶')ezplot(sin7)gtext('7阶')ezplot(sin9)gtext('9阶')title('y=sinx及其泰勒展式曲线')sin(pi/8)x=pi/8; zhi3=eval(sin3)zhi5=eval(sin5)zhi7=eval(sin7)zhi9=eval(sin9)图8-23的麦克劳林展式图的准确值：0.38268343236508978178。不同阶数的麦克劳林展开式的近似值如表8-3所示表8-3准确值表n3579sin（x）0.392699081698724154810.382605892675189018790.382683717505507559050.3826834317539123953表8-3准确值表由图形8-23和表8-3可知，阶数越高，近似精度越高。【例8-32】设g(x)是QUOTEQUOTE为周期函数，它在QUOTEQUOTE的表达式，将g(x)展开成傅里叶级数。因为g(x)是奇函数,所以它的傅里叶展开式中只含正弦项。计算系数： symskx bk=2*int(sin(k*x),x,0,pi)/pi输出 bk=(4*sin((pi*k)/2)^2)/(k*pi)再输入:clear;f='sign(sin(x))';x=-3*pi:0.1:3*pi;y1=eval(f);plot(x,y1,'r')pauseholdonforn=3:2:9fork=1:nbk=(4*sin((pi*k)/2)^2)/(pi*k);s(k,:)=bk*sin(k*x);ends=sum(s);plot(x,s)pauseholdonend运行结果如图8-24所示。可以看到,n越大,g(x)的傅里叶级数的前n项与g(x)越接近。图8-24g(x)展开傅里叶级数图【例8-33】计算e的近似值，精确到。e的值就是的展开式在处的函数值，即,其误差为．要精确到，须，故取，得。用MATLAB创建函数edejs.m文件，wch代表误差。functionzhi=edejs(wch)symsxsn=5;s=taylor(exp(x),n);x=1;zhi1=eval(s);n=n+1;symsxss=taylor(exp(x),n);x=1;zhi2=eval(s);while(abs(zhi1-zhi2)>wch)zhi1=zhi2;n=n+1;symsxss=taylor(exp(x),n);x=1;zhi2=eval(s);endzhi=zhi2;在MATLAB命令窗口输入formatlongjszhi=edejs(0.0000000001)执行结果为:jszhi=2.71828182845823【例8-34】计算ln2的近似值，要求误差不超过0.0001。方法1：在展开式：中，设x=1，得为了保证误差不超过，须取n＝10000项进行计算。这样做计算量太大了，我们必须用收敛较快的级数代替它。方法2：在的展开式中将x换成，得：两式相减，得到不含偶次幂的展开式：令，得。如果取前四项的和作为ln2的近似值，其误差。于是有。用MATLAB创建函数ln2js1.m文件实现方法1对ln2进行近似，函数ln2js2.m文件实现方法2对ln2进行近似，n代表泰勒展开阶数。functionzhi=ln2js1(n)%函数ln2js1.m文件symsxstt=log(1+x);s=taylor(t,’Order’,n);x=1;zhi=eval(s);functionzhi=ln2js2(n)%函数ln2js2.m文件symsxstt=log((1+x)/(1-x));s=taylor(t,’Order’,n);x=1/3;zhi=eval(s);在MATLAB命令窗口输入以下命令序列：formatlongln2js1(5)ln2js1(10)ln2js1(50)ln2js1(1000)log(2)ln2js2(5)ln2js2(10)ln2js2(50)执行结果比较如下:ln2的准确值为:0.693147180559945表8-4近似值表n方法5105010000方法10.7833333333333330.7833333333333330.7833333333333330.783333333333333方法20.6931471805599450.6930041152263370.6930041152263370.693004115226337观察表8-4数据得知：方法2收敛速度比方法1快得多。因此在近似计算中应恰当地选取近似公式，保证收敛的同时，注意收敛的速度。【例8-35】用MATLAB函数编程、二分法、切线法三种方法求方程的实根的近似值，使误差不超过。令,显然f(x)在内连续。图8-25f(x)的函数图像因为，故f(x)在内单调递增，至多有一个实根。由，知在[0,1]内有唯一的实根。取a=0,b=1,[0,1]即是一个隔离区间。先画出函数f(x)的图形,如图8-25，在MATLAB的命令窗口输入如下命令：f=@(x)x^3+1.1*x^2+0.9*x-1.4'fplot(@(x)x.^3+1.1.*x.^2+0.9.*x-1.4,[0,1])gridon方法1：在MATLAB的命令窗口输入如下命令：f=@(x)x^3+1.1*x^2+0.9*x-1.4;fzero(f,1)运行结果：ans=0.670657310725810方法2：二分法求方程的近似解。求方程f＝0的近似解的函数erff定义如下：functiony=erff(f,a,b,wch)％其中f是符号函数，a,b是隔离区间的左右端点值，wch是要求的误差，函数返回近似解和迭代次数。x=a;fa=eval(f);x=(a+b)/2;fx=eval(f);n=0;while(b-a)>wchn=n+1;if(fx*fa>0)a=x;elseb=x;endx=a;fa=eval(f);x=(a+b)/2;fx=eval(f);endy(1)=a;y(2)=n;在MATLAB的命令窗口输入如下命令：f='x^3+1.1*x^2+0.9*x-1.4';a=0;b=1;wch=0.0001;erf(f,a,b,wch)运行结果为：ans=0.67065429687500014.000000000000000二分法迭代14次，近似解为0.670654296875000。方法3：切线法求方程的近似解。求方程f＝0的近似解的函数qxf定义如下：functiony=qxf(f,a,b,wch)％f是符号函数a，b是隔离区间的左右端点值，wch是要％求的误差，函数返回近似解和迭代次数。d1f=diff(f);d2f=diff(f);x=(a+b)/2;f1=eval(d1f);f2=eval(d2f);n=1;iff1*f2>0z(1)=b;elsez(1)=a;endx=z(n);z(n+1)=z(n)-eval(f)/eval(d1f);whileabs(z(n+1)-z(n))>wchn=n+1;x=z(n);z(n+1)=z(n)-eval(f)/eval(d1f);endy(1)=z(n+1);y(2)=n;在MATLAB的命令窗口输入如下命令：f='x^3+1.1*x^2+0.9*x-1.4';a=0;b=1;wch=0.0001;qxf(f,a,b,wch)运行结果为：ans=0.67074.0000切线法计算共迭代4次，近似解为0.6707。比较结果可知切线法的收敛速度快。【例8-36】解常微分方程。symsy(x)eqn=diff(y,x,2)+diff(y,x,1)==x*cos(2*x);y=dsolve(eqn)解得结果为：y=C1+(13*cos(2*x))/100+(4*sin(2*x))/25-(x*cos(2*x))/5+(x*sin(2*x))/10+C2*exp(-x)-1/4原文件没有，后来添加的原文件没有，后来添加的【例8-37】求常微分方程满足的特解。symsy(x)Dy2=diff(y,x,2);Dy=diff(y,x);eqn=diff(y,x,3)-diff(y,x,2)==x;cond=[y(1)==8,Dy(1)==7,Dy2(1)==4];y=dsolve(eqn,cond)解得结果为:y=(5*x)/2+6*exp(-1)*exp(x)-x^2/2-x^3/6+1/6【例8-38】用Euler方法求解常数微分方程初值问题。并将数值解和该问题的解析解比较：y(x)=QUOTEQUOTEx1+x2x1+x2Euler方法的具体格式如下：程序实现：h=0.2;y(1)=0.2;x=0.2:h:3;forn=1:14xn=x(n);yn=y(n);y(n+1)=yn+h*(yn/xn-2*yn*yn);endx0=0.2:h:3;y0=x0./(1+x0.^2);plot(x0,y0,x,y,x,y,'o')结果如下图8-27。图8-27Euler方法h=0.2,xn=nh,(n=0,1,2…,15),f(x,y)=y/x–2y2，计算中取f(0,0)=1。计算结果如下：xn y(xn) yn yn-y(xn) 0.0 0 0 00.2 0.1923 0.2000 0.00770.4 0.3448 0.3840 0.03920.6 0.4412 0.5170 0.07580.8 0.4878 0.5824 0.0946 1.0 0.5000 0.5924 0.09241.2 0.4918 0.5705 0.07871.4 0.4730 0.5354 0.06241.6 0.4494 0.4972 0.0478 1.8 0.4245 0.4605 0.03592.0 0.4000 0.4268 0.02682.2 0.3767 0.3966 0.0199 2.