数值积分的matlab实现

实验10 数值积分

实验目的：

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

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

0 d sin x x x

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

近似值，称为数值积分。

数值微分简介

设函数()y f x =在*

x 可导，则其导数为

h

x f h x f x f h )

()(lim )(**0*

-+='→ （）

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

h

x f h x f x f )

()()(***

-+≈' （）

表 10-1

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

x 的差分，分母称为自变量在*

x 的差分，所以右端项又称为差商。数值微分即用差商近似代替微商。常用的差商公式为：

000()()

()2f x h f x h f x h

+--'≈

（）

h y y y x f 243)(2

100-+-≈

' （）

h

y y y x f n

n n n 234)(12+-≈

'-- （）

其误差均为2

()O h ，称为统称三点公式。

数值微分的MATLAB 实现

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

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

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

解：建立三点公式的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=[,,,,];y=[,,,,];diff3(x,y)

运行得各点的导数值为：，，，，。所以()y f x =在=x ,,处的导数值分别为，和。

对于高阶导数，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 函数文件实现如下：

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 基于正弦函数sin y x 的数据点，利用三点公式和三次样条插值分别求导，并与解析所求得的导数进行比较。

解:编写M 脚本文件如下：

h=*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 =，error2 =，生成图形见图。

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

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

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

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

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

深度（m）0

温度（℃）

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

1．问题的分析

湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。对于给定的数据，可以利用数值微分计算各深度的温度变化值，从而得到温度变化最大的深度，但考虑到所给的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得到更准确的结果。