科学工程计算与matlab编程3

第三章微积分问题的计算机求解•微积分问题的解析解•函数的级数展开与级数求和问题求解•数值微分•数值积分问题•曲线积分与曲面积分的计算3.1微积分问题的解析解3.1.1极限问题的解析解•单变量函数的极限–格式1：L=limit(fun,x,x0)–格式2：L=limit(fun,x,x0,‘left’或‘right’)•例:试求解极限问题symsxab;f=x*(1+a/x)^x*sin(b/x);%输入的函数为符号的函数，所以其中的加减乘除都是通常的加减乘除。L=limit(f,x,inf)L=exp(a)*bInf表示无穷大•例:求解单边极限问题symsx;limit((exp(x^3)-1)/(1-cos(sqrt(x-sin(x)))),x,0,'right')ans=12•在(-0.1,0.1)区间绘制出函数曲线：x=-0.1:0.001:0.1;y=(exp(x.^3)-1)./(1-cos(sqrt(x-sin(x))));Warning:Dividebyzero.(TypewarningoffMATLAB:divideByZerotosuppressthiswarning.)plot(x,y,'-',[0],[12],‘o’)%表示将[0，12这个点放o](Matlat7.0Type“warningofflast”)•多变量函数的极限：–格式：L1=limit(limit(f,x,x0),y,y0)或L1=limit(limit(f,y,y0),x,x0)如果x0或y0不是确定的值，而是另一个变量的函数，如x-g(y)，则上述的极限求取顺序不能交换。•例：求出二元函数极限值symsxya;f=exp(-1/(y^2+x^2))…*sin(x)^2/x^2*(1+1/y^2)^(x+a^2*y^2);L=limit(limit(f,x,1/sqrt(y)),y,inf)L=exp(a^2)3.1.2函数导数的解析解•函数的导数和高阶导数–格式：y=diff(fun,x)%求导数y=diff(fun,x,n)%求n阶导数•例：一阶导数：symsx;f=sin(x)/(x^2+4*x+3);f1=diff(f,x);pretty(f1)cos(x)sin(x)(2x+4)-----------------------------------222x+4x+3(x+4x+3)原函数及一阶导数图：x1=0:.01:5;%定义一个数组y=subs(f,x,x1);%将x换成x1y1=subs(f1,x,x1);%将y1中的x换成x1plot(x1,y,x1,y1,’:’)更高阶导数：tic,diff(f,x,100);tocelapsed_time=4.6860f1=cos(x)/(x^2+4*x+3)-sin(x)/(x^2+4*x+3)^2*(2*x+4)f1=cos(x)/(x^2+4*x+3)-sin(x)/(x^2+4*x+3)^2*(2*x+4)f1=cos(x)/(x^2+4*x+3)-sin(x)/(x^2+4*x+3)^2*(2*x+4)•原函数4阶导数f4=diff(f,x,4);pretty(f4)2sin(x)cos(x)(2x+4)sin(x)(2x+4)------------+4--------------------12-----------------22223x+4x+3(x+4x+3)(x+4x+3)3sin(x)cos(x)(2x+4)cos(x)(2x+4)+12----------------24-----------------+48----------------222423(x+4x+3)(x+4x+3)(x+4x+3)42sin(x)(2x+4)sin(x)(2x+4)sin(x)+24------------------72-----------------+24---------------252423(x+4x+3)(x+4x+3)(x+4x+3)•多元函数的偏导：–格式：f=diff(diff(f,x,m),y,n)或f=diff(diff(f,y,n),x,m)•例：求其偏导数并用图表示。symsxyz=(x^2-2*x)*exp(-x^2-y^2-x*y);zx=simple(diff(z,x))zx=-exp(-x^2-y^2-x*y)*(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y)zy=diff(z,y)zy=(x^2-2*x)*(-2*y-x)*exp(-x^2-y^2-x*y)•直接绘制三维曲面[x,y]=meshgrid(-3:.2:3,-2:.2:2);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);surf(x,y,z),axis([-33-22-0.71.5])contour(x,y,z,30),holdon%绘制等值线zx=-exp(-x.^2-y.^2-x.*y).*(-2*x+2+2*x.^3+x.^2.*y-4*x.^2-2*x.*y);zy=-x.*(x-2).*(2*y+x).*exp(-x.^2-y.^2-x.*y);%偏导的数值解quiver(x,y,zx,zy)%绘制引力线•例symsxyz;f=sin(x^2*y)*exp(-x^2*y-z^2);df=diff(diff(diff(f,x,2),y),z);df=simple(df);pretty(df)22222-4zexp(-xy-z)(cos(xy)-10cos(xy)yx+42422422sin(xy)xy+4cos(xy)xy-sin(xy))•多元函数的Jacobi矩阵:–格式：J=jacobian(Y,X)%注意其顺序其中，X是自变量构成的向量，Y是由各个函数构成的向量。•例：试推导其Jacobi矩阵symsrthetaphi;x=r*sin(theta)*cos(phi);y=r*sin(theta)*sin(phi);z=r*cos(theta);J=jacobian([x;y;z],[rthetaphi])J=[sin(theta)*cos(phi),r*cos(theta)*cos(phi),-r*sin(theta)*sin(phi)][sin(theta)*sin(phi),r*cos(theta)*sin(phi),r*sin(theta)*cos(phi)][cos(theta),-r*sin(theta),0]•隐函数的偏导数：•格式：F=-diff(f,xj)/diff(f,xi)0ijjijxdfffdxxxx•例：•要求的是y对x的偏导数，而不是z对x或者Y的偏导数，所以是隐函数的求导。symsxy;f=(x^2-2*x)*exp(-x^2-y^2-x*y);pretty(-simple(diff(f,x)/diff(f,y)))322-2x+2+2x+xy-4x-2xy------------------------------------------x(x-2)(2y+x)3.1.3积分问题的解析解•不定积分的推导：–格式：F=int(fun,x)•例：用diff()函数求其一阶导数，再积分，检验是否可以得出一致的结果。symsx;y=sin(x)/(x^2+4*x+3);y1=diff(y);y0=int(y1);pretty(y0)%对导数积分sin(x)sin(x)-1/2------+1/2------x+3x+1•对原函数求4阶导数，再对结果进行4次积分y4=diff(y,4);y0=int(int(int(int(y4))));pretty(simple(y0))sin(x)------------2x+4x+3•例:证明symsax;f=simple(int(x^3*cos(a*x)^2,x))f=1/16*(4*a^3*x^3*sin(2*a*x)+2*a^4*x^4+6*a^2*x^2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a^4f1=x^4/8+(x^3/(4*a)-3*x/(8*a^3))*sin(2*a*x)+...(3*x^2/(8*a^2)-3/(16*a^4))*cos(2*a*x);simple(f-f1)%求两个结果的差ans=-3/16/a^4%得到的结果是常数，所以得到证明。•定积分与无穷积分计算:–格式：I=int(f,x,a,b)–格式：I=int(f,x,a,inf)•例：symsx;I1=int(exp(-x^2/2),x,0,1.5)％无解I1=1/2*erf(3/4*2^(1/2))*2^(1/2)*pi^(1/2)vpa(I1,70)ans=1.085853317666016569702419076542265042534236293532156326729917229308528I2=int(exp(-x^2/2),x,0,inf)I2=1/2*2^(1/2)*pi^(1/2)2/2()xfxe202()xterfxedt•多重积分问题的MATLAB求解•例：symsxyz;f0=-4*z*exp(-x^2*y-z^2)*(cos(x^2*y)-10*cos(x^2*y)*y*x^2+...4*sin(x^2*y)*x^4*y^2+4*cos(x^2*y)*x^4*y^2-sin(x^2*y));f1=int(f0,z);f1=int(f1,y);f1=int(f1,x);f1=simple(int(f1,x))f1=exp(-x^2*y-z^2)*sin(x^2*y)f2=int(f0,z);f2=int(f2,x);f2=int(f2,x);f2=simple(int(f2,y))f2=2*exp(-x^2*y-z^2)*tan(1/2*x^2*y)/(1+tan(1/2*x^2*y)^2)simple(f1-f2)ans=0顺序的改变使化简结果不同于原函数，但其误差为0，表明二者实际完全一致。这是由于积分顺序不同，得不出实际的最简形式。•例：symsxyzint(int(int(4*x*z*exp(-x^2*y-z^2),x,0,1),y,0,pi),z,0,pi)ans=(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2))*pi^2*hypergeom([1],[2],-pi^2)%Ei指数积分函数，为试错函数，没有解析解vpa(ans,60)ans=3.108079402085412722834614647671385210191423063170218634835881EulerEi(,)Ei(,),ztneulergammanznzetdt为常数；为指数积分，即＝该函数无解析解，但可求其数值解。3.2函数的级数展开与级数求和问题求解•3.2.1Taylor幂级数展开•3.2.2Fourier级数展开•3.2.3级数求和的计算3.2.1Taylor幂级数展开3.2.1.1单变量函数的Taylor幂级数展开例:symsx;f=sin(x)/(x^2+4*x+3);y1=taylor(f,x,9);pretty(y1)22333444087530676515273738645981/3x-4/9x+--x-----x+------x-------x+----------x----------x5481972072901224720918540taylor(f,x,9,2)ans=1/15*sin(2)+(1/15*cos(2)-8/225*sin(2))*