各种数值积分方法的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+1

error('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:n

Temp=zeros(1,n);

Temp(p)=delt;

Z(:,p)=(fun(x',a+Temp)-fun(x',a))./delt;

end

D=y'-fun(x',a);

%计算迭代所得参数增量向量

delta=(Z'*Z+lambda*eye(n))\(Z'*D);

if norm(y'-fun(x',a+delta'))

a=a+delta';

lambda=lambda/v; %衰减阻尼项

flag=1; else

lambda=lambda*v; %若当前迭代不收敛，则增大阻尼 flag=0; end end

%将函数表达式转换为可直接使用的函数句柄 syms t

funchd=@(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)dx

x1

x0+∫f(x)dx

x2

x1

+⋯+∫f(x)dx

x n

x 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−1

i=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==4

xi=linspace(varargin{2},varargin{3},varargin{4}+1);

yi=varargin{1}(xi);

n=varargin{4};

else

error('Input error, please help for the m file !')