1、第八章 积分与数值积分,一、符号积分,求符号积分函数:int,格式:int(f,x,a,b),功能:计算定积分,格式:int(f,x),功能:计算不定积分,使用int函数之前,先用syms声明x是符号变量,例8.1 求积分,syms x y1=1/(1+x4); y2=(x*exp(x)/(1+x)2; y3=1/(x2+2*x+3); fy1=int(y1) fy2=int(y2,0,1) fy3=int(y3,-inf,+inf),二、数值积分,在科学研究和工程技术中,经常遇到积分的计算,虽然有些函数的不定积分可以求出其初等函数表示式,但有更多的函数,它们的不定积分不是初等函数,这样就无法
2、利用牛顿莱布尼兹公式求出其定积分,如,甚至经常遇到只知道函数在一些离散点的值,但函数表达式未知的情况,在上述情况下就必须以数值方法求定积分的近似值。用数值方法求定积分的近似值,通常称为数值积分。,常用的数值积分算法有: 牛顿-柯特斯积分及其复化格式,自适用积分公式,龙贝格求积公式,高斯求积公式。Matlab中提供了求数值积分的函数:(1) quad函数:基于变步长辛普森法计算积分。该函数调用格式:I,n=quad(fname,a,b,Tol,trace)其中:fname是被积函数名 a,b是积分上下限 Tol是精度控制值,省却时取0.001 Trace:控制是否显示展现积分过程,取0不展现 I
3、:积分值 n:被积函数调用次数,如:求积分,ac=(x)sin(x)./xs=quad(ac,pi/4,pi/2),ac = (x)sin(x)./xs = 0.6118,(2) trapz函数:用梯形法计算积分,适用于被积函数为离散数据时,求函数的定积分。该函数调用格式:I=trapz(x,y),如:求积分,clc,clearformat longac=(x)sin(x)./xs1=quad(ac,pi/4,pi/2)x1=pi/4:pi/50:pi/2;y1=ac(x1);s2=trapz(x1,y1)x2=pi/4:pi/100:pi/2;y2=ac(x2);s3=trapz(x2,y2
4、),ac = (x)sin(x)./xs1 = 0.611786287055481s2 = 0.591535799623818s3 = 0.611773186493237,例8.2 神舟六号载人飞船与2005年10月12日发射升空后,在太空中运行115小时32分,绕地球飞行76圈,其先在轨道倾角42.4度、近地点高度200公里、远地点高度347公里的椭圆轨道上运行5圈,实施变轨后,进入343公里圆形轨道。试计算神州六号载人飞船在轨飞行公里数。,解 建立坐标系,使坐标原点位于椭圆中心,x轴位于椭圆长轴上。 s1=200, s2=347, 地球半径r=6371,椭圆长半轴,椭圆半焦距为c:,椭圆短
5、半轴:,椭圆方程:,已知 s1=200, s2=347, r=6371,s1=200, s2=347, r=6371,椭圆方程:,椭圆周长:,s1=200, s2=347, r=6371,s1=200;s2=347; r=6371;r2=343; a=(2*r+s1+s2)/2; c=(s2-s1)/2; k=c2/a2; g=(t)sqrt(1-k*cos(t).2); L=quad(g,0,pi/2) Lch=4*a*L*5+2*pi*(r+r2)*71,得:L1 =1.5707 Lch =3.2039e+006,例8.3 (面积的估算),测得某人工湖边界上n个点的坐标(见表),估算出该人
6、工湖的面积,表:边界上点实测坐标,作出人工湖的平面图:,clc,clear,load hqdatax=A(1,:),A(4,:),A(7,:);y1=A(2,:),A(5,:),A(8,:);y2=A(3,:),A(6,:),A(9,:);fill(x,y1,g)hold onfill(x,y2,g)hold off,(1)利用梯形法求面积,S=trapz(x,y2-y1); Mianji1=vpa(S,6),得:Mianji1 = 73330.8,(2)利用样条函数积分法 y=y2-y1; sc=spline(x,y); js=fnint(sc); %样条函数求积分 S=ppval(js,2
7、,493.1); mianji2=vpa(S(2),6),得:Mianji2 = 73586.0,例 8.4 在某山区的一个矩形区域内,测得一些地面点的高程(见下表),假设所观测的矩形区域内均是陆地,地表可近似看做一张光滑曲面,A(0,300,420.2),B(1200,300,332.5),C(800,1000,407.5)为地表上三个点,现需在A与B,A与C点之间的直线走向上沿地表铺设通讯线路,估计该通信线路的长度。,区域0,1200;0,1000内观测点处的高程表 (单位:米),记地表曲面方程为z=f(x,y),由于A与B的直线走向平行于x轴,,沿该走向的地表曲线LAB为:,曲线LAB的
8、弧长为:,问题求解: (1)通过插值法求出函数z=g(x) (2)通过数值积分求出该地表曲线的长度,地表高程数据存于文件:data804,clear,clcload data804x=0:100:1200;y=0:100:1000;g=z(4,:);pp=csape(x,g);dp=fnder(pp);gx=ppval(dp,x);hwf=sqrt(1+gx.2);s=trapz(x,hwf);dAB=vpa(s,8),得:dAB =1211.9855,A与B点的连线(y=300)的数据位于地面高程矩阵的第四行上,,记地表曲面方程为z=f(x,y),A(0,300,420.2)与C(800,1
9、000,407.5)的直线走向与x轴正向夹角为u,则:,记水平面上A1(0,300)与C1(800,1000)连线上任一点P(x,y)到点A1的距离为t,有,沿该走向的地表曲线LAC为:,沿该走向的地表曲线LAC为:,曲线LAC的弧长为:,load data804x=0:100:1200;y=0:100:1000;t=0:100*sqrt(113);xx=8/sqrt(113)*t;yy=300+7/sqrt(113)*t;x0,y0=meshgrid(x,y);h=interp2(x0,y0,z,xx,yy,spline);pp1=csape(t,h); %样条插值求函数h(t)dp1=fn
10、der(pp1); %求样条插值函数h(t)的导数ht=fnval(dp1,t);hwf=sqrt(1+ht.2); %求出被积函数s=trapz(t,hwf); %梯形法数值积分,求出弧长dAC=vpa(s,8),得:dAC =1069.3188,曲线LAC的弧长为:,画图:X,Y=meshgrid(x,y);xq=1:20:1200;yq=1:20:1000;Xq,Yq=meshgrid(xq,yq);Zq=interp2(X,Y,z,Xq,Yq);surf(Xq,Yq,Zq)hold onplot3(x,300+0*x,g,linewidth,2) %画曲线弧ABplot3(xx,yy,h,linewidth,2) %画曲线弧ACplot3(0,1200,800,300,300,1000,420.2,332.5,407.5,rp,linewidth,8) %画点A,B,C hold off,三、数值二重积分 I=dblquad(f,a,b,c,d,tol,trace) 求f(x,y)在a,bc,d区域上的二重积分。,例8.5 计算二重积分:,其中,f=(x,y)exp(-x.2/2).*sin(x.2+y) I=dblquad(f,-2,2,-1,1),得:I =1.5745,