中科院-matlab在科学计算中的应用优秀PPT.ppt

《中科院-matlab在科学计算中的应用优秀PPT.ppt》由会员分享，可在线阅读，更多相关《中科院-matlab在科学计算中的应用优秀PPT.ppt（109页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第三章 微积分问题的计算机求解微积分问题的解析解微积分问题的解析解函数的级数绽开与级数求和问题求解函数的级数绽开与级数求和问题求解数值微分数值微分数值积分问题数值积分问题曲线积分与曲面积分的计算曲线积分与曲面积分的计算13.1 微积分问题的解析解 3.1.1 极限问题的解析解单变量函数的极限格式1：L=limit(fun,x,x0)格式2：L=limit(fun,x,x0,left 或 right)2例:试求解极限问题 syms x a b;f=x*(1+a/x)x*sin(b/x);L=limit(f,x,inf)L=exp(a)*b例:求解单边极限问题 syms x;limit(exp(x

2、3)-1)/(1-cos(sqrt(x-sin(x),x,0,right)ans=123在(-0.1,0.1)区间绘制出函数曲线：x=-0.1:0.001:0.1;y=(exp(x.3)-1)./(1-cos(sqrt(x-sin(x);Warning:Divide by zero.(Type warning off MATLAB:divideByZero to suppress this warning.)?plot(x,y,-,0,12,o)4多变量函数的极限：格式：L1=limit(limit(f,x,x0),y,y0)或 L1=limit(limit(f,y,y0),x,x0)假如x0

3、 或y0不是确定的值，而是另一个变量的函数，如x-g(y)，则上述的极限求取依次不能交换。5例：求出二元函数极限值 syms x y a;f=exp(-1/(y2+x2)*sin(x)2/x2*(1+1/y2)(x+a2*y2);L=limit(limit(f,x,1/sqrt(y),y,inf)L=exp(a2)63.1.2 函数导数的解析解函数的导数和高阶导数格式：y=diff(fun,x)%求导数 y=diff(fun,x,n)%求n阶导数例：一阶导数：syms x;f=sin(x)/(x2+4*x+3);f1=diff(f);pretty(f1)7 cos(x)sin(x)(2 x+4

4、)-2 2 2 x +4 x+3 (x +4 x+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.68608原函数4阶导数 f4=diff(f,x,4);pretty(f4)2 sin(x)cos(x)(2 x+4)sin(x)(2 x+4)-+4-12-2 2 2 2 3 x +4 x+3 (x +4 x+3)(x +4 x+3)3 sin(x)cos(x)(2 x+4)cos(x)(2 x+4)+12-24-+

5、48-2 2 2 4 2 3 (x +4 x+3)(x +4 x+3)(x +4 x+3)4 2 sin(x)(2 x+4)sin(x)(2 x+4)sin(x)+24-72-+24-2 5 2 4 2 3 (x +4 x+3)(x +4 x+3)(x +4 x+3)9多元函数的偏导：格式：f=diff(diff(f,x,m),y,n)或 f=diff(diff(f,y,n),x,m)例：求其偏导数并用图表示。syms x y z=(x2-2*x)*exp(-x2-y2-x*y);zx=simple(diff(z,x)zx=-exp(-x2-y2-x*y)*(-2*x+2+2*x3+x2*y-

6、4*x2-2*x*y)10 zy=diff(z,y)zy=(x2-2*x)*(-2*y-x)*exp(-x2-y2-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(-3 3-2 2-0.7 1.5)11 contour(x,y,z,30),hold on%绘制等值线 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.*

7、y);%偏导的数值解 quiver(x,y,zx,zy)%绘制引力线12例 syms x y z;f=sin(x2*y)*exp(-x2*y-z2);df=diff(diff(diff(f,x,2),y),z);df=simple(df);pretty(df)2 2 2 2 2 -4 z exp(-x y-z)(cos(x y)-10 cos(x y)y x +4 2 4 2 2 4 2 2sin(x y)x y+4 cos(x y)x y -sin(x y)13多元函数的Jacobi矩阵:格式：J=jacobian(Y,X)其中，X是自变量构成的向量，Y是由各个函数构成的向量。14例：试推导

8、其 Jacobi 矩阵 syms r theta phi;x=r*sin(theta)*cos(phi);y=r*sin(theta)*sin(phi);z=r*cos(theta);J=jacobian(x;y;z,r theta phi)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 15隐函数的偏导数：格式：F=-diff(f,

9、xj)/diff(f,xi)16例：syms x y;f=(x2-2*x)*exp(-x2-y2-x*y);pretty(-simple(diff(f,x)/diff(f,y)3 2 2 -2 x+2+2 x +x y-4 x -2 x y -x(x-2)(2 y+x)173.1.3 积分问题的解析解不定积分的推导：格式：F=int(fun,x)例：用diff()函数求其一阶导数，再积分，检验是否可以得出一样的结果。syms x;y=sin(x)/(x2+4*x+3);y1=diff(y);y0=int(y1);pretty(y0)%对导数积分 sin(x)sin(x)-1/2-+1/2-x+

10、3 x+118对原函数求对原函数求4 4 阶导数，再对结果进行阶导数，再对结果进行4 4次积分次积分 y4=diff(y,4);y0=int(int(int(int(y4);pretty(simple(y0)sin(x)-2 x +4 x+319例:说明明 syms a x;f=simple(int(x3*cos(a*x)2,x)f=1/16*(4*a3*x3*sin(2*a*x)+2*a4*x4+6*a2*x2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a4 f1=x4/8+(x3/(4*a)-3*x/(8*a3)*sin(2*a*x)+.(3*

11、x2/(8*a2)-3/(16*a4)*cos(2*a*x);simple(f-f1)%求两个结果的差ans=-3/16/a420定积分与无穷积分计算:格式：I=int(f,x,a,b)格式：I=int(f,x,a,inf)21例：syms x;I1=int(exp(-x2/2),x,0,1.5)无解无解I1=1/2*erf(3/4*2(1/2)*2(1/2)*pi(1/2)vpa(I1,70)ans=I2=int(exp(-x2/2),x,0,inf)I2=1/2*2(1/2)*pi(1/2)22多重积分问题的MATLAB求解例：syms x y z;f0=-4*z*exp(-x2*y-z2

12、)*(cos(x2*y)-syms x y z;f0=-4*z*exp(-x2*y-z2)*(cos(x2*y)-10*cos(x2*y)*y*x2+.10*cos(x2*y)*y*x2+.4*sin(x2*y)*x4*y2+4*cos(x2*y)*x4*y2-sin(x2*y);4*sin(x2*y)*x4*y2+4*cos(x2*y)*x4*y2-sin(x2*y);f1=int(f0,f1=int(f0,z z);f1=int(f1,);f1=int(f1,y y);f1=int(f1,);f1=int(f1,x x););f1=simple(int(f1,f1=simple(int(f

13、1,x x)f1=f1=exp(-x2*y-z2)*sin(x2*y)exp(-x2*y-z2)*sin(x2*y)23 f2=int(f0,z);f2=int(f2,x);f2=int(f2,x);f2=simple(int(f2,y)f2=2*exp(-x2*y-z2)*tan(1/2*x2*y)/(1+tan(1/2*x2*y)2)?f2=sin(x2*y)/exp(y*x2+z2)simple(f1-f2)ans=0 依次的变更使化简结果不同于原函数，但其误差为0，表明二者实际完全一样。这是由于积分依次不同，得不出实际的最简形式。24例：syms x y z int(int(int(4

14、*x*z*exp(-x2*y-z2),x,0,1),y,0,pi),z,0,pi)ans=(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2)*pi2*hypergeom(1,2,-pi2)?Ei(n,z)为指数积分，无解析解，但可求其数值解：vpa(ans,60)ans=25263.2 函数的级数绽开与 级数求和问题求解3.2.1 Taylor 幂级数绽开3.2.2 Fourier 级数绽开3.2.3 级数求和的计算273.2.1 Taylor 幂级数绽开 3.2.1.1 单变量函数的 Taylor 幂级数绽开2829例:syms x;f=sin(x)/(x2+4*

15、x+3);y1=taylor(f,x,9);pretty(y1)23 34 4087 3067 515273 386459 1/3 x-4/9 x2 +-x3 -x4 +-x5 -x6 +-x7 -x8 54 81 9720 7290 1224720 91854030 taylor(f,x,9,2)ans=syms a;taylor(f,x,5,a)%结果较冗长，显示从略ans=sin(a)/(a2+3+4*a)+(cos(a)-sin(a)/(a2+3+4*a)*(4+2*a)/(a2+3+4*a)*(x-a)+(-sin(a)/(a2+3+4*a)-1/2*sin(a)-(cos(a)*a

16、2+3*cos(a)+4*cos(a)*a-4*sin(a)-2*sin(a)*a)/(a2+3+4*a)2*(4+2*a)/(a2+3+4*a)*(x-a)2+31例:对y=sinx进行Taylor幂级数绽开,并视察不同阶次的近似效果。x0=-2*pi:0.01:2*pi;y0=sin(x0);syms x;y=sin(x);plot(x0,y0,r-.),axis(-2*pi,2*pi,-1.5,1.5);hold on for n=8:2:16 p=taylor(y,x,n),y1=subs(p,x,x0);line(x0,y1)endp=x-1/6*x3+1/120*x5-1/5040

17、*x7p=x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9p=x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11p=x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11+1/6227020800*x13 32p=333.2.1.2 多变量函数的Taylor 幂级数绽开多变量函数 在的Taylor幂级数的绽开34例：？syms x y;f=(x2-2*x)*exp(-x2-y2-x*y);F=maple(mtaylor,f,x,y,8)F=mtaylor

18、(x2-2*x)*exp(-x2-y2-x*y),x,y,8)35 maple(readlib(mtaylor);读库，把函数调入内存 F=maple(mtaylor,f,x,y,8)F=-2*x+x2+2*x3-x4-x5+1/2*x6+1/3*x7+2*y*x2+2*y2*x-y*x3-y2*x2-2*y*x4-3*y2*x3-2*y3*x2-y4*x+y*x5+3/2*y2*x4+y3*x3+1/2*y4*x2+y*x6+2*y2*x5+7/3*y3*x4+2*y4*x3+y5*x2+1/3*y6*x syms a;F=maple(mtaylor,f,x=1,y=a,3);F=maple

19、(mtaylor,f,x=a,3)F=(a2-2*a)*exp(-a2-y2-a*y)+(a2-2*a)*exp(-a2-y2-a*y)*(-2*a-y)+(2*a-2)*exp(-a2-y2-a*y)*(x-a)+(a2-2*a)*exp(-a2-y2-a*y)*(-1+2*a2+2*a*y+1/2*y2)+exp(-a2-y2-a*y)+(2*a-2)*exp(-a2-y2-a*y)*(-2*a-y)*(x-a)2363.2.2 Fourier 级数绽开37function A,B,F=fseries(f,x,n,a,b)if nargin=3,a=-pi;b=pi;endL=(b-a)/

20、2;if a+b,f=subs(f,x,x+L+a);end变量区域互换A=int(f,x,-L,L)/L;B=;F=A/2;%计算a0for i=1:n an=int(f*cos(i*pi*x/L),x,-L,L)/L;bn=int(f*sin(i*pi*x/L),x,-L,L)/L;A=A,an;B=B,bn;F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L);endif a+b,F=subs(F,x,x-L-a);end 换回变量区域38例:syms x;f=x*(x-pi)*(x-2*pi);A,B,F=fseries(f,x,6,0,2*pi)A=0,0,0,

21、0,0,0,0 B=-12,3/2,-4/9,3/16,-12/125,1/18 F=12*sin(x)+3/2*sin(2*x)+4/9*sin(3*x)+3/16*sin(4*x)+12/125*sin(5*x)+1/18*sin(6*x)39例:syms x;f=abs(x)/x;%定义方波信号 xx=-pi:pi/200:pi;xx=xx(xx=0);xx=sort(xx,-eps,eps);%剔除零点 yy=subs(f,x,xx);plot(xx,yy,r-.),hold on%绘制出理论值并保持坐标系 for n=2:20 a,b,f1=fseries(f,x,n),y1=sub

22、s(f1,x,xx);plot(xx,y1)end40a=0,0,0b=4/pi,0f1=4/pi*sin(x)a=0,0,0,0 b=4/pi,0,4/3/pif1=4/pi*sin(x)+4/3/pi*sin(3*x)413.2.3 级数求和的计算是在符号工具箱中供应的42例:计算 format long;sum(2.0:63)%数值计算ans=1.844674407370955e+019 sum(sym(2).0:200)%或 syms k;symsum(2k,0,200)把2定义为符号量可使计算更精确更精确ans=syms k;symsum(2k,0,200)ans=43例:试求解无穷

23、级数的和 syms n;s=symsum(1/(3*n-2)*(3*n+1),n,1,inf)%接受符号运算工具箱s=1/3 m=1:10000000;s1=sum(1./(3*m-2).*(3*m+1);%数值计算方法,双精度有效位16,“大数吃小数”,无法精确 format long;s1%以长型方式显示得出的结果s1=0.3333333222216544例:求解 syms n x s1=symsum(2/(2*n+1)*(2*x+1)(2*n+1),n,0,inf);simple(s1)%对结果进行化简，MATLAB 6.5 及以前版本因本身 bug 化简很麻烦ans=log(2*x+1

24、)2)(1/2)+1)/(2*x+1)2)(1/2)-1)%实际应为log(x+1)/x)45例:求 syms m n;limit(symsum(1/m,m,1,n)-log(n),n,inf)ans=eulergamma vpa(ans,70)%显示 70 位有效数字ans=46符号函数计算器 单变量符号函数计算器 Taylor 靠近计算器 47单变量符号函数计算器（1/3）在叮嘱窗口中执行 funtool 即可调出单变量符号函数计算器。单变量符号函数计算器用于对单变量函数进行操作，可以对符号函数进行化简、求导、绘制图形等。该工具的界面如图所示。函数 f 的图形窗口 函数 g 的图形窗口 限

25、制窗口 48单变量符号函数计算器（2/3）1输入框的功能如图：函数函数 f 的编辑框的编辑框函数函数 g 的编辑框的编辑框显示绘制显示绘制 f 和和 g 的图像的的图像的 x 区间区间用于修改用于修改 f 的常数因子的常数因子0函数 f 自身的操作函数 f 与常数 a 的操作函数 f 与函数 g 的操作系统操作49单变量符号函数计算器（3/3）单变量符号函数计算器应用示例在 f 函数输入栏中输入 cos(x3)cos(x3)(1+x2)在 g 函数输入栏中输入(1+x2)点击50Taylor 靠近计算器 Taylor 靠近计算器用于实现函数的 taylor 靠近。在叮嘱窗口中输入 taylor

26、tool，调出Taylor 靠近计算器，界面及功能如图。输入待靠近的函数 输入拟合函数的阶数 级数的绽开点，默认为 0 输入拟合区间 51MAPLE 函数的调用 maple 函数的运用 mfun 函数的运用 52maple 函数的运用maple 是符号工具箱中的一个通用叮嘱，运用它可以实现对 MAPLE 中大部分函数的调用。其运用格式为：r=maple(statement)，其中 statement 为符合 MAPLE 语法的可执行语句的字符串，该叮嘱将 statement 传递给 MAPLE，该叮嘱的输出结果也符合 MAPLE 的语法；r=maple(function,arg1,arg2,.

27、)，该函数调用引号中的函数，并接受指定的参数，相当于 MAPLE 语句 function(arg1,arg2,.)；r,status=maple(.)，返回函数的运行状态，假如函数运行成功，则 status 为 0，r 为运行结果；假如函数运行失败，则 status 为一个正数，r 为相应的错误信息；maple(traceon)或者 maple trace on，输出 MAPLE 函数运行中的全部子表达式和运行结果；maple(traceoff)或 maple trace off，不显示中间过程。53mfun 函数的运用 mfun 函数用于对 maple 函数进行数字评估。该函数的调用格式为：

28、Y=mfun(function,par1,par2,par3,par4)。该语句对指定的数学函数进行评估。其中 function 指定待评估的函数，par1、par2 等为 function 的参数，为待评估的数值，其维数有 function 函数的参数类型确定。在该语句中最多可以设置四个参数，最终一个参数可以为矩阵。用户可以通过 help mfunlist 查看 MATLAB 中 mfun 可以调用的函数列表，另外，可以通过 mhelp function 查看指定函数的相关信息。543.3 数值微分 x-h x x+hBCAT f(x)55 3.3.1 数值微分算法向前差商公式：向后差商公式

29、56两种中心公式：5758593.3.2 中心差分方法及其 MATLAB 实现 function dy,dx=diff_ctr(y,Dt,n)yx1=y 0 0 0 0 0;yx2=0 y 0 0 0 0;yx3=0 0 y 0 0 0;yx4=0 0 0 y 0 0;yx5=0 0 0 0 y 0;yx6=0 0 0 0 0 y;switch n case 1 dy=(-diff(yx1)+7*diff(yx2)+7*diff(yx3)-diff(yx4)/(12*Dt);L0=3;case 2 dy=(-diff(yx1)+15*diff(yx2)-15*diff(yx3)+diff(yx

30、4)/(12*Dt2);L0=3;数值计算diff(X)表示数组X相邻两数的差60 case 3 dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+.7*diff(yx5)-diff(yx6)/(8*Dt3);L0=5;case 4 dy=(-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*diff(yx4)-11*diff(yx5)+diff(yx6)/(6*Dt4);L0=5;end dy=dy(L0+1:end-L0);dx=(1:length(dy)+L0-2-(n2)*Dt;调用格式：y为 等距实测数据

31、，dy为得出的导数向量，dx为相应的自变量向量，dy、dx的数据比y短。61例：求导数的解析解,再用数值微分求取原函数的14 阶导数，并和解析解比较精度。h=0.05;x=0:h:pi;syms x1;y=sin(x1)/(x12+4*x1+3);%求各阶导数的解析解与比照数据 yy1=diff(y);f1=subs(yy1,x1,x);yy2=diff(yy1);f2=subs(yy2,x1,x);yy3=diff(yy2);f3=subs(yy3,x1,x);yy4=diff(yy3);f4=subs(yy4,x1,x);62 y=sin(x)./(x.2+4*x+3);%生成已知数据点

32、y1,dx1=diff_ctr(y,h,1);subplot(221),plot(x,f1,dx1,y1,:);y2,dx2=diff_ctr(y,h,2);subplot(222),plot(x,f2,dx2,y2,:)y3,dx3=diff_ctr(y,h,3);subplot(223),plot(x,f3,dx3,y3,:);y4,dx4=diff_ctr(y,h,4);subplot(224),plot(x,f4,dx4,y4,:)求最大相对误差：norm(y4-f4(4:60)./f4(4:60)ans=3.5025e-004633.3.3 用插值、拟合多项式的求导数基本思想：当已知

33、函数在一些离散点上的函数值时，该函数可用插值或拟合多项式来近似，然后对多项式进行微分求得导数。选取x=0旁边的少量点进行多项式拟合或插值g(x)在x=0处的k阶导数为64通过坐标变换用上述方法计算随意x点处的导数值令将g(x)写成z的表达式导数为可干脆用 拟合节点 得到系数 d=polyfit(x-a,y,length(xd)-1)65例：数据集合如下：xd:0 0.2000 0.4000 0.6000 0.8000 1.000 yd:0.3927 0.5672 0.6982 0.7941 0.8614 0.9053计算x=a=0.3处的各阶导数。xd=0 0.2000 0.4000 0.60

34、00 0.8000 1.000;yd=0.3927 0.5672 0.6982 0.7941 0.8614 0.9053;a=0.3;L=length(xd);d=polyfit(xd-a,yd,L-1);fact=1;for k=1:L-1;fact=factorial(k),fact;end deriv=d.*factderiv=1.8750 -1.3750 1.0406 -0.9710 0.6533 0.637666建立用拟合（插值）多项式计算各阶导数的poly_drv.mfunction der=poly_drv(xd,yd,a)m=length(xd)-1;d=polyfit(xd-

35、a,yd,m);c=d(m:-1:1);去掉常数项fact(1)=1;for i=2:m;fact(i)=i*fact(i-1);endder=c.*fact;例：xd=0 0.2000 0.4000 0.6000 0.8000 1.000;yd=0.3927 0.5672 0.6982 0.7941 0.8614 0.9053;a=0.3;der=poly_drv(xd,yd,a)der=0.6533 -0.9710 1.0406 -1.3750 1.8750673.3.4 二元函数的梯度计算格式：若z矩阵是建立在等间距的形式生成的网格基础上，则实际梯度为68例：计算梯度，绘制引力线图：x,

36、y=meshgrid(-3:.2:3,-2:.2:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);fx,fy=gradient(z);fx=fx/0.2;fy=fy/0.2;contour(x,y,z,30);hold on;quiver(x,y,fx,fy)%绘制等高线与引力线图69绘制误差曲面: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);surf(x,y,abs(fx-zx);axis(-3 3-2 2 0,0.

37、08)figure;surf(x,y,abs(fy-zy);axis(-3 3-2 2 0,0.11)建立一个新图形窗口70为削减误差，对网格加密一倍：x,y=meshgrid(-3:.1:3,-2:.1:2);z=(x.2-2*x).*exp(-x.2-y.2-x.*y);fx,fy=gradient(z);fx=fx/0.1;fy=fy/0.1;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);surf(x,y,abs(fx-zx);axi

38、s(-3 3-2 2 0,0.02)figure;surf(x,y,abs(fy-zy);axis(-3 3-2 2 0,0.06)713.4 数值积分问题 4.3.1 由给定数据进行梯形求积72Sum(2*y(1:end-1,:)+diff(y).*diff(x)/273格式：S=trapz(x,y)例：x1=0:pi/30:pi;y=sin(x1)cos(x1)sin(x1/2);x=x1 x1 x1;S=sum(2*y(1:end-1,:)+diff(y).*diff(x)/2S=1.9982 0.0000 1.9995 S1=trapz(x1,y)%得出和上述完全一样的结果S1=1.9

39、982 0.0000 1.999574例：画图 x=0:0.01:3*pi/2,3*pi/2;%这样赋值能确保 3*pi/2点被包含在内 y=cos(15*x);plot(x,y)%求取理论值 syms x,A=int(cos(15*x),0,3*pi/2)A=1/1575随着步距h的减小，计算精度渐渐增加：h0=0.1,0.01,0.001,0.0001,0.00001,0.000001;v=;for h=h0,x=0:h:3*pi/2,3*pi/2;y=cos(15*x);I=trapz(x,y);v=v;h,I,1/15-I;end format long;vv=0.0000100000

40、0000 0.06666666654167 0.00000000012500 0.00000100000000 0.06666666666542 0.00000000000125 763.4.2 单变量数值积分问题求解 梯形公式格式：(变步长）(Fun:函数的字符串变量）y=quad(Fun,a,b)y=quadl(Fun,a,b)%求定积分 y=quad(Fun,a,b,)y=quadl(Fun,a,b,)%限定精度的定积分求解，默认精度为106。后面函数算法更精确，精度更高。77例：第三种：匿名函数(MATLAB 7.0)其次种：inline 函数第一种，一般函数方法78函数定义被积函数：

41、y=quad(c3ffun,0,1.5)y=0.9661用 inline 函数定义被积函数：f=inline(2/sqrt(pi)*exp(-x.2),x);y=quad(f,0,1.5)y=0.9661运用符号工具箱：syms x,y0=vpa(int(2/sqrt(pi)*exp(-x2),0,1.5),60)y0=y=quad(f,0,1.5,1e-20)%设置高精度，但该方法失效79提高求解精度：y=quadl(f,0,1.5,1e-20)y=0.9661 abs(y-y0)ans=.6402522848913892e-16 format long 16位精度 y=quadl(f,0,

42、1.5,1e-20)y=80例：求解绘制函数:x=0:0.01:2,2+eps:0.01:4,4;y=exp(x.2).*(x2);y(end)=0;x=eps,x;y=0,y;fill(x,y,g)为削减视觉上的误差，对端点与间断点（有跳动）进行处理。81调用quad():f=inline(exp(x.2).*(x2)./(4-sin(16*pi*x),x);I1=quad(f,0,4)I1=57.76435412500863调用quadl():I2=quadl(f,0,4)I2=57.76445016946768 syms x;I=vpa(int(exp(x2),0,2)+int(80/(

43、4-sin(16*pi*x),2,4)I=823.4.3 Gauss求积公式为使求积公式得到较高的代数精度对求积区间a,b，通过变换有83以n=2的高斯公式为例：function g=gauss2(fun,a,b)h=(b-a)/2;c=(a+b)/2;x=h*(-0.7745967)+c,c,h*0.7745967+c;g=h*(0.55555556*(gaussf(x(1)+gaussf(x(3)+0.88888889*gaussf(x(2);function y=gaussf(x)y=cos(x);gauss2(gaussf,0,1)ans=0.8415843.4.4 双重积分问题的数值

44、解矩形区域上的二重积分的数值计算格式：矩形区域的双重积分：y=dblquad(Fun,xm,xM,ym,yM)限定精度的双重积分：y=dblquad(Fun,xm,xM,ym,yM，)85例：求解 f=inline(exp(-x.2/2).*sin(x.2+y),x,y);y=dblquad(f,-2,2,-1,1)y=1.5744931897449486随意区域上二元函数的数值积分(调用工具箱NIT），该函数指定依次先x后y.87例 fh=inline(sqrt(1-x.2/2),x);%内积分上限 fl=inline(-sqrt(1-x.2/2),x);%内积分下限 f=inline(ex

45、p(-x.2/2).*sin(x.2+y),y,x);%交换依次的被积函数 y=quad2dggen(f,fl,fh,-1/2,1,eps)y=88解析解方法：syms x y i1=int(exp(-x2/2)*sin(x2+y),y,-sqrt(1-x2/2),sqrt(1-x2/2);int(i1,x,-1/2,1)Warning:Explicit integral could not be found.In D:MATLAB6p5toolboxsymbolicsymint.m at line 58 ans=int(2*exp(-1/2*x2)*sin(x2)*sin(1/2*(4-2*

46、x2)(1/2),x=-1/2.1)vpa(ans)ans=89例：计算单位圆域上的积分：先把二重积分转化：syms x y i1=int(exp(-x2/2)*sin(x2+y),x,-sqrt(1-y.2),sqrt(1-y.2);Warning:Explicit integral could not be found.In D:MATLAB6p5toolboxsymbolicsymint.m at line 5890对x是不行积的，故调用解析解方法不会得出结果，而数值解求解不受此影响。fh=inline(sqrt(1-y.2),y);%内积分上限 fl=inline(-sqrt(1-y.

47、2),y);%内积分下限 f=inline(exp(-x.2/2).*sin(x.2+y),x,y);%交换依次的被积函数 I=quad2dggen(f,fl,fh,-1,1,eps)Integral did not converge-singularity likelyI=0.53686038269795913.4.5 三重定积分的数值求解格式:I=triplequad(Fun,xm,xM,ym,yM,zm,zM，,quadl)其中quadl为具体求解一元积分的数值函数,也可选用quad或自编积分函数,但调用格式要与quadl一样。92例：triplequad(inline(4*x.*z.*

48、exp(-x.*x.*y-z.*z),x,y,z),0,1,0,pi,0,pi,1e-7,quadl)ans=1.7328933.5 曲线积分与曲面积分的计算3.5.1 曲线积分及MATLAB求解第一类曲线积分 起源于对不匀整分布的空间曲线总质量的求取.设空间曲线L的密度函数为f(x,y,z),则其总质量 其中s为曲线上某点的弧长,又称这类曲线积分为对弧长的曲线积分.94数学表示 若弧长表示为95例：syms t;syms a positive;x=a*cos(t);y=a*sin(t);z=a*t;I=int(z2/(x2+y2)*sqrt(diff(x,t)2+diff(y,t)2+dif

49、f(z,t)2),t,0,2*pi)I=8/3*pi3*a*2(1/2)pretty(I)3 1/2 8/3 pi a 296例：x=0:.001:1.2;y1=x;y2=x.2;plot(x,y1,x,y2)绘出两条曲线 syms x;y1=x;y2=x2;I1=int(x2+y22)*sqrt(1+diff(y2,x)2),x,0,1);I2=int(x2+y12)*sqrt(1+diff(y1,x)2),x,1,0);I=I2+I1I=-2/3*2(1/2)+349/768*5(1/2)+7/512*log(-2+5(1/2)973.5.1.2 其次类曲线积分又称对坐标的曲线积分，起源于

50、变力沿曲线 移动时作功的探讨曲线 亦为向量,若曲线可以由参数方程表示则两个向量的点乘可由这两个向量干脆得出.98例：求曲线积分 syms t;syms a positive;x=a*cos(t);y=a*sin(t);F=(x+y)/(x2+y2),-(x-y)/(x2+y2);ds=diff(x,t);diff(y,t);I=int(F*ds,t,2*pi,0)%正向圆周 I=2*pi99例：syms x;y=x2;F=x2-2*x*y,y2-2*x*y;ds=1;diff(y,x);I=int(F*ds,x,-1,1)I=-14/151003.5.2曲面积分与MATLAB语言求解3.5.2