各种数值积分方法的MATLAB程序和分析--以一个实际问题为例

各种数值积分方法的MATLAB 程序和分析 --以一个实际问题为例问题回顾从某条河的横截面上获得与河岸不同距离(y,m)处的深度数据(H,m)如下，确定这条河的平问题分析由于所得数据为离散点，并且数据点的间隔不均匀，因此先对离散数据进行连续化处理，再对所得函数进行数值积分。

离散数据点的分布图如下：一般来说，不考虑泥沙沉降因素的河底，用变分法来求取切力最小的河床面，可以知道椭圆是最理想的截面形状，然而根据上图可知，河床底部明显存在泥沙沉降。

并且该河段可能处于流向曲率较大的河段，导致内河床(x <5)泥沙沉降较多，河床较浅；而外河床(x >5)部分由于水流速度大于泥沙启动速度，河床底部泥沙被水流携带，因而较深。

根据泥沙启动概率与流速的关系可将泥沙沉降量与河岸距离近似处理为三次关系，如下图：河岸距离 m河道深度 m泥沙沉降量与离河岸距离的关系离河岸距离泥沙沉降量（负值为携带量)由于泥沙启动速度为止，那么将泥沙沉降和携带量达到平衡的地方作为拟合参数处理。

在理想河床截面形状的基础上以上述三次关系修正形状，则回归的函数模型如下：H=α1√1−α2(y−5)2−α3y(y−α4)(y−10)上式中前一项为理想河床截面的椭圆方程，后一项为泥沙沉降量的修正项，先利用非线性回归对上述数据进行拟合。

非线性回归采用Levenberg–Marquardt算法，其函数fitnonpoly.m 源代码如下：%本函数使用Levenberg–Marquardt算法计算待求数据点在要求函数类下的最佳拟合函数%输入：% x为待拟合数据点的横坐标，y为待拟合数据点的纵坐标% fun为自变量和待求参数的函数句柄，格式为fun(x,a)，其中a为参数向量% n为待求参数向量个数，即a长度% tol为结果要求精度%输出：% funchd为所得的非线性拟合函数，可直接使用function funchd=fitnonpoly(x,y,fun,n,tol)m=length(x);%输入检查if m~=length(y)error('The length of x must equal that of y !')elseif m<=n+1error('The length of x must be larger than n+1 !');end%变量创建a=zeros(1,n)+0.1;Z=zeros(m,n);delt=tol/10; %数值微分的扰动项delta=Inf*ones(1,n); %迭代的参数增量初始化v=5; %阻尼项衰减比lambda=10^-2; %阻尼项初始化flag=0;while norm(delta,Inf)>tol||flag==0 %迭代停止条件为参数增量小于要求精度%计算迭代矩阵方程for p=1:nTemp=zeros(1,n);Temp(p)=delt;Z(:,p)=(fun(x',a+Temp)-fun(x',a))./delt;endD=y'-fun(x',a);%计算迭代所得参数增量向量delta=(Z'*Z+lambda*eye(n))\(Z'*D);if norm(y'-fun(x',a+delta'))<norm(D)%若当前迭代收敛，则更新参数向量a=a+delta';lambda=lambda/v; %衰减阻尼项flag=1; elselambda=lambda*v; %若当前迭代不收敛，则增大阻尼 flag=0; end end%将函数表达式转换为可直接使用的函数句柄 syms tfunchd=@(t)fun(t,a); %依次显示所求得的参数值 disp('The parameter in turn is'); disp(a); end运行附件中脚本文件scriofint.m 可以得到拟合后的函数，拟合结果如下：非线性回归所得函数如下：H =1.782√1−0.3995(y −5)2−0.02014y(y −2.735)(y −10)可见泥沙沉降和携带量达到动态平衡的地方里河岸2.735m ,下面对河岸截面积进行数值积分求取。

河道深度与和河岸距离的关系(离散数据)河岸距离 m河道深度 m河道深度拟合残差河岸距离 m河道深度的残差 m复合梯形公式提高梯形法则的计算精度的一种方式是将[a,b]分为若干个子区间，在每个子区间上使用梯形法则，然后对所有子区间求和，得到整个区间上的积分值。

这样的积分公式成为复合梯形公式。

将积分区间n等分，积分区间内的n+1个基点为(x0,x1,x2,…,x n)，每个子区间的宽度为：ℎ=b−a n若取a=x0，b=x n，则总积分可表示为：I=∫f(x)dxx1x0+∫f(x)dxx2x1+⋯+∫f(x)dxx nx n−1将梯形法则应用于每一个积分项，得I=ℎf(x0)+f(x1)2+ℎf(x1)+f(x2)2+⋯+ℎf(x n−1)+f(x n)2合并同类项并分组，得I=(b−a)f(x0)+2∑f(x i)n−1i=1+f(x n)2n在上式中，对分子部分的f(x)的系数求和，再除以2n，结果正好是1。

因此，平均高度表示函数值的一种加权平均。

其中，每个内部点所对应的权值都是端点f(x0)和f(x n)所对应的权值的两倍。

复合梯形法求解数值积分的m函数文件compoundtrapz.m源代码如下：%本函数利用复合梯形公式求解数值积分%输入模式一：% y=compoundtrapz(xi,yi)% xi和yi为积分区间的自变量和因变量值，xi元素之间须等距，y为积分数值%输入模式二：% y=compoundtrapz(f,a,b,n)% f为待积函数句柄，a和b为积分下限和上限，n为积分区间等分数function y=compoundtrapz(varargin)if nargin==2 %输入判断xi=varargin{1};yi=varargin{2};n=length(xi)-1;elseif nargin==4xi=linspace(varargin{2},varargin{3},varargin{4}+1);yi=varargin{1}(xi);n=varargin{4};elseerror('Input error, please help for the m file !')y=(xi(end)-xi(1))*(sum(yi)+sum(yi(2:n)))/(2*n); %梯形公式 end运行附件中的scriofint.m 脚本文件可以得到如下结果： 河床横截面面积S 1=21.6023m 2，平均深度2。

16023m为了更好地了解到计算精度同区间等分数和舍入误差这两者的关系，本算例中使用单精度计算，并以I (f )−T 2n (f )≈13(T 2n (f )−T n (f ))作为误差估计。

运行脚本文件可以得到单精度下区间等分数与近似绝对误差的关系图。

结果为：从上图中可以看出，当区间等分数增大时，积分的数值结果与的近似绝对误差的对数近似呈指数关系下降，并且下降速率逐渐放缓，最终在等分区间数达到1000个以后，积分的数值结果的绝对误差不再下降，而是停滞在10−7.7左右。

这是由于当等分数增大的时候，区间两端点之间越来越趋近于一个直线段，以梯形逼近每一区间的面积就越来越精确，以至于当区间等分数增大的时候，积分的数值结果绝对误差越来越小。

然而由于本算例是使用单精度来计算，因此受限于舍入误差，单精度舍入误差必定是积分数值结果的下界，因此从上图可以看出，结果相对误差无法低于单精度下的机器精度10−8，即满足10−7.7<10−8。

为了量化截断误差与等分区间数的关系，下面对梯形公式的截断误差作详细推导。

根据一次插值及其误差项：f(a +α(b −a ))=f (a )+∆f (a )α+f ′′(ξ)2α(α−1)ℎ2 0≤α≤1则可将积分写成I =∫[f (a )+∆f (a )α+f ′′(ξ)2α(α−1)ℎ2]ba dx为了分析便利，根据α=(x −a)/ℎ，有dx =ℎdα因h =b −a ，积分区间由[a,b]变换为[0,1]，于是上式可表示为x 104单精度下复合梯形公式区间等分数与计算精度的关系(对数图)区间等分数 n近似绝对误差的对数 l o g 10(e p s a )I =ℎ∫[f (a )+∆f (a )α+f ′′(ξ)2α(α−1)ℎ2]1dα假设当h 很小时，f ′′(ξ)近似为常数，计算上式中的积分得I =ℎ[αf (a )+α22∆f (a )+(α36−α24)f ′′(ξ)ℎ2]01代入积分限，有I =ℎ[f (a )+∆f (a )2]−112f ′′(ξ)ℎ3 由∆f (a )=f (b )−f(a)，得到结果I =ℎf (a )+f(b)2−112f ′′(ξ)ℎ3 则可见截断误差e =112|f ′′(ξ)|ℎ3当使用复合型梯形公式时，由于积分区间被等分为n 个子区间，因此截断误差如下：e =112∑|f ′′(ξi )|ni=1(b −a n)3假设对于所有的i 都有|f ′′(ξi )|≤M ，那么复合梯形公式的截断误差为e =M (b −a )312n 2对于被积函数在积分区间内的二阶导数平均值可直接求导取两端之差的数值方法得到M =1.2300，另外令b −a =1代入到截断误差和等分区间关系式中并取对数： log 10e =log 10M (b −a )312n 2=log 101230.012n 2运行附件的scriofint.m 脚本文件亦可得如下关系图：从上图可以看出估计误差和截断误差随着区间等分数变化的总体趋势相同，在等份区间小于5000时，近似绝对误差和截断误差近似相等，这是由于在区间等分数较小的时候，舍入误差相对于截断误差来说很小，因此复合梯形公式的截断误差可以很好地估计实际绝对误差。

然而当区间等分数继续增大的时候，舍入误差随着区间等分数的增大而积累，使得相对x 104单精度下复合梯形公式区间等分数与近似绝对误差和截断误差的关系(对数图)区间等分数 n误差的对数 l o g 10(e p s )于总体误差来说，舍入误差变得不可忽略，因而在截断误差继续减小的时候，实际近似误差停滞在单精度的机器精度之上。

复合辛普森法除了对梯形法则的求积区间进行细化之外，另一种提高积分计算精度的方法是采用高阶多项式来连接数据点。

比方说在f(a)和f(b)之间添加一个中点，用抛物线将这三点连接起来。