数值积分的matlab实现

《数值积分的matlab实现》由会员分享，可在线阅读，更多相关《数值积分的matlab实现（10页珍藏版）》请在装配图网上搜索。

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

2、 - f (x *)h10.1)如果函数y二f(x)以列表形式给出(见表io-i),则其精确值无法求得，但可由下式求得 其近似值f (x * + h) - f (x *)f (x*)沁(10.2)h表 10-1xxxxxy0y-eiy12ynyny c2一般的，步长h越小，所得结果越精确。(10.2)式右端项的分子称为函数y二f(x)在x*的差分，分母称为自变量在x*的差分，所以右端项又称为差商。数值微分即用差商近似代替微商。常用的差商公式为：f(x +h)- f (x -h)f (x )沁 oo(10.3)O2 h八 x)3 y o ：4 yi- y 2(10.4)02hrf()y 4 y+

3、 3 yL、f (x )沁 n-2n1n(10.5)n2 h其误差均为O(h2)，称为统称三点公式。10.2数值微分的MATLAB实现MATLAB 提供了一个指令求解一阶向前差分，其使用格式为：dx=diff(x)其中x是n维数组，dx为n-1维数组L 一x ,x 一x ,x 一x ，这样基于两点的数值导2132 n 1数可通过指令diff(x)/h实现。对于三点公式，读者可参考例1的M函数文件diff3.m。例1用三点公式计算y二f (x)在x = 1.0,1.2,1.4处的导数值，f (x)的值由下表给出。x1.01.11.21.31.4f (x)0.25000.22680.20660.1

4、8900.1736解：建立三点公式的M函数文件diff3.m如下：function f=diff3(x,y)n=length(x);h=x(2)-x(1); f(1)=(-3*y(1)+4*y(2)-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,

5、 -0.1890, -0.1650, -0.0014。所以 y 二 f (x)在 x = 1.0,1.2,1.4 处的导数值分别为-0.2470, -0.1890 和-0.0014。对于高阶导数， MATLAB 提供了几个指令借助于样条函数进行求导，详细使用步骤如下step1：对给定数据点(x,y),利用指令pp=spline(x,y)，获得三次样条函数数据pp, 供后面ppval等指令使用。其中，pp是一个分段多项式所对应的行向量，它包含此多项式 的阶数、段数、节点的横坐标值和各段多项式的系数。step2:对于上面所求的数据向量pp，利用指令breaks,coefs,m,n二unmkpp(p

6、p)进行 处理，生成几个有序的分段多项式pp。step3：对各个分段多项式pp的系数，利用函数ppval生成其相应导数分段多项式的系 数，再利用指令mkpp生成相应的导数分段多项式step4：将待求点xx代入此导数多项式，即得样条导数值。上述过程可建立M函数文件ppd.m实现如下：function dy=ppd(pp)breaks,coefs,m=unmkpp(pp);for i=1:mcoefsm(i,:)=polyder(coefs(i,:);enddy=mkpp(breaks,coefsm);于是，如果已知节点处的值x,y,可用下面指令计算xx处的导数dyy：pp=spline(x,y)

7、,dy=ppd(pp);dyy=ppval(dy,xx);例2基于正弦函数ysin x的数据点，利用三点公式和三次样条插值分别求导，并与解析所求得的导数进行比较。解:编写 M 脚本文件 bijiao.m 如下：h=0.1*pi;x=0:h:2*pi;y=sin(x);dy1=diff3(x,y); pp=spline(x,y);dy=ppd(pp);dy2=ppval(dy,x); z=cos(x);error1=norm(dy1-z),error2=norm(dy2-z)plot(x,dy1,k:,x,dy2,r-,x,z,b)运行得结果为： error1 =0.0666， error2 =

8、0.0025，生成图形见图 10.1。图 10.1 三点公式、三次样条插值与解析求导比较图显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15 可知利用 三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点 附近误差较大，其它点吻合的较好。10.3 应用示例：湖水温度变化问题问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的 温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡 对某个湖的水温进行观测得数据见表10-2。表10-2某湖的水温观测数据深度（m）02.34.99.113.718.322.

9、927.2温度（C）22.822.822.820.613.911.711.111.1试找出湖水温度变化最大的深度。1问题的分析 湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导 数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。对于给定的数据， 可以利用数值微分计算各深度的温度变化值，从而得到温度变化最大的深度，但考虑到所给的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值 以得到更准确的结果。2模型的建立及求解记湖水的深度为h (m),相应的温度为T (D,且有T = T(h)，并假定函数T(h)可对给定的数据进行三次样条

10、插值，并对其求导，得到T(h)的插值导函数；然后将给定的深度数据加密，搜索加密数据的导数值的绝对值，找出其最大值及其相应的深度，相应的MATLAB 指令如下：h=0 2.3 4.9 9.1 13.7 18.3 22.9 27.2;T=22.8 22.8 22.8 20.6 13.9 11.7 11.1 11.1; hh=0:0.1:27.2;pp=spline(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运行得导函数绝对值的最大值点

11、为：h =11.4,最大值为1.6139,即湖水在深度为11.4m时温度变化最大，如图10.2所示(黑点为温度变化最大的点)。考虑定积分Jbf (x)dx(10.6)a如果被积函数f (x)是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项式P (x)近似地代替被积函数f (x)，然后计算积分JbP (x)dx，得(10.6)式的近似值；n a n如果被积函数的原函数不是初等函数，则将积分区间进行细分，对每个小区间，用一个近似函数代替被积函数f (x)，然后积分得(10.6)式的近似值。这两种类型最终都可归结为函数f(x)在节点xk上的函数值f(叮的某种线性组合，即下面数值求积公式:1

12、bf(x)gHAkf(xk)或10.7)k=0Jbf (x)dx =工A f (x ) + Rf (10.8)a k kk =0其中R f为截断误差。此误差可用代数精度衡量，代数精度越高，误差越小；反之误 差越大。代数精度是用来衡量数值积分公式近似程度的办法，如果f (x)是一个次数不超过m的 代数多项式，(10.7)式等号成立；而当f (x)是一个m +1次多项式时，(10.7)式不能精 确成立，则称(10.7)式的代数精度为m。选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普森公式 和高斯公式。10.5数值积分的MATLAB实现MATLAB 提供了下面几个函数计算积分

13、，其使用格式分别为：(1) trapz(x)采用梯形公式计算积分(h = 1), x为f (k = 0,1, ,n)k(2) quad(fun,a,b, tol)采用自适应Simpson法计算积分(3) quadl(fun,a,b,tol) 采用自适应 Gauss-Lobatto 法计算积分 其中fun为被积函数；tol是可选项，表示绝对误差，a，b为积分的上、下限。例1分别利用梯形公式、Simpson公式和Gauss-Lobatto法计算J ;1 + x2dx，并与其0 精确值比较。解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求 得的数值解比较，其MATLAB指令

14、为：syms xxz0=simple(int(sqrt(1+xxA2) ,0,1)z=double(z0);z=vpa(z,8)x=0:0.01:l;y=sqrt(l+x.A2);z1=trapz(y)*0.01;z1=vpa(z1,8),err1=z-z1;err1=vpa(err1,8)z2=quad(sqrt(1+x.A2), 0,l);z2=vpa(z2,8),err2=z-z2;err2=vpa(err2,8)z3=quadl(sqrt(1+x.A2), 0,l);z3=vpa(z3,8),err3=z-z3;err3=vpa(err3,8)运行得精确值为远-lnG 2 -1) =

15、1.1477936，三种公式计算得数值积分值分别为1.1477995，1.1477935 和 1.1477936，其相应误差分别为-.59e-5, .1e-6 和 0.，由三者误 差可见，Gauss-Lobatto法计算最为精确，Simpson公式次之，梯形公式最差，但它也能精 确到小数点后 5 位数。例 2人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球表 面439km，远地点距地球表面2384km，地球半径为6371km，求该卫星的轨道长度。解：卫星轨道椭圆的参数方程为x = a cos t, y = b sin t(0 t 2兀),a, b分别是长、短半轴，则根据所

16、给数据知 a =6371+2384=8755, b =6371+439=6810。 由对弧长的曲线积分知识知，椭圆的长度为f a1L = 4J 2 (a2 sin21 + b2 cos21)2dt0上积分称为椭圆积分，它无法用解析方法计算，可用计算其数值解，编写M函数文件如下: function y=y(t)a=8755;b=6810;y=4*sqrt(aA2*sin(t).A2+bA2*cos(t).A2);在 MATLAB 指令窗中输入以下指令: l=quad(y,0,pi/2)运行得结果为：l=4.9090e+004。即卫星轨道长度为49090km。对于用列表形式给出的函数，上述方法不再

17、适用，可利用指令spline构造三次样条插 值函数，再计算积分，具体步骤可参考例2。例3在桥梁的一端每隔一段时间记录lmin有几辆车过桥，得到数据见表10-3：表10-3过桥车辆数据时间0:002:004:005:006:007:008:00车辆数/辆22025825时间9:0010:3011:3012:3014:0016:0017:00车辆数/辆12510127928时间18:0019:0020:0021:0022:0023:0024:00车辆数/辆2210911893试估计一天通过桥梁的车流量。解：记记录时刻为t时，相应的车辆数为x(t)，则所求车流量即为计算积分24x(t)dt， 0则在

18、 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=2,2,0,2,5,8,25,12,5,10,12,7,9,28,22,10,9,11,8,9,3;pp=spline(x,y);s1=quadl(fun,0,24,pp)%利用三次样条插值计算积分s2=trapz(x,y)%利用梯形公式计算积分其中 M 函数文件 fun.m 为：function vf=fun(x,pp) vf=ppval(pp,x);运行得三次样条插值计算所得车流量为212辆，梯形公式计算所得车流量为21

19、6辆。10.6 数值积分的建模示例：煤炭储量计算问题问题：某煤矿为估计其煤炭的储量，在该矿区内进行勘探，得到数据如表 10-4。试估 算出该矿区(1 X 4,1 y 5 )煤炭的储量。表10-4某煤矿勘探数据表编号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.58

20、19.9523.7315.3518.0116.291问题的分析问题给出了很多点对应的煤炭厚度，显然整个煤矿可以看作是一个巨大的曲顶柱体(见 图 10.3，此图经过插值得到)，而煤炭的储量即为此立体的体积。要计算此立体的体积，可 以利用插值得到曲顶柱体的顶面函数，再对其积分；也可以将此曲顶柱体分割成若干个细的 曲顶柱体，用数值方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。10001ooo图 10.3 煤炭厚度曲面图2模型的建立及求解以煤炭的厚度为三维空间中的z坐标建立空间坐标系。记煤炭的厚度为h ,则它是坐标x,y的二元函数，即z(x,y)，则由二重积分的知识可知，此煤矿的煤炭储

21、量W为其中 D 二(x, y) 11 x 4,1 y 5。W = U 9 (x, y)dxdyD10.9)由于函数9 (x, y)只给出了一些离散点上的函数值，无法直接计算上述二重积分，所以下面采用数值积分的方法计算其值。由数值积分的知识知，计算定积分有复合梯形公式为Jbf (x)dx 沁ah f (a) + f (b) + 2艺 f (x ) 2kk=110.10)其中h为步长，x (k = 0,1,k由( 10.9)式得,n) 为节点，且有 x 二 a + kh。kJd 9 (x, y)dy dx 二Jb g(x)dxa10.11)其中 g(x) = idp(x,y)dy，则由(10.10

22、)式可得10.12)W = Jbg(x)dx Q g(a) + g(b) + 2艺g(x ) a 2 j aj =1g (a) = Jd p (a, y )dy q p (a, c) +p (a, d) +2kcck =1g (b) = jcd P (b, y )dy Q 2 P (b, C)+P (b, d ) + 逻 P (b, yk)k=1g(x ) =Jdp(x ,y)dy q y p(x ,c) +p(x ,d) + 2艺p(x ,y )kk2kkk kcck =1所以有p(a,c) + p(a,d) +p(b,c) +p(b,d) + 2迟p(a,y ) +p(b,y ) k k=

23、1hhW q 4+ 2艺p(x ,c) +p(x ,d) + 2迟p(x ,y )10.13)jjj kj=1k=1hh=x y -5p (a, c) - 5p (a, d) - 5p (b, c) - 5p (b, d)4+ 2区p(a,y ) +p(b,y )1+ 2p(x ,c) +p(x ,d) +kkjjj kk=0j=0j=0 k=0考虑到给定的数据较少，由此产生的误差较大，所以利用插值后的数据计算(10.13) 式，相应的 MATLAB 计算指令如下：x=1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4*1000;y=l 2 3 4 5 1 2 3

24、4 5 1 2 3 4 5 1 2 3 4 5*1000; z=-13.72 25.80 8.47 25.27 22.32 15.47 21.33 14.49 24.83 26.19 23.28 26.48 29.14 12.04 14.58 19.95 23.73 15.35 18.01 16.29;hx=10;hy=10;cx=1000:hx:4000;cy=1000:hy:5000;X,Y= meshgrid(cx,cy);n=length(cx);m=length(cy);Z=griddata(x,y,z,X,Y,v4); surf(X,Y,Z)%插值%绘制图 10.3W=-hx*hy

25、*(-5*Z(1,1)-5*Z(1,n)-5*Z(m,1)-5*Z(m,n)+2*(sum(Z(1,:)+Z(m,:)+sum(Z(:,1)+Z(:,n) +4*sum(sum(Z)/4运行得W =2 . 52 42X 108，即煤矿的煤炭储量约为2 . 52 42X 108m3.实验任务：1.一个物体的运动距离D = D(t)如下表所示:tD(t)8.09.010.011.012.017.45321.46025.75230.30135.084(1)利用三点公式和三次样条插值方法分别计算此物体在各个时刻的运动速率V(t)；(2)与D(t)的解析式D(t)二70 + 7t + 70e-io比较计

26、算结果。2已知20世纪美国人口统计数据如表10-5，试计算1900-1990年各年份的人口相对 增长率，并以此分析美国在这些年的人口变化情况。表10-5 20世纪美国人口统计数据年份1900191019201930194019501960197019801990人口( x 106)76.0 92.0106.5123.2131.7150.7179.3204.0226.5251.43.计算由表10-6数据给出的积分I 1.5y(x)dx。0.1表 10-6x0.10.30.50.70.91.11.31.5y1.01001.08731.22981.41501.61361.79431.92841.99

27、504.已知某铸件为曲顶柱体，现要对一个铸件的曲顶面进行涂漆，单位表面的费用为 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