中科院-matlab在科学计算中的应用3

1第三章微积分问题的计算机求解•微积分问题的解析解•函数的级数展开与级数求和问题求解•数值微分•数值积分问题•曲线积分与曲面积分的计算23.1微积分问题的解析解3.1.1极限问题的解析解•单变量函数的极限–格式1：L=limit(fun,x,x0)–格式2：L=limit(fun,x,x0,‘left’或‘right’)3•例:试求解极限问题symsxab;f=x*(1+a/x)^x*sin(b/x);L=limit(f,x,inf)L=exp(a)*b•例:求解单边极限问题symsx;limit((exp(x^3)-1)/(1-cos(sqrt(x-sin(x)))),x,0,'right')ans=124•在(-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')5•多变量函数的极限：–格式：L1=limit(limit(f,x,x0),y,y0)或L1=limit(limit(f,y,y0),x,x0)如果x0或y0不是确定的值，而是另一个变量的函数，如x-g(y)，则上述的极限求取顺序不能交换。6•例：求出二元函数极限值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)73.1.2函数导数的解析解•函数的导数和高阶导数–格式：y=diff(fun,x)%求导数y=diff(fun,x,n)%求n阶导数•例：一阶导数：symsx;f=sin(x)/(x^2+4*x+3);f1=diff(f);pretty(f1)8cos(x)sin(x)(2x+4)-----------------------------------222x+4x+3(x+4x+3)原函数及一阶导数图：x1=0:.01:5;y=subs(f,x,x1);y1=subs(f1,x,x1);plot(x1,y,x1,y1,‘:’)更高阶导数：tic,diff(f,x,100);tocelapsed_time=4.68609•原函数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)10•多元函数的偏导：–格式：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)11zy=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])12contour(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)%绘制引力线13•例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))14•多元函数的Jacobi矩阵:–格式：J=jacobian(Y,X)其中，X是自变量构成的向量，Y是由各个函数构成的向量。15•例：试推导其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]16•隐函数的偏导数：–格式：F=-diff(f,xj)/diff(f,xi)17•例：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)183.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+119•对原函数求4阶导数，再对结果进行4次积分y4=diff(y,4);y0=int(int(int(int(y4))));pretty(simple(y0))sin(x)------------2x+4x+320•例:说明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^421•定积分与无穷积分计算:–格式：I=int(f,x,a,b)–格式：I=int(f,x,a,inf)22•例：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23•多重积分问题的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)24f2=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)???f2=sin(x^2*y)/exp(y*x^2+z^2)simple(f1-f2)ans=0顺序的改变使化简结果不同于原函数，但其误差为0，表明二者实际完全一致。这是由于积分顺序不同，得不出实际的最简形式。25•例：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(n,z)为指数积分，无解析解，但可求其数值解：vpa(ans,60)ans=3.1080794020854127228346146476713852101914230631702186348358826Mathematica计算多重积分问题（参考）计算多重积分22200()axxydxdy273.2函数的级数展开与级数求和问题求解•3.2.1Taylor幂级数展开•3.2.2Fourier级数展开•3.2.3级数求和的计算283.2.1Taylor幂级数展开3.2.1.1单变量函数的Taylor幂级数展开2930例:symsx;f=sin(x)/(x^2+4*x+3);y1=taylor(f,x,9);pretty(y1)2334408730675152733864591/3x-4/9x2+--x3-----x4+------x5-------x6+----------x7----------x8548197207290122472091854031taylor(f,x,9,2)ans=1/15*sin(2)+(1/15*cos(2)-8/225*sin(2))*(x-2)+(-127/6750*sin(2)-8/225*cos(2))*(x-2)^2+(23/6750*cos(2)+628/50625*sin(2))*(x-2)^3+(-15697/6075000*sin(2)+28/50625*cos(2))*(x-2)^4+(203/6075000*cos(2)+6277/11390625*sin(2))*(x-2)^5+(-585671/2733750000*sin(2)-623/11390625*cos(2))*(x-2)^6+(262453/19136250000*cos(2)+397361/51