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

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

1、1,第三章 微积分问题的计算机求解,微积分问题的解析解 函数的级数展开与级数求和问题求解 数值微分 数值积分问题 曲线积分与曲面积分的计算,2,3.1 微积分问题的解析解 3.1.1 极限问题的解析解,单变量函数的极限 格式1： L= limit( fun, x, x0) 格式2： L= limit( fun, x, x0, left 或 right),3,例: 试求解极限问题 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(x3)-1)/(1-cos(s

2、qrt(x-sin(x),x,0,right) ans = 12,4,在(-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),5,多变量函数的极限： 格式： L1=limit(limit(f,x,x0),y,y0) 或 L1=limit(limit(f,y,y0), x,x0

3、) 如果x0 或y0不是确定的值，而是另一个变量的函数，如x-g(y)，则上述的极限求取顺序不能交换。,6,例：求出二元函数极限值 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),7,3.1.2 函数导数的解析解,函数的导数和高阶导数 格式： y=diff(fun,x) %求导数 y= diff(fun,x,n) %求n阶导数 例： 一阶导数： syms x; f=sin(x)/(x2+4*x+3); f1=diff(f); pr

4、etty(f1),8,cos(x) sin(x) (2 x + 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); toc elapsed_time = 4.6860,9,原函数4阶导数 f4=diff(f,x,4); pretty(f4) 2 sin(x) cos(x) (2 x + 4) sin(x) (2 x + 4) - + 4 - - 12 -

5、 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 - + 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),10,多元函数的偏导：

6、 格式： 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-4*x2-2*x*y),11, 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

7、-x.*y); surf(x,y,z), axis(-3 3 -2 2 -0.7 1.5),12, 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.*y); % 偏导的数值解 quiver(x,y,zx,zy) % 绘制引力线,13,例 syms x y z; f=sin(x2*y)*exp(-x2*y-z2); df=diff(diff(diff(f,x,2),y),z); d

8、f=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 2 sin(x y) x y+ 4 cos(x y) x y - sin(x y),14,多元函数的Jacobi矩阵: 格式：J=jacobian(Y,X) 其中，X是自变量构成的向量，Y是由各个函数构成的向量。,15,例： 试推导其 Jacobi 矩阵 syms r theta phi; x=r*sin(theta)*cos(phi); y=r*sin(theta)*sin(phi); z=r*c

9、os(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 ,16,隐函数的偏导数： 格式：F=-diff(f,xj)/diff(f,xi),17,例： syms x y; f=(x2-2*x)*exp(-x2-y2-x*y); pretty

10、(-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),18,3.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 + 3 x + 1,19,对原函数求4

11、 阶导数，再对结果进行4次积分 y4=diff(y,4); y0=int(int(int(int(y4); pretty(simple(y0) sin(x) - 2 x + 4 x + 3,20,例:说明 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*x2/(8*a2)-3/(16*a4)*

12、cos(2*a*x); simple(f-f1) % 求两个结果的差 ans = -3/16/a4,21,定积分与无穷积分计算: 格式： I=int(f,x,a,b) 格式： I=int(f,x,a,inf),22,例： 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 = 1.085853317666016569702419076542265042534236293532156326729917229308528 I2=int(exp(-x2/2),x,0,i

13、nf) I2 = 1/2*2(1/2)*pi(1/2),23,多重积分问题的MATLAB求解 例： syms x y z; f0=-4*z*exp(-x2*y-z2)*(cos(x2*y)-10*cos(x2*y)*y*x2+. 4*sin(x2*y)*x4*y2+4*cos(x2*y)*x4*y2-sin(x2*y); f1=int(f0,z);f1=int(f1,y);f1=int(f1,x); f1=simple(int(f1,x) f1 = exp(-x2*y-z2)*sin(x2*y),24, f2=int(f0,z); f2=int(f2,x); f2=int(f2,x); f2=

14、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，表明二者实际完全一致。这是由于积分顺序不同，得不出实际的最简形式。,25,例： syms x y z int(int(int(4*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)*pi

15、2*hypergeom(1,2,-pi2) ? Ei(n,z)为指数积分，无解析解，但可求其数值解： vpa(ans,60) ans = 3.10807940208541272283461464767138521019142306317021863483588,26,27,3.2 函数的级数展开与 级数求和问题求解,3.2.1 Taylor 幂级数展开 3.2.2 Fourier 级数展开 3.2.3 级数求和的计算,28,3.2.1 Taylor 幂级数展开 3.2.1.1 单变量函数的 Taylor 幂级数展开,29,30,例: syms x; f=sin(x)/(x2+4*x+3); y

16、1=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 918540,31, taylor(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 +(-156

17、97/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/5125781250*sin(2)*(x-2)7 +(-875225059/34445250000000*sin(2)-131623/35880468750*cos(2)*(x-2)8 syms a; taylor(f,x,5,a) % 结

18、果较冗长，显示从略 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)*a2+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+,32,例:对y=sinx进行Taylor幂级数展开,并观察不同阶次的近似效果。 x0=-2*pi:0.01:2*pi; y0=sin(x0); syms x; y=sin(x); p

19、lot(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) end p = x-1/6*x3+1/120*x5-1/5040*x7 p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9 p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11 p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916

20、800*x11+1/6227020800*x13,33,p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11+1/6227020800*x13-1/1307674368000*x15,34,3.2.1.2 多变量函数的Taylor 幂级数展开,多变量函数 在 的Taylor幂级数的展开,35,例：？ syms x y; f=(x2-2*x)*exp(-x2-y2-x*y); F=maple(mtaylor,f,x,y,8) F = mtaylor(x2-2*x)*exp(-x2-y2-x*y),x, y,8),36, maple(

21、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(mtaylor,f,x=a,3) F = (a2-2*a)*exp(

22、-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)2,37,3.2.2 Fourier 级数展开,38,function A,B,F=fseries(f,x,n,a,b) if nargin=3, a=-pi; b=pi; end L=(b-a)/2; if a+b, f=subs(f,x,x

23、+L+a); end变量区域互换 A=int(f,x,-L,L)/L; B=; F=A/2; %计算a0 for 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); end if a+b, F=subs(F,x,x-L-a); end 换回变量区域,39,例: syms x; f=x*(x-pi)*(x-2*pi); A,B,F=fseries(f,x,6,0,2*pi) A = 0, 0

24、, 0, 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),40,例: 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:

25、20 a,b,f1=fseries(f,x,n), y1=subs(f1,x,xx); plot(xx,y1) end,41,a = 0, 0, 0 b = 4/pi, 0 f1 = 4/pi*sin(x) a = 0, 0, 0, 0 b = 4/pi, 0, 4/3/pi f1 = 4/pi*sin(x)+4/3/pi*sin(3*x) ,42,3.2.3 级数求和的计算,是在符号工具箱中提供的,43,例:计算 format long; sum(2.0:63) %数值计算 ans = 1.844674407370955e+019 sum(sym(2).0:200) % 或 syms k;

26、symsum(2k,0,200) 把2定义为符号量可使计算更精确 ans = 3213876088517980551083924184682325205044405987565585670602751 syms k; symsum(2k,0,200) ans = 3213876088517980551083924184682325205044405987565585670602751,44,例:试求解无穷级数的和 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

27、-2).*(3*m+1);%数值计算方法,双精度有效位16,“大数吃小数”,无法精确 format long; s1 % 以长型方式显示得出的结果 s1 = 0.33333332222165,45,例:求解 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)2)(1/2)+1)/(2*x+1)2)(1/2)-1) %实际应为log(x+1)/x),46,例:求 syms m n; limit(symsu

28、m(1/m,m,1,n)-log(n),n,inf) ans = eulergamma vpa(ans, 70) % 显示 70 位有效数字 ans = .5772156649015328606065120900824024310421593359399235988057672348848677,47,符号函数计算器,单变量符号函数计算器 Taylor 逼近计算器,48,单变量符号函数计算器（1/3）,在命令窗口中执行 funtool 即可调出单变量符号函数计算器。单变量符号函数计算器用于对单变量函数进行操作，可以对符号函数进行化简、求导、绘制图形等。该工具的界面如图所示。,函数 f 的图形窗

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

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

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

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

33、前差商公式： 向后差商公式,57,两种中心公式：,58,59,60,3.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)+

34、15*diff(yx2)- 15*diff(yx3) +diff(yx4)/(12*Dt2);L0=3; 数值计算diff(X)表示数组X相邻两数的差,61,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)

35、; dx=(1:length(dy)+L0-2-(n2)*Dt; 调用格式： y为 等距实测数据， dy为得出的导数向量， dx为相应的自变量向量，dy、dx的数据比y短 。,62,例： 求导数的解析解,再用数值微分求取原函数的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);

36、yy4=diff(yy3); f4=subs(yy4,x1,x);,63, y=sin(x)./(x.2+4*x+3); % 生成已知数据点 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,:)

37、 求最大相对误差： norm(y4- f4(4:60)./f4(4:60) ans = 3.5025e-004,64,3.3.3 用插值、拟合多项式的求导数,基本思想：当已知函数在一些离散点上的函数值时，该函数可用插值或拟合多项式来近似，然后对多项式进行微分求得导数。 选取x=0附近的少量点 进行多项式拟合或插值 g(x)在x=0处的k阶导数为,65,通过坐标变换用上述方法计算任意x点处的导数值 令 将g(x)写成z的表达式 导数为 可直接用 拟合节点 得到系数 d=polyfit(x-a,y,length(xd)-1),66,例：数据集合如下： xd: 0 0.2000 0.4000 0.6

38、000 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.6000 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.*fact deriv = 1.8750 -1.3750 1.0406

39、 -0.9710 0.6533 0.6376,67,建立用拟合（插值）多项式计算各阶导数的poly_drv.m function der=poly_drv(xd,yd,a) m=length(xd)-1; d=polyfit(xd-a,yd,m); c=d(m:-1:1); 去掉常数项 fact(1)=1;for i=2:m; fact(i)=i*fact(i-1);end der=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; d

40、er=poly_drv(xd,yd,a) der = 0.6533 -0.9710 1.0406 -1.3750 1.8750,68,3.3.4 二元函数的梯度计算,格式： 若z矩阵是建立在等间距的形式生成的网格基础上，则实际梯度为,69,例： 计算梯度，绘制引力线图： x,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) %绘制等高线与 引力线图,

41、70,绘制误差曲面: 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.08) figure; surf(x,y,abs(fy-zy); axis(-3 3 -2 2 0,0.11) 建立一个新图形窗口,71,为减少误差，对网格加密一倍： x,y=meshgrid(-3:.1:3,-2:.1:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y)

42、; 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); axis(-3 3 -2 2 0,0.02) figure; surf(x,y,abs(fy-zy); axis(-3 3 -2 2 0,0.06),72,3.4 数值积分问题 4.3.1 由给定数据进行梯形求积,73,Sum(2*y(1:end-1,:)+diff

43、(y).*diff(x)/2,74,格式： 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)/2 S = 1.9982 0.0000 1.9995 S1=trapz(x1,y) % 得出和上述完全一致的结果 S1 = 1.9982 0.0000 1.9995,75,例： 画图 x=0:0.01:3*pi/2, 3*pi/2; % 这样赋值能确保 3*pi/2点被包含在内 y=cos(15*x); plot(x,y) % 求取理

44、论值 syms x, A=int(cos (15*x),0,3*pi/2) A = 1/15,76,随着步距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; v v = 0.10000000000000 0.05389175150076 0.01277491516591 0.01000000000000 0.066541695465

45、84 0.00012497120083 0.00100000000000 0.06666541668004 0.00000124998663 0.00010000000000 0.06666665416667 0.00000001250000 0.00001000000000 0.06666666654167 0.00000000012500 0.00000100000000 0.06666666666542 0.00000000000125,77,3.4.2 单变量数值积分问题求解,梯形公式 格式：(变步长）(Fun:函数的字符串变量） y=quad(Fun,a,b) y=quadl(Fun

46、,a,b) % 求定积分 y=quad(Fun,a,b, ) y=quadl(Fun,a,b, ) %限定精度的定积分求解，默认精度为106。后面函数算法更精确，精度更高。,78,例：,第三种：匿名函数(MATLAB 7.0),第二种：inline 函数,第一种，一般函数方法,79,函数定义被积函数： 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(p

47、i)*exp(-x2),0,1.5),60) y0 = .966105146475310713936933729949905794996224943257461473285749 y=quad(f,0,1.5,1e-20) % 设置高精度，但该方法失效,80,提高求解精度： y=quadl(f,0,1.5,1e-20) y = 0.9661 abs(y-y0) ans = .6402522848913892e-16 format long 16位精度 y=quadl(f,0,1.5,1e-20) y = 0.96610514647531,81,例：求解 绘制函数: 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) 为减少视觉上的误 差，对端点与间断点 （有跳跃）进行处理。,82,调用quad( ): f=inline(exp(x.2).*(x2)./(4-sin(16*pi*x)