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

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

同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。

此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分⎰10 d sin x x x。

这时我们一般考虑用数值方法计算其近似值，称为数值积分。

10.1 数值微分简介设函数()y f x =在*x 可导，则其导数为hx f h x f x f h )()(lim )(**0*-+='→ （10.1）如果函数()y f x =以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得其近似值hx f h x f x f )()()(***-+≈' （10.2）表 10-1一般的，步长h 越小，所得结果越精确。

（10.2）式右端项的分子称为函数()y f x =在*x 的差分，分母称为自变量在*x 的差分，所以右端项又称为差商。

数值微分即用差商近似代替微商。

常用的差商公式为：000()()()2f x h f x h f x h +--'≈（10.3）hy y y x f 243)(2100-+-≈' （10.4）hy y y x f nn n n 234)(12+-≈'-- （10.5）其误差均为2()O h ，称为统称三点公式。

10.2 数值微分的MATLAB 实现MATLAB 提供了一个指令求解一阶向前差分，其使用格式为： dx=diff(x) 其中x 是n 维数组，dx 为1n -维数组[]21321,,,n x x x x x x ---,这样基于两点的数值导数可通过指令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-1f(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，-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:mcoefsm(i,:)=polyder(coefs(i,:)); enddy=mkpp(breaks,coefsm);于是，如果已知节点处的值x,y ，可用下面指令计算xx 处的导数dyy ：pp=spline(x,y),dy=ppd(pp);dyy=ppval(dy,xx);例2 基于正弦函数sin y 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 某湖的水温观测数据深度（m ） 0 2.3 4.9 9.1 13.7 18.3 22.9 27.2 温度（℃）22.822.822.820.613.911.711.111.1试找出湖水温度变化最大的深度。

1．问题的分析湖水的温度可视为关于深度的函数，于是湖水温度的变化问题便转化为温度函数的导数问题，显然导函数的最大绝对值所对应的深度即为温度变化最大的深度。

对于给定的数据，可以利用数值微分计算各深度的温度变化值，从而得到温度变化最大的深度，但考虑到所给的数据较少，由此计算的深度不够精确，所以采用插值的方法计算加密深度数据的导数值，以得到更准确的结果。

2．模型的建立及求解记湖水的深度为h （m ），相应的温度为T （℃），且有)(h T T =，并假定函数)(h T 可导。

对给定的数据进行三次样条插值，并对其求导，得到)(h T 的插值导函数；然后将给定的深度数据加密，搜索加密数据的导数值的绝对值，找出其最大值及其相应的深度，相应的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 数值积分简介考虑定积分()d baf x x ⎰（10.6）如果被积函数()f x 是以列表形式给出，则其求解思想同数值微分类似，即用逼近多项式()n P x 近似地代替被积函数()f x ，然后计算积分()d bn aP x x ⎰，得（10.6）式的近似值；如果被积函数的原函数不是初等函数，则将积分区间进行细分，对每个小区间，用一个近似函数代替被积函数()f x ，然后积分得（10.6）式的近似值。

这两种类型最终都可归结为函数()f x 在节点k x 上的函数值()k f x 的某种线性组合，即下面数值求积公式：()d ()nbk k ak I f x x A f x ==≈∑⎰ 或 （10.7）[]0()d ()nbk k ak f x x A f x R f ==+∑⎰（10.8）其中[]R f 为截断误差。

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

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

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

10.5 数值积分的MATLAB 实现MATLAB 提供了下面几个函数计算积分，其使用格式分别为：（1）trapz(x) 采用梯形公式计算积分（1=h ）,x 为),,1,0(n k f k = （2）quad('fun',a,b,tol) 采用自适应Simpson 法计算积分（3）quadl('fun',a,b,tol) 采用自适应Gauss-Lobatto 法计算积分 其中fun 为被积函数；tol 是可选项，表示绝对误差，a ，b 为积分的上、下限。

例1分别利用梯形公式、Simpson 公式和Gauss-Lobatto 法计算⎰+12d 1x x ，并与其精确值比较。

解：先对积分作符号运算，然后将其计算结果转换为数值型，再将其与这三种方法求得的数值解比较，其MATLAB 指令为：syms xxz0=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)运行得精确值为=--))12ln(2(211.1477936，三种公式计算得数值积分值分别为1.1477995，1.1477935和1.1477936，其相应误差分别为-.59e-5，.1e-6和0.，由三者误差可见，Gauss-Lobatto 法计算最为精确，Simpson 公式次之，梯形公式最差，但它也能精确到小数点后5位数。