数值积分与matlab求解省名师优质课获奖课件市赛课一等奖课件.ppt

MATLAB求解数值积分第1页 MATLAB求解连续函数积分第2页引言我们知道,若函数f(x)在区间a,b上连续且其原函数为F(x),则可用Newton-Leibnitz公式求定积分值,Newton-Leibnitz公式不论在理论上还是在处理实际问题上都起了很大作用，但它并不能完全处理定积分计算问题，因为积分学包括实际问题极为广泛，而且极其复杂，在实际计算中经常碰到以下三种情况：1 Matlab求解连续函数积分第3页(1)被积函数f(x)并不一定能够找到用初等函数有限形式表示原函数F(x)，比如：Newton-Leibnitz公式就无能为力了(2)还有被积函数f(x)原函数能用初等函数表示，但表示式太复杂，比如函数并不复杂，但积分后其表示式却很复杂，积分后其原函数F(x)为：第4页(3)被积函数f(x)没有详细解析表示式,其函数关系由表格或图形表示。对于这些情况,要计算积分准确值都是十分困难。由此可见,经过原函数来计算积分有它不足,因而研究一个新积分方法来处理Newton-Leibniz公式所不能或极难处理积分问题,这时需要用数值解法来建立积分近似计算方法。将积分区间细分,在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分思想，用代数插值多项式去代替被积函数发f(x)进行积分是本章讨论数值积分主要内容。第5页1.1 定积分Matlab符号计算例1 由 y=sin x，y=cos x,x=-1/2,x=3/2所围成平面区域D.求平面区域D 面积S.解 输入作函数图形程序 x=-1:0.001:2;F1=sin(x);F2=cos(x);plot(x,F1,b-,x,F2,g-),axis(-1,pi/4+1,-1.3,1.3),xlabel(x),ylabel(y),title(y=sinx,y=cosx 和x=-0.5及x=1.5所围成平面区域图形)运行后屏幕显示图形.求平面区域D 面积S.输入计算面积S 程序 syms xf1=cos(x)-sin(x);f2=-f1;S1=int(f1,x,-0.5,pi/4);S2=int(f2,x,pi/4,1.5);S=S1+S2,Sj=double(S)运行后屏幕显示计算面积值 S 及其近似值Sj 以下S=2*2(1/2)+sin(1/2)-cos(1/2)-sin(3/2)-cos(3/2)Sj=1.36203791318826坐标调整第6页1.2 变限积分Matlab符号计算例2 已知 ,求F(x)解：输入程序：syms x tF1=int(exp(t)*sin(2+sqrt(t3),x,0);F2=int(exp(t)*sin(2+sqrt(t3),0,x2);Fi=F1+F2;dF=diff(Fi)运行后屏幕显示计算变限积分F(x)导数.第7页 建立数值积分公式路径比较多,其中最惯用有两种：（1）由积分中值定理可知，对于连续函数f(x)，在积分区间a,b内存在一点，使得即所求曲边梯形面积恰好等于底为(b-a)，高为矩形面积。不过点详细位置普通是未知,因而值也是未知,称为f(x)在区间a,b上平均高度。那么只要对平均高度提供一个算法，对应地就取得一个数值求积方法1.3 数值求积方法第8页 梯形公式 矩形公式按照这种思想，可结构出一些求积分值近似公式。比如 分别取 和则分别得到中矩形公式和梯形公式。第9页 Simpson公式矩形公式把a,b中点处函数值作为平均高度f()近似值而取得一个数值积分方法。在这三个公式中,梯形公式把f(a),f(b)加权平均值作为平均高度f()近似值而取得一个数值积分方法。第10页Simpson公式是以函数f(x)在a,b,(a+b)/2这三点函数值f(a),f(b),加权平均值似值而取得一个数值积分方法。作为平均高度f()近（2）先用某个简单函数近似迫近f(x),用代替原被积函数f(x)，即以此结构数值算法。从数值计算角度考虑,函数应对f(x)有充分迫近程度,而且轻易计算其积分。因为多项式能很好地迫近连续函数,且又轻易计算积分,所以将选取为插值多项式,这么f(x)积分就能够用其插值多项式积分来近似代替第11页(一)用函数 trapz 计算定积分调用格式一：Z=trapz(Y)调用格式二：Z=trapz(X,Y)调用格式三：Z=trapz(X,Y,DIM)或 trapz(Y,DIM)(二)用函数 cumtrapz 计算定积分调用格式一：Z=cumtrapz(Y)调用格式二：Z=cumtrapz(X,Y)调用格式三：Z=cumtrapz(X,Y,DIM)或 cumtrapz(Y,DIM)1.3.1 梯形公式Matlab程序第12页梯形数值积分MATLAB 主程序function T=rctrap(fun,a,b,m)n=1;h=b-a;T=zeros(1,m+1);x=a;T(1)=h*(feval(fun,a)+feval(fun,b)/2;for i=1:mh=h/2;n=2*n;s=0;for k=1:n/2x=a+h*(2*k-1);s=s+feval(fun,x);endT(i+1)=T(i)/2+h*s;endT=T(1:m);第13页例3 用 MATLAB 函数trapz 和 cumtrapz 分别计算 准确到10-4，并与矩形公式比较.解：将0,/2分成20 等份，步长为/40，输入程序以下（注意trapz(y)是单位步长,trapz(y)*h=trapz(x,y)）：h=pi/40;x=0:h:pi/2;y=exp(-x).*sin(x);z1=sum(y(1:20)*h,z2=sum(y(2:21)*h,z=(z1+z2)/2z3=trapz(y)*h,z3h=trapz(x,y),z3c=cumtrapz(y)*h,运行后屏幕显示用矩形公式（9.3），（9.4）计算结果z1、z2 和二者平均数z、函数trapz和cumtrapz 分别计算结果z3、z3c.第14页（一）函数sum 调用格式调用格式一：sum(X)调用格式二：sum(X,DIM)（二）函数cumsum 调用格式调用格式一：cumsum(X)调用格式二：cumsum(X,DIM)1.3.2 矩形公式Matlab程序第15页例4 用 MATLAB 函数sum 和 cumsum 及矩形公式计算 ，并与准确值比较.解：将0,/2分成20 等份，步长为/40，输入程序以下（注意sum 和 cumsum使用方法）h=pi/40;x=0:h:pi/2;y=exp(-x).*sin(x);z1=sum(y(1:20)*h,z2=sum(y(2:21)*h,z=cumsum(y);z11=z(20)*h,z12=(z(21)-z(1)*h,运行后屏幕显示计算结果分别以下z1=z2=z11=z12=0.3873 0.4036 0.3873 0.4036求定积分准确值，输入程序 syms xF=int(exp(-x)*sin(x),x,0,pi/2)Fs=double(F),wz1=abs(Fs-z1),wz2=abs(Fs-z2)运行后屏幕显示定积分准确值Fs 和用矩形公式计算结果绝对误差wz1、wz2 分别以下F=Fs=1/2*(-1+exp(pi)(1/2)/exp(pi)(1/2)0.3961wz1=wz2=0.0088 0.0075第16页调用格式一：quad(fun,a,b)调用格式二：quad(fun,a,b,tol)调用格式三：Q,FCNT=quad(.)调用格式四：quad(fun,a,b,tol,TRACE)调用格式五：quad(fun,a,b,tol,TRACE,P1,P2,)复合辛普森（Simpson）数值积分MATLAB 主程序function y=comsimpson(fun,a,b,n)z1=feval(fun,a)+feval(fun,b);m=n/2;h=(b-a)/(2*m);x=a;z2=0;z3=0;x2=0;x3=0;for k=2:2:2*mx2=x+k*h;z2=z2+2*feval(fun,x2);endfor k=3:2:2*mx3=x+k*h;z3=z3+4*feval(fun,x3);endy=(z1+z2+z3)*h/3;1.3.3 辛普森公式Matlab程序第17页例5 用 comsimpson.m 和quad.m 分别计算定积分 ，取精度为10-4，并与准确值比较.解：输入程序 Q1,FCNT14=quad(fun,0,1,1.e-4,3),Q2=comsimpson(fun,0,1,10000)syms xfi=int(exp(-x.2)./2)./(sqrt(2*pi),x,0,1);Fs=double(fi)wQ1=double(abs(fi-Q1),wQ2=double(abs(fi-Q2)运行后屏幕显示I 准确值Fs，用comsimpson.m 和quad.m 分别计算I 近似值Q2、Q1和迭代次数FCNT14，取精度分别为104，Q2、Q1 分别与准确值Fs 绝对误差wQ2,wQ1第18页 MATLAB求解离散函数积分第19页 1、S=trapz(X,Y)等距梯形法求数值积分 调用格式一：Z=trapz(Y)调用格式二：Z=trapz(X,Y)调用格式三：Z=trapz(X,Y,DIM)或 trapz(Y,DIM)2、S=cumsum(Y)欧拉法求数值积分 调用格式一：cumsum(X)调用格式二：cumsum(X,DIM)3、对于离散点积分，先对其拟合取得函数表示式，再作为连续被积函数求积分。2 Matlab求解离散函数积分第20页例6 求 解：（1）输入程序 syms xF=limit(1/(sqrt(1-x2),x,1,left)运行后屏幕显示结果为F=inf 即当x 1-时，被积函数 （2）输入程序 syms xF1=int(1/(sqrt(1-x2),x,0,1),LimF1=double(F1)运行后屏幕显示结果为F1=LimF1=1/2*pi 1.5708 第21页例7 求解：（1）因为被积函数 在-1,1上除x=0外连续，输入程序 syms xF=limit(3*x5-8)/(x2),x,0)运行后屏幕显示结果为F=-inf（2）输入程序 syms xF1=int(3*x5-8)/(x2),x,-1,0),F2=int(3*x5-8)/(x2),x,0,1),F=F1+F2运行后屏幕显示结果为F1=-inf F2=35/4 F=-inf第22页例8 讨论反常积分 敛散性.解：（1）因为被积函数1/xq在(0,1上连续，在x=0无定义，输入 syms xLF05=limit(1/(x0.5),x,0,right),LF1=limit(1/x,x,0,right)LF2=limit(1/(x2),x,0,right)运行后屏幕显示结果为LF05=Inf LF1=Inf LF2=Inf（2）输入程序 syms xF05=int(1/(x0.5),x,0,1),F1=int(1/x,x,0,1),F2=int(1/(x2),x,0,1)运行后屏幕显示结果为F05=2 F1=Inf F2=Inf第23页MATLAB求解多重积分第24页3.1 二重积分数值求解 使用MATLAB提供dblquad函数就能够直接求出普通域上二重定积分数值解。该函数调用格式为：I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在a,bc,d区域上二重定积分。参数tol，trace使用方法与函数quad完全相同。3 Matlab求解多重积分第25页例9 用MATLAB 函数dblquad求直径为8 半球体积V球，误差为10-4.解:在 MATLAB 工作窗口输入以下MATLAB 程序 a=-4;b=4;c=-4;d=4;V1=dblquad(inline(sqrt(max(42-(x.2+y.2),0),a,b,c,d,1.e-4,quadl)V=dblquad(inline(sqrt(42-(x.2+y.2).*(x.2+y.2 x,y=meshgrid(-2:0.01:2);z1=8-(x.2+y.2);figure(1)meshc(x,y,z1)hold onx=-2:0.01:2;r=2;x,y,z=cylinder(r,30)mesh(x,y,z)hold offtitle(由旋转抛物面z=8-(x2+y2)，圆柱面 x2+y2=4和z=0 所围成积分区域V)figure(2)contour(x,y,z,10)title(由z=8-(x2+y2)，圆柱面 x2+y2=4和z=0 所围成区域V在x0y面上投影区域Dxy)运行后屏幕显示图形.第29页（2）确定积分限.输入程序 syms x y zf1=(z=8-(x2+y2);f2=(x2+y2=4);x,y,z=solve(f1,f2,x,y,z)运行后屏幕显示旋转抛物面和圆柱面交线以下x=y=z=(4-y2)(1/2)y 4-(4-y2)(1/2)y 4（3）输入计算程序 syms x y zf=x+exp(y)+sin(z);z1=0;z2=8-(x2+y2);x1=-sqrt(4-y2);x2=sqrt(4-y2);jfz=int(f,z,z1,z2);jfx=int(jfz,x,x1,x2);jfy=int(jfx,y,-2,2);jf2=double(jfy)运行后屏幕显示以下Warning:Explicit integral could not be found.In C:MATLAB6p5p1toolboxsymbolicsymint.m at line 58jf2=1.216650998803250e+002第30页Matlab多重积分（二）用MATLAB 数值计算三重积分调用格式一：Q3=triplequad(FUN,a,b,c,d,p,q)调用格式二：Q3=triplequad(FUN,a,b,c,d,p,q,tol)调用格式三：Q3=triplequad(FUN,a,b,c,d,p,q,tol,QUADL)调用格式四：Q3=triplequad(FUN,a,b,c,d,p,q,tol,MYQUADF)调用格式五：Q3=triplequad(FUN,a,b,c,d,p,q,tol,QUADL,P1,P2,.)调用格式六：Q3=triplequad(FUN,a,b,c,d,p,q,，,P1,P2,.)与 Q3=triplequad(FUN,a,b,c,d,p,q,1.e-6,QUAD,P1,P2,.)相同.第31页例11 分别用MATLAB 函数triplequad调用格式二和三计算 值，取误差限为tol=10-4，并将计算结果与准确值比较.其中V是三维长方体区域2x2,2y2,0z4.第32页解：建立并保留被积函数M 文件function u=integrnd1(x,y,z)u=x+exp(y)+sin(z);在MATLAB 工作窗口输入程序 a=-2;b=2;c=-2;d=2;p=0;q=4;Q3=triplequad(integrnd1,a,b,c,d,p,q,1.e-4)QL3=triplequad(integrnd1,a,b,c,d,p,q,1.e-4,quadl)syms x y zf=x+exp(y)+sin(z);jfz=int(f,z,p,q);jfy=int(jfz,y,c,d);jfx=int(jfy,x,a,b);I3=double(jfx),Juewu3=abs(I3-Q3)JuewuL3=abs(I3-QL3)第33页运行后屏幕显示以下Q3=1.425178451284647e+002,QL3=1.425178451284647e+002I3=1.425178309849224e+002,Juewu3=1.414354233020276e-005JuewuL3=1.414354233020276e-005第34页Thank You！第35页