整理数值积分的matlab实现.docx

1、整理数值积分的matlab实现实验10数值积分实验目的：1了解数值积分的基本原理；2.熟练掌握数值积分的 MATLAB实现；3会用数值积分方法解决一些实际问题。实验内容：积分是数学中的一个基本概念， 在实际问题中也有很广泛的应用。 同微分一样，在微积分中，它也是通过极限定义的， 由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。 此外有些函数虽然有解析式， 但其原函数不是初等函数， 所i Sin X以仍然得不到积分的精确值， 如不定积分 SUdX。这时我们一般考虑用数值方法计算其LO X近似值，称为数值积分。10.1数值微分简介设函数y = f()在X*可导，则其导数

2、为如果函数y = f (X)以列表形式给出(见表 10-1),则其精确值无法求得，但可由下式求得其近似值* f(x*+h)-f(x*)h表 10-1(10.2)XX0X1 X2 Xnyy。y1 y2 y般的，步长h越小，所得结果越精确。(10.2 )式右端项的分子称为函数 y = f (X)在X的差分，分母称为自变量在 X的差分，所以右端项又称为 差商。数值微分即用差商近似代替微商。常用的差商公式为:(10.3 )(10.4 )其误差均为0(h )，称为统称三点公式。10.2数值微分的MATLA实现MATLA醍供了一个指令求解一阶向前差分，其使用格式为：dx=diff(x)其中X是n维数组，d

3、x为n-1维数组&2 - X1,X3 - X2,IH, Xn - X1丨，这样基于两点的数值导数可通过指令diff(x)/h实现。对于三点公式，读者可参考例 1的M函数文件diff3.m 。例1用三点公式计算 y = f(x)在X =1.0,1.2,1.4 处的导数值，f(x)的值由下表给出。X1.0 1.1 1.2 1.3 1.4f(x)0.2500 0.2268 0.2066 0.1890 0.1736解：建立三点公式的 M函数文件diff3.m 如下：fun Cti On f=diff3(x,y)n=Ie ngth(x);h=x(2)-x(1);f(1)=(-3*y(1)+4*y (2)

4、-y(3)(2*h);for j=2: n-1 f(j)=(y(j+1)-y(j-1)/(2*h);endf(n )=(y( n-2)-4*y( n-1)+3*y( n)(2*h);在MATLAB旨令窗中输入指令：x=1.0,1.1,1.2,1.3,1.4;y=0.2500,0.2268,0.2066,0.1890,0.1736;diff3(x,y)运行得各点的导数值为：-0.2470 ,-0.2170 , -0.1890 , -0.1650 , -0.0014。所以 f(X)在 1.0,1.2,1.4 处的导数值分别为-0.2470 , -0.1890 和-0.0014。对于高阶导数，MAT

5、LA醍供了几个指令借助于样条函数进行求导，详细使用步骤如下：step1 :对给定数据点(x,y ),利用指令PP=SPline(x,y),获得三次样条函数数据 PP ,供后面ppval等指令使用。其中，PP是一个分段多项式所对应的行向量，它包含此多项式 的阶数、段数、节点的横坐标值和各段多项式的系数。step2 :对于上面所求的数据向量 pp,利用指令breaks,coefs,m,n=unmkpp(pp) 进行处理，生成几个有序的分段多项式 PPOstep3 :对各个分段多项式 PP的系数，利用函数ppval生成其相应导数分段多项式的系数，再利用指令mkpp生成相应的导数分段多项式step4

6、:将待求点XX代入此导数多项式，即得样条导数值。上述过程可建立 M函数文件ppd.m实现如下：fun Cti on dy=ppd(pp)breaks,coefs,m=u nmkpp(pp);for i=1:mCOefSm(i,:)=POIyder(COefS(i,:);enddy=mkpp(breaks,coefsm);于是，如果已知节点处的值 x,y ,可用下面指令计算 XX处的导数dyy:PP=SPIi ne(x,y),dy=ppd(pp);dyy=ppval(dy,xx);例2基于正弦函数y = Sin X的数据点，利用三点公式和三次样条插值分别求导，并与 解析所求得的导数进行比较。解：

7、编写M脚本文件bijiao.m 如下：h=0.1*pi;x=0:h:2*pi;y=si n(x);dy1=diff3(x,y);PP=SPIi ne(x,y);dy=ppd(pp);dy2=ppval(dy,x);Z=COS(x);error1= no rm(dy1-z),error2=n orm(dy2-z)PIOt(X,dy1, k: ,x,dy2, r- ,x,z, b)图10.1三点公式、三次样条插值与解析求导比较图显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图 2.15可知利用三次样条插值求导所得曲线与解析求导曲线基本重合， 而三点公式在极值点附近和两个端点附近误差较

8、大，其它点吻合的较好。10.3应用示例：湖水温度变化问题问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的 温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。 对某个湖的水温进行观测得数据见表 10-2。表10-2某湖的水温观测数据深度（m）02.34.99.113.718.322.927.2温度C）22.822.822.820.613.911.711.111.1试找出湖水温度变化最大的深度。1 问题的分析湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深

9、度。 对于给定的数据，可以利用数值微分计算各深度的温度变化值， 从而得到温度变化最大的深度，但考虑到所给的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值， 以得到更准确的结果。2 模型的建立及求解记湖水的深度为h (m),相应的温度为T (C),且有T=T(h),并假定函数T(h)可 导。对给定的数据进行三次样条插值，并对其求导，得到 T(h)的插值导函数；然后将给定的深度数据加密，搜索加密数据的导数值的绝对值， 找出其最大值及其相应的深度， 相应的MATLAB旨令如下：h=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=22.8 22.8

10、 22.8 20.6 13.9 11.7 11.1 11.1; hh=0:0.1:27.2;PP=SPli ne(h,T);dT=ppd(pp);dTT=ppval(dT,hh);dTTmax,i=max(abs(dTT),hh(i)plot(hh,dTT, b ,hh(i),dTT(i), r. ),grid On运行得导函数绝对值的最大值点为： h=11.4 ,最大值为1.6139,即湖水在深度为11.4m时温度变化最大，如图 10.2所示(黑点为温度变化最大的点)。10.4数值积分简介考虑定积分bf (x)dx (10.6 )a如果被积函数f (x)是以列表形式给出，则其求解思想同数值微

11、分类似，即用逼近多项b式PI(X)近似地代替被积函数 f (x),然后计算积分.P1(x)dx ,得(10.6 )式的近似值；a如果被积函数的原函数不是初等函数， 则将积分区间进行细分，对每个小区间，用一个近似函数代替被积函数 f (x),然后积分得(10.6 )式的近似值。这两种类型最终都可归结为函 数f(X)在节点Xk上的函数值f (Xk)的某种线性组合，即下面数值求积公式:(10.7 )b nI = J f (x)dx 止 Akf(Xk) ak=0b n(10.8 )【f(x)dx= Akf(Xk)+Rf】 ak0其中Rf 为截断误差。此误差可用代数精度衡量，代数精度越高，误差越小；反之

12、误 差越大。代数精度是用来衡量数值积分公式近似程度的办法， 如果f (X)是一个次数不超过 m的代数多项式，(10.7 )式等号成立；而当 f (x)是一个m 1次多项式时，(10.7)式不能精确成立，则称(10.7 )式的代数精度为 m。选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普森公式 和高斯公式。10.5数值积分的MATLA实现MATLA醍供了下面几个函数计算积分，其使用格式分别为：(1) trapz(x) 采用梯形公式计算积分(h=1),x为fk(k=0,1,，n)(2) quad(fun,a,b,tol) 采用自适应 SimPSOn法计算积分(3) quad

13、l(fun,a,b,tol) 采用自适应 Gauss-Lobatto 法计算积分其中fun为被积函数；tol是可选项，表示绝对误差， a，b为积分的上、下限。1 2例1分别利用梯形公式、SimPSOn公式和Gauss-Lobatto法计算 .V x2dx ,并与其#0精确值比较。解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求 得的数值解比较，其 MATLAB旨令为：SymS XXZO=Simple(i nt(sqrt(1+xx2),0,1)z=double(z0);Z=VPa(Z,8)x=09.01Jy=sqrt(1+x.2);z1=trapz(y)*0.01;Z

14、I=VPa(Z1,8),err1=z-z1;err1=vpa(err1,8)z2=quad(sqrt(1+x.2), 0,1);z2=vpa(z2,8),err2=z-z2;err2=vpa(err2,8)z3=quadl(sqrt(1+x.2), 0,1);z3=vpa(z3,8),err3=z-z3;err3=vpa(err3,8)运行得精确值为 Ic 2 -1 n(.2 -1) =1.1477936 ,三种公式计算得数值积分值分别为21.1477995 , 1.1477935 和 1.1477936 ,其相应误差分别为 -.59e-5 , .1e-6 和 0,由三者误 差可见，Gauss

15、-Lobatto法计算最为精确，SimPSOn公式次之，梯形公式最差，但它也能精 确到小数点后5位数。例2人造地球卫星轨道可视为平面上的椭圆。 我国第一颗人造地球卫星近地点距地球表面439km,远地点距地球表面 2384km,地球半径为6371km,求该卫星的轨道长度。解：卫星轨道椭圆的参数方程为 X =acost, y =bsin t(0 一 t 一 2二),a,b分别是长、短半轴，则根据所给数据知 a =6371+2384=8755, b =6371+439=6810。由对弧长的曲线积分知识知，椭圆的长度为 1L = 4 0 (a2 Sin2t b2 cos2 t)2dt上积分称为椭圆积分

16、，它无法用解析方法计算，可用计算其数值解，编写 M函数文件如下:fun Cti on y=y(t)a=8755;b=6810;y=4*sqrt(a2*si n( t).2+b2*cos(t).2);在MATLAB旨令窗中输入以下指令：l=quad(y,O,pi/2)运行得结果为：l=4.9090e+004。即卫星轨道长度为 4909QkmO对于用列表形式给出的函数，上述方法不再适用，可利用指令 spl ine构造三次样条插值函数，再计算积分，具体步骤可参考例 2o例3在桥梁的一端每隔一段时间记录 Imin有几辆车过桥，得到数据见表 10-3 :表10-3过桥车辆数据时间0:002:004:00

17、5:006:007:008:00车辆数/辆22025825时间9:0010:3011:3012:3014:0016:0017:00车辆数/辆12510127928时间18:0019:0020:0021:0022:0023:0024:00车辆数/辆2210911893试估计一天通过桥梁的车流量。24解：记记录时刻为t时，相应的车辆数为 x(t),则所求车流量即为计算积分 O x(t)dt ,则在MATLAB旨令窗中输入下面指令：x=0,2,4,5,6,7,8,9,10.3,11.3,12.3,14,16,17,18,19,20,21,22,23,24; y=220,2,5,8,25,12,5,1

18、0,12,7,9,28,22,10,9,11,8,9,3;PP=SPIi ne(x,y);s1=quadl(fu n ,0,24,pp) % 利用三次样条插值计算积分s2=trapz(x,y) %利用梯形公式计算积分其中M函数文件fun.m为：fun Cti on vf=fu n( x,pp)Vf=PPVal(pp,x);运行得三次样条插值计算所得车流量为 212辆，梯形公式计算所得车流量为 216辆。10.6数值积分的建模示例：煤炭储量计算问题问题：某煤矿为估计其煤炭的储量，在该矿区内进行勘探，得到数据如表 10-4。试估算出该矿区(1乞X乞4,1乞y乞5 )煤炭的储量。表10-4某煤矿勘探

19、数据表编号12345678910X坐标(km)1111122222y坐标(km)1234512345煤炭厚度h ( m)13.7225.808.4725.2722.3215.4721.3314.4924.8326.19编号11121314151617181920X坐标(km)3333344444y坐标(km)1234512345煤炭厚度h ( m)23.2826.4829.1412.0414.5819.9523.7315.3518.0116.291问题的分析问题给出了很多点对应的煤炭厚度， 显然整个煤矿可以看作是一个巨大的曲顶柱体 (见图10.3 ,此图经过插值得到)，而煤炭的储量即为此立体的

20、体积。要计算此立体的体积，可 以利用插值得到曲顶柱体的顶面函数， 再对其积分；也可以将此曲顶柱体分割成若干个细的曲顶柱体，用数值方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。图10.3煤炭厚度曲面图2 模型的建立及求解以煤炭的厚度为三维空间中的 Z坐标建立空间坐标系。记煤炭的厚度为h,则它是坐标x, y的二元函数，即ZhF:(x, y),则由二重积分的知识可知，此煤矿的煤炭储量 W为(10.9 )W= (X,y)dxdyD其中 D =( X, y) |1 乞 X m 4,1 乞 y 5。由于函数 (Xl y)只给出了一些离散点上的函数值，无法直接计算上述二重积分，所以F面采用数

21、值积分的方法计算其值。由数值积分的知识知，计算定积分有复合梯形公式为h n 1f (x)dx ： - f (a) f (b) 2 f (Xk)2 ki其中h为步长，Xk (k = 0,11 I In)为节点，且有X= a kh。由(10.9 )式得(10.11)LPCcP(X,y)dy IdX= a g(x)dxd其中 g() = (x, y)dy ,则由(10.10)式可得所以有hy m-1g(xQ = C W(Xk,y)dFcP(Xk,c)N(Xk,d)+2 CP(Xk,yQ2 k=1实验任务:1.一个物体的运动距离 D= D(t)如下表所示:t8.09.010.011.012.0D(t)

22、17.45321.46025.75230.30135.084(1)利用三点公式和三次样条插值方法分别计算此物体在各个时刻的运动速率 V(t)；(2)与D(t)的解析式D(t)二70 7t 70e/1比较计算结果。2.已知20世纪美国人口统计数据如表 10-5 ,试计算1900-1990年各年份的人口相对增长率，并以此分析美国在这些年的人口变化情况。表10-5 20世纪美国人口统计数据年份1900 1910 1920 1930 1940 1950 1960 1970 1980 1990人口(X 106)76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 2

23、26.5 251.41.53 计算由表10-6数据给出的积分 OIy(X)dx。表 10-6X0.10.30.50.70.91.11.31.5y1.01001.08731.22981.41501.61361.79431.92841.99504 已知某铸件为曲顶柱体，现要对一个铸件的曲顶面进行涂漆，单位表面的费用为 100元，该铸件曲顶面数据如表 10-7所示，试估算该项处理的费用。表10-7铸件曲顶面数据编号12345678910X坐标000001.221.221.221.221.22y坐标02.144.286.428.5602.144.286.428.56Z坐标1.201.201.201.201.201.201.700.332.200.35编号11121314151617181920X坐标2.442.442.442.442.443.663.663.663.663.66y坐标02.144.286.428.5602.144.286.428.56Z坐标1.200.330.351.252.091.202.201.240.201.11