实验10 数值积分
实验目的：
1．了解数值积分的基本原理；
2．熟练掌握数值积分的MATLAB实现；
3．会用数值积分方法解决一些实际问题。
实验内容：
积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来
直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所sinx 1?xd这时我们一般考虑用数值方法计算其如不定积分。以仍然得不到积分的精确值，x0
近似值，称为数值积分。10.1 数值微分简介
*x)f(x?y在设函数可导，则其导数为**)x?f(x(?h)f*?lim?(xf)）（10.1
h0h?y?f(x)以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得如果函数其近似值
**)(fx?h)f(x?*??xf)( 10.2）（h表 10-1
x
……n
hy?f(x)在（10.2）式右端项的分子称为函数一般的，步长越小，所得结果越精确。**xx的差分，所以右端项又称为差商。的差分，分母称为自变量在数值微分即用差商近似代替微商。常用的差商公式为：
f(x?h)?f(x?h)?00?(x)f（10.3 ）
02h?3y?4y?y?210?)(fx（10.4）0h2．
y?4y?3y?nn?2n?1?x)f(（10.5）n2h2)hO(，称为统称三点公式。其误差均为10.2 数值微分的MATLAB实现
MATLAB提供了一个指令求解一阶向前差分，其使用格式为：
dx=diff(x)
??nx?,xxx,x?,x?1?n,这样基于两点的数值导dx维数组，为其中x是维数组132n21数可通过指令diff(x)/h实现。对于三点公式，读者可参考例1的M函数文件diff3.m。
x?f()(xx)y?f的值由下表给用三点公式计算处的导数值，在1.0,1.2,1.41例
出。
1.0 1.1 1.2 1.3 1.40.2500 0.2268 0.2066 0.1890 0.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);
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)
y?f(x)所以，-0.0014。，运行得各点的导数值为：-0.2470，-0.2170-0.1890，
-0.1650x?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);
y?sinx的数据点，利用三点公式和三次样条插值分别求导，并与基于正弦函数例2解析所求得的导数进行比较。
解:编写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))'b''r--',x,z,plot(x,dy1,'k:',x,dy2,运行得结果为：error1 =0.0666，error2 =0.0025，生成图形见图10.1。
三点公式、三次样条插值与解析求导比较图图10.1
显然利用三次样条插值求导所得误差比三点公式求导小很多，同时由图2.15可知利用三次样条插值求导所得曲线与解析求导曲线基本重合，而三点公式在极值点附近和两个端点附近误差较大，其它点吻合的较好。
10.3 应用示例：湖水温度变化问题
问题：湖水在夏天会出现分层现象，其特点是接近湖面的水的温度较高，越往下水的温度越低。这种现象会影响水的对流和混合过程，使得下层水域缺氧，导致水生鱼类死亡。对某个湖的水温进行观测得数据见表10-2。
表10-2 某湖的水温观测数据
试找出湖水温度变化最大的深度。
1．问题的分析
湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。对于给定的数据，但考虑到所给从而得到温度变化最大的深度，可以利用数值微分计算各深度的温度变化值，
的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得到更准确的结果。
2．模型的建立及求解
ThT T(h)T(h)可（℃），相应的温度为，且有记湖水的深度为，并假定函数（m）导。
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)
'),grid on
plot(hh,dTT,',hh(i),dTT(i),'b'r.h=11.4，最大值为运行得导函数绝对值的最大值点为：1.6139,即湖水在深度为11.4m时温度变化最大，如图10.2所示（黑点为温度变化最大的点）。
