数值积分的matlab实现

实验10 数值积分

实验目的：

1．了解数值积分的基本原理； 2．熟练掌握数值积分的MATLAB 实现； 3．会用数值积分方法解决一些实际问题。

实验内容：

积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分近似值，称为数值积分。

sin x

⎰ 0x d x 。这时我们一般考虑用数值方法计算其

1

10.1 数值微分简介

设函数y =f (x ) 在x 可导，则其导数为

*

f (x *+h ) -f (x *)

f '(x ) =lim （10.1）

h →0h

*

如果函数y =f (x ) 以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得其近似值

f (x *+h ) -f (x *)

f '(x ) ≈ （10.2）

h

*

表 10-1

一般的，步长h 越小，所得结果越精确。（10.2）式右端项的分子称为函数y =f (x ) 在

x *的差分，分母称为自变量在x *的差分，所以右端项又称为差商。数值微分即用差商近似

代替微商。常用的差商公式为：

f '(x 0) ≈

f (x 0+h ) -f (x 0-h )

（10.3）

2h

f '(x 0) ≈

-3y 0+4y 1-y 2

（10.4）

2h

f '(x n ) ≈

其误差均为O (h 2) ，称为统称三点公式。

y n -2-4y n -1+3y n

（10.5）

2h

10.2 数值微分的MATLAB 实现

MATLAB 提供了一个指令求解一阶向前差分，其使用格式为： dx=diff(x)

其中x 是n 维数组，dx 为n -1维数组[x 2-x 1, x 3-x 2, , x n -x 1], 这样基于两点的数值导数可通过指令diff(x)/h实现。对于三点公式，读者可参考例1的M 函数文件diff3.m 。

例1 用三点公式计算y =f (x ) 在x =1.0,1.2,1.4处的导数值，f (x ) 的值由下表给

解：建立三点公式的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); end

f(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。所以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(pp)进行处理，生成几个有序的分段多项式pp 。

step3：对各个分段多项式pp 的系数，利用函数ppval 生成其相应导数分段多项式的系数，再利用指令mkpp 生成相应的导数分段多项式

step4：将待求点xx 代入此导数多项式，即得样条导数值。 上述过程可建立M 函数文件ppd.m 实现如下：

function dy=ppd(pp)

[breaks,coefs,m]=unmkpp(pp);

for i=1:m

coefsm(i,:)=polyder(coefs(i,:)); end

dy=mkpp(breaks,coefsm);

于是，如果已知节点处的值x,y ，可用下面指令计算xx 处的导数dyy ：

pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx);

例2 基于正弦函数y sin 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 =0.0025，生成图形见图10.1。

图10.1 三点公式、三次样条插值与解析求导比较图

显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15可知利用三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点附近误差较大，其它点吻合的较好。

10.3 应用示例：湖水温度变化问题

问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。对某个湖的水温进行观测得数据见表10-2。

表10-2 某湖的水温观测数据

试找出湖水温度变化最大的深度。 1．问题的分析

湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。对于给定的数据，可以利用数值微分计算各深度的温度变化值，从而得到温度变化最大的深度，但考虑到所给

的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得到更准确的结果。

2．模型的建立及求解

记湖水的深度为h （m ），相应的温度为T （℃），且有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 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

运行得导函数绝对值的最大值点为：h =11.4，最大值为1.6139, 即湖水在深度为11.4m 时温度变化最大，如图10.2所示（黑点为温度变化最大的点）。

图10.2 湖水温度变化曲线图

10.4 数值积分简介

考虑定积分

⎰

b

a

f (x )d x （10.6）

如果被积函数f (x ) 是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项式P n (x ) 近似地代替被积函数f (x ) ，然后计算积分

⎰

b

a

P n (x )d x ，得（10.6）式的近似值；

如果被积函数的原函数不是初等函数，则将积分区间进行细分，对每个小区间，用一个近似函数代替被积函数f (x ) ，然后积分得（10.6）式的近似值。这两种类型最终都可归结为函数f (x ) 在节点x k 上的函数值f (x k ) 的某种线性组合，即下面数值求积公式： I =

⎰

b

a

f (x )d x ≈∑A k f (x k ) 或 （10.7）

k =0

n

⎰

b

a

f (x )d x =∑A k f (x k ) +R [f ] （10.8）

k =0

n

其中R [f ]为截断误差。此误差可用代数精度衡量，代数精度越高，误差越小；反之误差越大。

代数精度是用来衡量数值积分公式近似程度的办法，如果f (x ) 是一个次数不超过m 的代数多项式，（10.7）式等号成立；而当f (x ) 是一个m +1次多项式时，（10.7）式不能精确成立，则称（10.7）式的代数精度为m 。

选取不同的近似函数，可产生不同的数值求积公式，常见的有：梯形公式、辛普森公式和高斯公式。

10.5 数值积分的MATLAB 实现

MATLAB 提供了下面几个函数计算积分，其使用格式分别为：

（1）trapz(x) 采用梯形公式计算积分（h =1）,x 为f k (k =0, 1, , n ) （2）quad('fun',a,b,tol) 采用自适应Simpson 法计算积分

（3）quadl('fun',a,b,tol) 采用自适应Gauss-Lobatto 法计算积分 其中fun 为被积函数；tol 是可选项，表示绝对误差，a ，b 为积分的上、下限。

例1分别利用梯形公式、Simpson 公式和Gauss-Lobatto 法计算

⎰

1

+x 2d x ，并与其

精确值比较。

解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求得的数值解比较，其MATLAB 指令为：

syms xx

z0=simple(int('sqrt(1+xx^2)',0,1)) z=double(z0);z=vpa(z,8) x=0:0.01:1;y=sqrt(1+x.^2);

z1=trapz(y)*0.01;z1=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)

运行得精确值为

1

(2-ln(2-1)) =1.1477936，三种公式计算得数值积分值分别为2

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 分别是长、短半轴，则根据所给数据知a =6371+2384=8755,b =6371+439=6810。

由对弧长的曲线积分知识知，椭圆的长度为

π

L =4

function y=y(t) a=8755;b=6810;

y=4*sqrt(a^2*sin(t).^2+b^2*cos(t).^2);

⎰

20

(a sin t +b cos t ) d t

2222

12

上积分称为椭圆积分，它无法用解析方法计算，可用计算其数值解，编写M 函数文件如下：

在MATLAB 指令窗中输入以下指令：

l=quad('y',0,pi/2)

运行得结果为：l=4.9090e+004。即卫星轨道长度为49090km 。

对于用列表形式给出的函数，上述方法不再适用，可利用指令spline 构造三次样条插值函数，再计算积分，具体步骤可参考例2。

例3 在桥梁的一端每隔一段时间记录1min 有几辆车过桥，得到数据见表10-3：

表10-3 过桥车辆数据

试估计一天通过桥梁的车流量。

解：记记录时刻为t 时，相应的车辆数为x (t ) ，则所求车流量即为计算积分则在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) %利用梯形公式计算积分

⎰

24

x (t ) d t ，

其中M 函数文件fun.m 为：

function vf=fun(x,pp) vf=ppval(pp,x);

运行得三次样条插值计算所得车流量为212辆，梯形公式计算所得车流量为216辆。

10.6 数值积分的建模示例：煤炭储量计算问题

问题：某煤矿为估计其煤炭的储量，在该矿区内进行勘探，得到数据如表10-4。试估算出该矿区（1≤x ≤4, 1≤y ≤5）煤炭的储量。

表10-4 某煤矿勘探数据表

1．问题的分析

问题给出了很多点对应的煤炭厚度，显然整个煤矿可以看作是一个巨大的曲顶柱体（见图10.3，此图经过插值得到），而煤炭的储量即为此立体的体积。要计算此立体的体积，可以利用插值得到曲顶柱体的顶面函数，再对其积分；也可以将此曲顶柱体分割成若干个细的曲顶柱体，用数值方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。

图10.3 煤炭厚度曲面图 2．模型的建立及求解

以煤炭的厚度为三维空间中的z 坐标建立空间坐标系。记煤炭的厚度为h ，则它是坐标

x , y 的二元函数，即z =ϕ(x , y ) ，则由二重积分的知识可知，此煤矿的煤炭储量W 为

W =其中D ={(x , y ) |1≤x ≤4, 1≤y ≤5}。

由于函数ϕ(x , y ) 只给出了一些离散点上的函数值，无法直接计算上述二重积分，所以下面采用数值积分的方法计算其值。

由数值积分的知识知，计算定积分有复合梯形公式为

⎰⎰ϕ(x , y ) d x d y （10.9）

D

⎰

b

a

n -1

h

f (x ) d x ≈[f (a ) +f (b ) +2∑f (x k ) ] （10.10）

2k =1

其中h 为步长，x k (k =0, 1, , n ) 为节点，且有x k =a +kh 。

由（10.9）式得

W =

⎰

b

a

⎡d ϕ(x , y ) d y ⎤d x =b g(x)d x （10.11）

⎰a ⎢⎥⎣⎰c ⎦

其中g (x ) =

⎰ϕ(x , y ) d y ，则由（10.10）式可得

c

d

W =而

g (a ) =

⎰

d

b

a

n -1

h x

g (x ) d x ≈[g (a ) +g (b ) +2∑g (x j ) ] （10.12）

2j =1

⎰ϕ(a , y ) d y ≈

c d c

h y 2h y 2

[ϕ(a , c ) +ϕ(a , d ) +2∑ϕ(a , y k ) ]

k =1

m -1

g (b ) =⎰ϕ(b , y ) d y ≈

d

[ϕ(b , c ) +ϕ(b , d ) +2∑ϕ(b , y k ) ]

k =1

m -1

g (x k ) =⎰ϕ(x k , y ) d y ≈

c

h y 2

[ϕ(x k , c ) +ϕ(x k , d ) +2∑ϕ(x k , y k ) ]

k =1

m -1

所以有

W ≈

h x h y 4

[ϕ(a , c ) +ϕ(a , d ) +ϕ(b , c ) +ϕ(b , d ) +2∑[ϕ(a , y k ) +ϕ(b , y k ) ]

k =1

m -1k =1

m -1

+2∑[ϕ(x j , c ) +ϕ(x j , d ) +2∑ϕ(x j , y k ) ]]

j =1

n -1

=

h x h y 4

m k =0

（10.13）

[-5ϕ(a , c ) -5ϕ(a , d ) -5ϕ(b , c ) -5ϕ(b , d )

n

n

m

+2∑[ϕ(a , y k ) +ϕ(b , y k ) ]+2∑[ϕ(x j , c ) +ϕ(x j , d )]+4∑∑ϕ(x j , y k )

j =0

j =0k =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=[1 2 3 4 5 1 2 3 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.3

W=-hx*hy*(-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.5242×10，即煤矿的煤炭储量约为2.5242×10m .

8

83

实验任务：

1．一个物体的运动距离D =D (t ) 如下表所示：

（1）利用三点公式和三次样条插值方法分别计算此物体在各个时刻的运动速率V (t ) ；

（2）与D (t ) 的解析式D (t ) =-70+7t +70e -t /10比较计算结果。

2．已知20世纪美国人口统计数据如表10-5，试计算1900-1990年各年份的人口相对增长率，并以此分析美国在这些年的人口变化情况。

表10-5 20世纪美国人口统计数据

3．计算由表10-6数据给出的积分

⎰

1. 5

0. 1

y (x ) d x 。

表10-6

4．已知某铸件为曲顶柱体，现要对一个铸件的曲顶面进行涂漆，单位表面的费用为100元，该铸件曲顶面数据如表10-7所示，试估算该项处理的费用。

表10-7 铸件曲顶面数据