数值积分与求解

数值积分与求解第一页，共三十二页，2022年，8月28日

MATLAB求解连续函数积分第二页，共三十二页，2022年，8月28日引言我们知道,若函数f(x)在区间[a,b]上连续且其原函数为F(x),则可用Newton-Leibnitz公式求定积分的值,Newton-Leibnitz公式无论在理论上还是在解决实际问题上都起了很大作用，但它并不能完全解决定积分的计算问题，因为积分学涉及的实际问题极为广泛，而且极其复杂，在实际计算中经常遇到以下三种情况：1Matlab求解连续函数积分第三页，共三十二页，2022年，8月28日

(1)被积函数f(x)并不一定能够找到用初等函数的有限形式表示的原函数F(x)，例如：Newton-Leibnitz公式就无能为力了(2)还有被积函数f(x)的原函数能用初等函数表示，但表达式太复杂，例如函数并不复杂，但积分后其表达式却很复杂，积分后其原函数F(x)为：第四页，共三十二页，2022年，8月28日(3)被积函数f(x)没有具体的解析表达式,其函数关系由表格或图形表示。对于这些情况,要计算积分的准确值都是十分困难的。由此可见,通过原函数来计算积分有它的局限性,因而研究一种新的积分方法来解决Newton-Leibniz公式所不能或很难解决的积分问题,这时需要用数值解法来建立积分的近似计算方法。

将积分区间细分,在每一个小区间内用简单函数代替复杂函数进行积分，这就是数值积分的思想，用代数插值多项式去代替被积函数发f(x)进行积分是本章讨论数值积分的主要内容。

第五页，共三十二页，2022年，8月28日1.1定积分的Matlab符号计算例1由y=sinx，y=cosx,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的程序>>symsxf1=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坐标调整第六页，共三十二页，2022年，8月28日1.2变限积分的Matlab符号计算例2已知,求F′(x)解：输入程序：>>symsxtF1=int(exp(t)*sin(2+sqrt(t^3)),x,0);F2=int(exp(t)*sin(2+sqrt(t^3)),0,x^2);Fi=F1+F2;dF=diff(Fi)运行后屏幕显示计算变限积分F(x)的导数.第七页，共三十二页，2022年，8月28日

建立数值积分公式的途径比较多,其中最常用的有两种：（1）由积分中值定理可知，对于连续函数f(x)，在积分区间[a,b]内存在一点ξ，使得即所求的曲边梯形的面积恰好等于底为(b-a)，高为的矩形面积。但是点ξ的具体位置一般是未知的,因而的值也是未知的,称为f(x)在区间[a,b]上的平均高度。那么只要对平均高度提供一种算法，相应地就获得一种数值求积方法1.3数值求积方法第八页，共三十二页，2022年，8月28日①梯形公式②矩形公式按照这种思想，可构造出一些求积分值的近似公式。例如分别取和则分别得到中矩形公式和梯形公式。第九页，共三十二页，2022年，8月28日③

Simpson公式

矩形公式把[a,b]的中点处函数值作为平均高度f()的近似值而获得的一种数值积分方法。

在这三个公式中,梯形公式把f(a),f(b)的加权平均值

作为平均高度f()的近似值而获得的一种数值积分方法。

第十页，共三十二页，2022年，8月28日Simpson公式是以函数f(x)在a,b,(a+b)/2这三点的函数值f(a),f(b),的加权平均值似值而获得的一种数值积分方法。作为平均高度f()的近（2）先用某个简单函数近似逼近f(x),用代替原被积函数f(x)，即以此构造数值算法。从数值计算的角度考虑,函数应对f(x)有充分的逼近程度,并且容易计算其积分。由于多项式能很好地逼近连续函数,且又容易计算积分,因此将选取为插值多项式,这样f(x)的积分就可以用其插值多项式的积分来近似代替第十一页，共三十二页，2022年，8月28日(一)用函数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程序第十二页，共三十二页，2022年，8月28日梯形数值积分的MATLAB主程序functionT=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;fori=1:mh=h/2;n=2*n;s=0;fork=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);第十三页，共三十二页，2022年，8月28日例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.第十四页，共三十二页，2022年，8月28日（一）函数sum的调用格式调用格式一：sum(X)调用格式二：sum(X,DIM)（二）函数cumsum的调用格式调用格式一：cumsum(X)调用格式二：cumsum(X,DIM)1.3.2矩形公式的Matlab程序第十五页，共三十二页，2022年，8月28日例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.38730.40360.38730.4036求定积分的精确值，输入程序>>symsxF=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.00880.0075第十六页，共三十二页，2022年，8月28日调用格式一：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主程序functiony=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;fork=2:2:2*mx2=x+k*h;z2=z2+2*feval(fun,x2);endfork=3:2:2*mx3=x+k*h;z3=z3+4*feval(fun,x3);endy=(z1+z2+z3)*h/3;1.3.3辛普森公式的Matlab程序第十七页，共三十二页，2022年，8月28日例5用comsimpson.m和quad.m分别计算定积分，取精度为10-4

，并与精确值比较.解：输入程序>>[Q1,FCNT14]=quad(@fun,0,1,1.e-4,3),Q2=comsimpson(@fun,0,1,10000)symsxfi=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第十八页，共三十二页，2022年，8月28日

MATLAB求解离散函数积分第十九页，共三十二页，2022年，8月28日

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、对于离散点的积分，先对其拟合获得函数表达式，再作为连续被积函数求积分。2Matlab求解离散函数积分第二十页，共三十二页，2022年，8月28日例6求解：（1）输入程序>>symsxF=limit(1/(sqrt(1-x^2)),x,1,'left')运行后屏幕显示结果为F=inf即当x→1-时，被积函数→∞（2）输入程序>>symsxF1=int(1/(sqrt(1-x^2)),x,0,1),LimF1=double(F1)运行后屏幕显示结果为F1=LimF1=1/2*pi1.5708第二十一页，共三十二页，2022年，8月28日例7求解：（1）因为被积函数在[-1,1]上除x=0外连续，输入程序>>symsxF=limit((3*x^5-8)/(x^2),x,0)运行后屏幕显示结果为F=-inf（2）输入程序>>symsxF1=int((3*x^5-8)/(x^2),x,-1,0),F2=int((3*x^5-8)/(x^2),x,0,1),F=F1+F2运行后屏幕显示结果为F1=-infF2=35/4F=-inf第二十二页，共三十二页，2022年，8月28日例8讨论反常积分

的敛散性.解：（1）因为被积函数1/xq在(0,1]上连续，在x=0无定义，输入>>symsxLF05=limit(1/(x^0.5),x,0,'right'),LF1=limit(1/x,x,0,'right')LF2=limit(1/(x^2),x,0,'right')运行后屏幕显示结果为LF05=InfLF1=InfLF2=Inf（2）输入程序>>symsxF05=int(1/(x^0.5),x,0,1),F1=int(1/x,x,0,1),F2=int(1/(x^2),x,0,1)运行后屏幕显示结果为F05=2F1=InfF2=Inf第二十三页，共三十二页，2022年，8月28日MATLAB求解多重积分第二十四页，共三十二页，2022年，8月28日3.1二重积分的数值求解

使用MATLAB提供的dblquad函数就可以直接求出一般域上二重定积分的数值解。该函数的调用格式为：

I=dblquad(f,a,b,c,d,tol,trace)该函数求f(x,y)在[a,b]×[c,d]区域上的二重定积分。参数tol，trace的用法与函数quad完全相同。3Matlab求解多重积分第二十五页，共三十二页，2022年，8月28日例9用MATLAB函数dblquad求直径为8的半球的体积V球，误差为10-4.解:在MATLAB工作窗口输入下列MATLAB程序>>a=-4;b=4;c=-4;d=4;V1=dblquad(inline('sqrt(max(4^2-(x.^2+y.^2),0))'),a,b,c,d,1.e-4,@quadl)V=dblquad(inline('sqrt(4^2-(x.^2+y.^2)).*(x.^2+y.^2<=4^2)'),a,b,c,d,1.e-4)symstwbjh=sqrt(4^2-(w.^2+t.^2));y1=-sqrt(4^2-t.^2);y2=sqrt(4^2-t.^2);jfx=int(bjh,w,y1,y2);jfy=int(jfx,t,a,b);I2=double(jfy),JuewuL1=abs(I2-V1)Juewu1=abs(I2-V),ezplot('x^2+y^2-16',[-5,5]);axisequal第二十六页，共三十二页，2022年，8月28日运行后屏幕显示如下V1=V=1.340434607608882e+0021.340477821376317e+002Warning:Explicitintegralcouldnotbefound.I2=JuewuL1=Juewu1=0.00649558446713第二十七页，共三十二页，2022年，8月28日3.2三重积分的数值求解(一)用MATLAB符号计算三重积分（1）画出积分区域V的草图.输入程序（2）确定积分限.输入程序（3）输入计算程序例10计算，其中积分区域V是由旋转抛物面Z=8-x2-y2，圆柱x2+y2=4和z=0所围成的空间闭区域.第二十八页，共三十二页，2022年，8月28日解：（1）画出积分区域V的草图.输入程序>>[x,y]=meshgrid(-2:0.01:2);z1=8-(x.^2+y.^2);figure(1)meshc(x,y,z1)holdonx=-2:0.01:2;r=2;[x,y,z]=cylinder(r,30)mesh(x,y,z)holdofftitle('由旋转抛物面z=8-(x^2+y^2)，圆柱面x^2+y^2=4和z=0所围成的积分区域V')figure(2)contour(x,y,z,10)title('由z=8-(x^2+y^2)，圆柱面x^2+y^2=4和z=0所围成区域V在x0y面上的投影区域Dxy')运行后屏幕显示图形.第二十九页，共三十二页，2022年，8月28日（2）确定积分限.输入程序>>symsxyzf1=('z=8-(x^2+y^2)');f2=('x^2+y^2=4');[x,y,z]=solve(f1,f2,x,y,z)运行后屏幕显示旋转抛物面和圆柱面的交线如下x=y=z=[(4-y^2)^(1/2)][