数计算数学

辛普森积分

字数统计: 1.6k字   |   阅读时长: 7分
,

摘要: 辛普森积分的原理、代码、例题

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：潮汐朝夕
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

$0 背景

梯形法则和辛普森法则都是数值积分的方法，用于计算定积分，它们都是牛顿-柯特斯公式的特例。

用牛顿-柯特斯公式求定积分对于次数较高的多项式有很大的误差，也就是龙格现象。次数较高的时候应该用高斯积分法。

$1 定积分与数值积分

数值积分

其中 $f(x)$ 是一个函数并且不容易求原函数。计算这个定积分可以把区间 $[a, b]$ 分成很多小区间，然后把小区间的积分再求和。

对于小区间 $[a, b]$，选择一个容易计算原函数的函数 $g(x)$ 拟合 $f(x)$ 并求 $g(x)$ 在 $[a, b]$ 上的定积分，将所有小区间上 $g(x)$ 的积分值加起来得到所求积分值的近似值。

数值积分的作用：

  1. 不容易求原函数的场景
  2. 只能知道积分函数在某些特定点取值的场景
  3. 积分区域是曲面等高维流形的时候，牛顿莱布尼茨公式不再适用

代数精度

如果数值求积公式对于任何不高于m次的代数多项式都准确成立，而对m+1次代数多项式不准确成立，则称该求积公式具有m次代数精确度，简称代数精度。

$2 梯形法则

对于区间 $$h[a, b] 用直线拟合 $f(x)$，然后求直线的积分。

将被积函数近似为直线，被积部分近似为梯形。用 $f(a)$ 和 $f(b)$ 即可拟合出直线并求出积分值。

分成多个小区间

将被积区间分成多个小区间，分别求近似积分值再加起来：

可以写成：

其中 $k = 0,1,…,n,x_{k}=a+k\frac{b-a}{n}$。

$3 辛普森法则

对于区间 $[a, b]$ 用二次函数拟合 $f(x)$，然后求二次函数的积分：

将被积函数近似为二次函数。用 $f(a)$, $f(b)$ 和 $f((a+b)/2)$ 即可拟合出二次函数并求出积分值。

辛普森公式的推导

$f(x)$ 为被积函数，$g(x) = Ax^{2} + Bx + C$ 为拟合的函数。$g(x)$ 是比较好求出原函数的 $G(x)$ 的，因此 $\int_{a}^{b}f(x)\mathrm{d}x \approx \int_{a}^{b}g(x)\mathrm{d}x = G(b) - G(a)$。

分成多个小区间

将被积区间分成多个小区间，分别求近似积分值再加起来：

其中 $h = \frac{b-a}{n}$, $x_{k} = a + kh$。

代码 : 辛普森公式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
struct Function
{
    ...
    double operator()(double x) const
    {
        ...
    }
};

double simpson(Function& f, double a, double b)
{
    double mid = (a + b) / 2;
    return (f(a) + 4 * f(mid) + f(b)) * (b - a) / 6;
}

自适应辛普森积分

辛普森公式计算定积分需要注意精度问题。区间分少了精度会不够，区间分多了时间成本增加。自适应辛普森积分的目的就是自动控制分割的小区间的大小。

算法: 二分地递归，当满足精度需要后终止递归。

代码（C++，模板）

1
2
3
4
5
6
7
8
9
10
11
12
const double EPS = 1e-10;

double adaptive_simpson(Function& f, double a, double b)
{
    double mid = (a + b) / 2;
    double l_ans = simpson(f, a, mid);
    double r_ans = simpson(f, mid, b);
    double ans = simpson(f, a, b);
    if(fabs(l_ans + r_ans - ans) < EPS)
        return ans;
    return adaptive_simpson(f, a, mid) + adaptive_simpson(f, mid, b);
}

题目 P4525 自适应辛普森法1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

const double EPS = 1e-9;

struct Function
{
    double a, b, c, d;
    Function(double a, double b, double c, double d):a(a),b(b),c(c),d(d){}
    double operator()(double x) const
    {
        return (c * x + d) / (a * x + b);
    }
};

double simpson(Function &f, double a, double b)
{
    // f(x) 在 [a, b] 上的积分
    double mid = (a + b) / 2;
    return (f(a) + 4 * f(mid) + f(b)) * (b - a) / 6;
}

double adaptive_simpson(Function& f, double a, double b)
{
    double mid = (a + b) / 2;
    double l_ans = simpson(f, a, mid);
    double r_ans = simpson(f, mid, b);
    double ans = simpson(f, a, b);
    if(fabs(l_ans + r_ans - ans) < EPS)
        return ans;
    return adaptive_simpson(f, a, mid) + adaptive_simpson(f, mid, b);
}

int main()
{
    double a, b, c, d, L, R;
    cin >> a >> b >> c >> d >> L >> R;
    Function f(a, b, c, d);
    double ans = adaptive_simpson(f, L, R);
    cout << std::fixed << std::setprecision(6);
    cout << ans + 1e-9 << endl;
}

$4 牛顿-柯特斯公式

要求积分 $\int_{a}^{b}f(x)\mathrm{d}x$，其中 $f(x)$ 不好求原函数。用一个 n 次多项式 $g(x)$ 拟合 $f(x)$ 将原问题转化为求 n 次多项式 $g(x)$ 的定积分，而 n 多项式是好求原函数的。

整个过程可以理解为先用 n + 1 个点 $f(x_{0}), f(x_{0}), …, f(x_{0})$ 进行插值，求得对应 $f(x)$ 的拉格朗日多项式，然后对该 n 次多项式求积分，该积分为 $\int_{a}^{b}f(x)\mathrm{d}x$ 的近似。

其中 k=0,1,…,n, $w_{i}$ 是常数, $x_{k} = a + k\frac{b-a}{n}$。

梯形法则和辛普森法则是 n = 1,2 的情况。

牛顿-柯特斯公式对于次数较高的多项式误差很大，也就是龙格现象，此时可以考虑龙贝格求积公式，高斯积分法，蒙特卡洛法。

Share