辛普森积分

  |  

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

【对数据分析、人工智能、金融科技、风控服务感兴趣的同学,欢迎关注我哈,阅读更多原创文章】
我的网站:潮汐朝夕的生活实验室
我的公众号:潮汐朝夕
我的知乎:潮汐朝夕
我的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 梯形法则

对于区间 [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)dx \approx \int_{a}^{b}g(x)dx = 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;
}

自适应辛普森积分

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

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

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)dx$,其中 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)dx$ 的近似。

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

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

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


Share