数值积分：梯形积分公式和Simpson积分公式

Posted on

Background

在实际工程中许多函数的原函数并不容易求出，因此对它们使用Newton-Lebniz公式是相当困难的，因此衍生出数值积分，基本思想就是对于一个黎曼可积分的函数，不断细分积分区间对小面积求和得到定积分的近似值；

但是过于细致的区间划分会导致大量资源的消耗，因此这篇文章旨在介绍一种启发式的办法来使得误差控制在我们接收范围内同时减少对计算资源的消耗；

由于对数值分析的了解还不够，不少内容只能留坑了；

启发式的策略主要有以下几种：

  1. 误差估算与细分
    在自适应梯形积分中，可以利用Richardson外推或其他误差估计技术估算当前的误差。如果对某个区间的误差估计超出了预期阈值，只有这些区间会被进一步细分，从而减少不必要的计算。
  2. 局部适应性
    对于有急剧变化或高度非线性的函数，使用相同数量的区间可能会导致某些部分的误差大于其他部分。在这种情况下，可以设计算法以更高的密度对这些特定区域进行采样，同时对函数较为平缓的部分进行较少的采样。同时，也应该自适应的调整步长，在函数变化剧烈的地方采用较小的步长，而在平缓区域使用较大的步长。
  3. 多步骤递增
    不是一开始就使用非常高的划分数，可以从一个较小的区间数开始，并逐步增加，直到满足精确度要求，以此来优化计算时间。
  4. 停止准则的动态调整
    动态调整误差阈值，可以在算法运行过程中根据已经计算的结果和剩余区间来适应性地调整。

Content

梯形积分

在matlab中有一个现成的函数可供调用：
Q=trapz(x,y)表示y在x细分下的梯形积分，注意Q求出来并不是真实精确的积分结果；

梯形公式将图像线性近似$f(x)=f(a)+\frac{f(b)-f(a)}{b-a}(x-a)$，并用线下方的梯形面积替代积分面积：
$$
\int_{a}^{b}f(x)dx\approx(b-a)\frac{f(a)+f(b)}{2}
$$

其中如果将区间细化为$a=x_0<x1<…<x_{N-1}<x_N=b$,那么计算结果为
$$
\int_{a}^{b}f(x)dx=\frac{b-a}{2N}(f(x_0)+\sum_{i=0}^{N-1}f(x_i)+f(x_N))
$$

$$
R_T(f)=-\frac{(b-a)^2}{12}f’’(\xi)
$$

用matlab实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
x=linspace(0,pi,3);
y=sin(x);
my_trapz_res=my_trapz(x,y)
trapz_res=trapz(x,y)

function res=my_trapz0(f,a,b)
res=(b-a)*(f(b)+f(a))/2;
end

function res=my_trapz(x,y)
%x,y传入的应该是一对长度相等的向量
n=length(x);
res=0;
for i=2:n
res=res+(x(i)-x(i-1))*(y(i)+y(i-1))/2;
end
end

从结果上看到，和官方提供的函数并无多大区别：

1
2
3
4
5
6
7
8
my_trapz_res =

    1.5708


trapz_res =

    1.5708

在自适应梯形积分中，可以利用Richardson外推或其他误差估计技术估算当前的误差。如果对某个区间的误差估计超出了预期阈值，只有这些区间会被进一步细分，从而减少不必要的计算。（留坑）

自适应Simpson积分公式

在Simpson积分中，将图像近似为一个二次函数$f(x)=Ax^2+Bx+C$,那么积分可以做如下近似：
$$
\int_{a}^{b}f(x)\approx\frac{A}{3}(b^3-a^3)+\frac{B}2(b^2-a^2)+C(b-a)
$$

$$
=(b-a)[\frac{A}{3}(b^2+a^2+ab)+\frac{B}{2}(b+a)+C]=(b-a)\frac{f(a)+f(b)+4f(\frac{a+b}{2})}6
$$

我们也采用变步长的启发式策略：

  • 对于一个极小的区间认为他的Simpson值加上误差容忍可以近似为积分值；
  • 二分区间分别计算左边区间的Simpson值和右边区间的Simpson值
  • 如果两值加起来与整个区间Simpson值差距不超过15倍得误差容忍度，就认为这个区间足够小，返回这个小区间得积分值；
  • 否则需要继续递归，返回区间的左边积分值加上右边积分值；

matlab程序如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
quad(@(x)sin(x),0,pi)
my_quad(@(x) sin(x),0,pi)


function res=my_simpson(f,a,b)
res=(b-a)*(f(a)+f(b)+4*f((a+b)/2))/6;
end

function res=my_quad(f,a,b,tol)
if nargin==3
    tol=10^(-6);
end
mid=(a+b)/2;
quad_all=my_simpson(f,a,b);
quad_l=my_simpson(f,a,mid);
quad_r=my_simpson(f,mid,b);
if(abs(quad_l+quad_r-quad_all)<15*tol) 
    res=quad_l+quad_r+(quad_l+quad_r-quad_all)/15; 
else 
    res=my_quad(f,a,mid,tol)+my_quad(f,mid,b,tol);
end
end

可以看到和官方提供的quad函数没有太大差距：

1
2
3
4
5
6
7
8
ans =

    2.0000


ans =

    2.0000

洛谷模板题

请注意，自适应simpson公式数值积分适用于精度要求低，被积函数平滑性较差的数值积分；

高精度Lobatto积分法

留坑，暂时找不到资料

Remark

这里大多数数值积分的介绍都缺少一点理论知识，等哪天补一补数值计算再回来填坑了；

其中的启发式算法还是很不错的；