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

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

《整理数值积分的matlab实现.docx》由会员分享，可在线阅读，更多相关《整理数值积分的matlab实现.docx（16页珍藏版）》请在冰豆网上搜索。

整理数值积分的matlab实现

实验10数值积分

实验目的：

1•了解数值积分的基本原理；

2.熟练掌握数值积分的MATLAB实现；

3•会用数值积分方法解决一些实际问题。

实验内容：

积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。

同微分一样，在《微

积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所

以常常不能用来直接进行积分。

此外有些函数虽然有解析式，但其原函数不是初等函数，所

iSinX

以仍然得不到积分的精确值，如不定积分SUdX。

这时我们一般考虑用数值方法计算其

LOX

近似值，称为数值积分。

10.1数值微分简介

设函数y=f（χ）在X*可导，则其导数为

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

其近似值

「*f（x*+h）-f（x*）

h

表10-1

（10.2）

X

X0

X1X2

Xn

y

y。

y1y2

y∏

般的，步长h越小，所得结果越精确。

（10.2）式右端项的分子称为函数y=f（X）在

X的差分，分母称为自变量在X的差分，所以右端项又称为差商。

数值微分即用差商近似

代替微商。

常用的差商公式为:

（10.3）

（10.4）

其误差均为0（h），称为统称三点公式。

10.2数值微分的MATLA实现

MATLA醍供了一个指令求解一阶向前差分，其使用格式为：

dx=diff（x）

其中X是n维数组，dx为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）的值由下表给

出。

X

1.01.11.21.31.4

f（x）

0.25000.22680.20660.18900.1736

解：

建立三点公式的M函数文件diff3.m如下：

funCtiOnf=diff3（x,y）

n=Iength（x）;h=x

（2）-x

（1）;

f

（1）=（-3*y

（1）+4*y

（2）-y（3））∕（2*h）;

forj=2:

n-1f（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。

所以^f（X）

在^1.0,1.2,1.4处的导数值分别为-0.2470,-0.1890和-0.0014。

对于高阶导数，MATLA醍供了几个指令借助于样条函数进行求导，详细使用步骤如下：

step1:

对给定数据点（x,y）,利用指令PP=SPline（x,y）,获得三次样条函数数据PP,

供后面ppval等指令使用。

其中，PP是一个分段多项式所对应的行向量，它包含此多项式的阶数、段数、节点的横坐标值和各段多项式的系数。

step2:

对于上面所求的数据向量pp,利用指令[breaks,coefs,m,n]=unmkpp（pp）进行

处理，生成几个有序的分段多项式PPO

step3:

对各个分段多项式PP的系数，利用函数ppval生成其相应导数分段多项式的系

数，再利用指令mkpp生成相应的导数分段多项式

step4:

将待求点XX代入此导数多项式，即得样条导数值。

上述过程可建立M函数文件ppd.m实现如下：

funCtiondy=ppd（pp）

[breaks,coefs,m]=unmkpp（pp）;

fori=1:

m

COefSm（i,:

）=POIyder（COefS（i,:

））;

end

dy=mkpp（breaks,coefsm）;

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

PP=SPIine（x,y）,dy=ppd（pp）;dyy=ppval（dy,xx）;

例2基于正弦函数y=SinX的数据点，利用三点公式和三次样条插值分别求导，并与解析所求得的导数进行比较。

解：

编写M脚本文件bijiao.m如下：

h=0.1*pi;x=0:

h:

2*pi;y=sin（x）;

dy1=diff3（x,y）;

PP=SPIine（x,y）;dy=ppd（pp）;dy2=ppval（dy,x）;

Z=COS（x）;

error1=norm（dy1-z）,error2=norm（dy2-z）

PIOt（X,dy1,'k:

',x,dy2,'r--',x,z,'b'）

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

显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15可知利用

三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点

附近误差较大，其它点吻合的较好。

10.3应用示例：

湖水温度变化问题

问题：

湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。

这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。

对某个湖的水温进行观测得数据见表10-2。

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

深度（m）

0

2.3

4.9

9.1

13.7

18.3

22.9

27.2

温度「C）

22.8

22.8

22.8

20.6

13.9

11.7

11.1

11.1

试找出湖水温度变化最大的深度。

1•问题的分析

湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导

数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。

对于给定的数据，

可以利用数值微分计算各深度的温度变化值，从而得到温度变化最大的深度，但考虑到所给

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

2•模型的建立及求解

记湖水的深度为h（m）,相应的温度为T（C）,且有T=T（h）,并假定函数T（h）可导。

对给定的数据进行三次样条插值，并对其求导，得到T（h）的插值导函数；然后将给定

的深度数据加密，搜索加密数据的导数值的绝对值，找出其最大值及其相应的深度，相应的

MATLAB旨令如下：

h=[02.34.99.113.718.322.927.2];T=[22.822.822.820.613.911.711.111.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.'）,gridOn

运行得导函数绝对值的最大值点为：

h=11.4,最大值为1.6139,即湖水在深度为11.4m

时温度变化最大，如图10.2所示（黑点为温度变化最大的点）。

10.4数值积分简介

考虑定积分

b

f（x）dx（10.6）

a

如果被积函数f（x）是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项

b

式PI（X）近似地代替被积函数f（x）,然后计算积分.P1（x）dx,得（10.6）式的近似值；

a

如果被积函数的原函数不是初等函数，则将积分区间进行细分，对每个小区间，用一个近似

函数代替被积函数f（x）,然后积分得（10.6）式的近似值。

这两种类型最终都可归结为函数f（X）在节点Xk上的函数值f（Xk）的某种线性组合，即下面数值求积公式:

（10.7）

bn

I=Jf（x）dx止ΣAkf（Xk）a

k=0

bn

（10.8）

【f（x）dx=∑Akf（Xk）+R[f】a

k⊂0

其中Rf］为截断误差。

此误差可用代数精度衡量，代数精度越高，误差越小；反之误差越大。

代数精度是用来衡量数值积分公式近似程度的办法，如果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）quadl（'fun',a,b,tol）采用自适应Gauss-Lobatto法计算积分

其中fun为被积函数；tol是可选项，表示绝对误差，a，b为积分的上、下限。

12"

例1分别利用梯形公式、SimPSOn公式和Gauss-Lobatto法计算..Vx2dx,并与其

#0

精确值比较。

解：

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

SymSXX

ZO=Simple（int（'sqrt（1+xx^2）',0,1））

z=double（z0）;Z=VPa（Z,8）

x=09.01Jy=sqrt（1+x.^2）;

z1=trapz（y）*0.01;ZI=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）

运行得精确值为Ic2-1n（'∙.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=acost,y=bsint（0一t一2二）,a,b分别是长、短半

轴，则根据所给数据知a=6371+2384=8755,b=6371+439=6810。

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

π1

L=40（a2Sin2tb2cos2t）2dt

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

funCtiony=y（t）

a=8755;b=6810;

y=4*sqrt（a^2*sin（t）.^2+b^2*cos（t）.^2）;

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

l=quad（'y',O,pi/2）

运行得结果为：

l=4.9090e+004。

即卫星轨道长度为4909QkmO

对于用列表形式给出的函数，上述方法不再适用，可利用指令spline构造三次样条插

值函数，再计算积分，具体步骤可参考例2o

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

表10-3过桥车辆数据

时间

0:

00

2:

00

4:

00

5:

00

6:

00

7:

00

8:

00

车辆数/辆

2

2

0

2

5

8

25

时间

9:

00

10:

30

11:

30

12:

30

14:

00

16:

00

17:

00

车辆数/辆

12

5

10

12

7

9

28

时间

18:

00

19:

00

20:

00

21:

00

22:

00

23:

00

24:

00

车辆数/辆

22

10

9

11

8

9

3

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

24

解：

记记录时刻为t时，相应的车辆数为x（t）,则所求车流量即为计算积分Ox（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,10,12,7,9,28,22,10,9,11,8,9,3];

PP=SPIine（x,y）;s1=quadl（'fun',0,24,[],[],pp）%利用三次样条插值计算积分

s2=trapz（x,y）%利用梯形公式计算积分

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

funCtionvf=fun（x,pp）

Vf=PPVal（pp,x）;

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

10.6数值积分的建模示例：

煤炭储量计算问题

问题：

某煤矿为估计其煤炭的储量，在该矿区内进行勘探，得到数据如表10-4。

试估

算出该矿区（1乞X乞4,1乞y乞5）煤炭的储量。

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

编号

1

2

3

4

5

6

7

8

9

10

X坐标（km）

1

1

1

1

1

2

2

2

2

2

y坐标（km）

1

2

3

4

5

1

2

3

4

5

煤炭厚度h（m）

13.72

25.80

8.47

25.27

22.32

15.47

21.33

14.49

24.83

26.19

编号

11

12

13

14

15

16

17

18

19

20

X坐标（km）

3

3

3

3

3

4

4

4

4

4

y坐标（km）

1

2

3

4

5

1

2

3

4

5

煤炭厚度h（m）

23.28

26.48

29.14

12.04

14.58

19.95

23.73

15.35

18.01

16.29

1问题的分析

问题给出了很多点对应的煤炭厚度，显然整个煤矿可以看作是一个巨大的曲顶柱体（见

图10.3,此图经过插值得到），而煤炭的储量即为此立体的体积。

要计算此立体的体积，可以利用插值得到曲顶柱体的顶面函数，再对其积分；也可以将此曲顶柱体分割成若干个细的

曲顶柱体，用数值方法计算这些细曲顶柱体的体积，再对其求和即得原曲顶柱体的体积。

图10.3煤炭厚度曲面图

2•模型的建立及求解

以煤炭的厚度为三维空间中的Z坐标建立空间坐标系。

记煤炭的厚度为h,则它是坐标

x,y的二元函数，即ZhF:

（x,y）,则由二重积分的知识可知，此煤矿的煤炭储量W为

（10.9）

W=（X,y）dxdy

D

其中D={（X,y）|1乞Xm4,1乞y<5}。

由于函数（Xly）只给出了一些离散点上的函数值，无法直接计算上述二重积分，所以

F面采用数值积分的方法计算其值。

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

hn~1

f（x）dx：

-[f（a）f（b）2f（Xk）]

2k±i

其中h为步长，Xk（k=0,11IIn）为节点，且有X^=akh。

由（10.9）式得

（10.11）

LPCcP（X,y）dyIdX=ag（x）dx

d

其中g（χ）=[®（x,y）dy,则由（10.10）式可得

所以有

hym-1

g（xQ=CW（Xk,y）d"F[cP（Xk,c）N（Xk,d）+2∑CP（Xk,yQ]

2k=1

实验任务:

1.一个物体的运动距离D=D（t）如下表所示:

t

8.0

9.0

10.0

11.0

12.0

D（t）

17.453

21.460

25.752

30.301

35.084

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

（2）与D（t）的解析式D（t）二―707t70e」/1°比较计算结果。

2.已知20世纪美国人口统计数据如表10-5,试计算1900-1990年各年份的人口相对

增长率，并以此分析美国在这些年的人口变化情况。

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

年份

1900191019201930194019501960197019801990

人口（X106）

76.092.0106.5123.2131.7150.7179.3204.0226.5251.4

1.5

3•计算由表10-6数据给出的积分OIy（X）dx。

表10-6

X

0.1

0.3

0.5

0.7

0.9

1.1

1.3

1.5

y

1.0100

1.0873

1.2298

1.4150

1.6136

1.7943

1.9284

1.9950

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

表10-7铸件曲顶面数据

编号

1

2

3

4

5

6

7

8

9

10

X坐标

0

0

0

0

0

1.22

1.22

1.22

1.22

1.22

y坐标

0

2.14

4.28

6.42

8.56

0

2.14

4.28

6.42

8.56

Z坐标

1.20

1.20

1.20

1.20

1.20

1.20

1.70

0.33

2.20

0.35

编号

11

12

13

14

15

16

17

18

19

20

X坐标

2.44

2.44

2.44

2.44

2.44

3.66

3.66

3.66

3.66

3.66

y坐标

0

2.14

4.28

6.42

8.56

0

2.14

4.28

6.42

8.56

Z坐标

1.20

0.33

0.35

1.25

2.09

1.20

2.20

1.24

0.20

1.11