Merge pull request #2526 from zyj2297349886/integral

简单来说，函数 $f(x)$ 在区间 $[l,r]$ 上的定积分 $\int_{l}^{r}f(x)\mathrm{d}x$ 指的是 $f(x)$ 在区间 $[l,r]$ 中与 $x$ 轴围成的区域的面积（其中 $x$ 轴上方的部分为正值， $x$ 轴下方的部分为负值）。

很多情况下，我们需要高效，准确地求出一个积分的近似值。下面介绍的 **自适应辛普森法** ，就是这样一种求数值积分的方法。

很多情况下，我们需要高效，准确地求出一个积分的近似值。下面介绍的 **辛普森法** ，就是这样一种求数值积分的方法。

## 自适应辛普森法

## 辛普森法

这个方法的思想是将被积区间分为若干小段，每段套用二次函数的积分公式进行计算。

    \end{aligned}

    $$

根据这个辛普森公式，我们先介绍一种普通的辛普森积分法。

### 普通辛普森法

1743 年，这种方法发表于托马斯·辛普森的一篇论文中。

#### 描述

给定一个自然数 $n$ ，将区间 $[l, r]$ 分成 $2n$ 个等长的区间 $x$ 。

 $x_i = l + i h, ~~ i = 0 \ldots 2n,$  $h = \frac {r-l} {2n}.$ 

我们就可以计算每个小区间 $[x_ {2i-2}, x_ {2i}]$ , $i = 1\ldots n$ 的积分值，将所有区间的积分值相加即为总积分。

对于 $[x_ {2i-2}, x_ {2i}]$ , $i = 1\ldots n$ 的一个区间，选其中的三个点 $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$ 就可以构成一条抛物线从而得到一个函数 $P(x)$ ，这个函数存在且唯一。计算原函数在该区间的积分值就变成了计算新的二次函数 $P(x)$ 在该段区间的积分值。这样我们就可以利用辛普森公式来近似计算它。

 $\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3}$ 

将其分段求和即可得到如下结论：

 $\int_l ^ r f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3}$ 

#### 误差

我们直接给出结论，普通辛普森法的误差为：

$$

-\tfrac{1}{90} \left(\tfrac{r-l}{2}\right)^5 f^{(4)}(\xi)

$$

其中 $\xi$ 是位于区间 $[l,r]$ 的某个值。

#### 实现

```cpp

const int N = 1000 * 1000;

double simpson_integration(double a, double b) {

  double h = (b - a) / N;

  double s = f(a) + f(b);

  for (int i = 1; i <= N - 1; ++i) {

    double x = a + h * i;

    s += f(x) * ((i & 1) ? 4 : 2);

  }

  s *= h / 3;

  return s;

}

```

### 自适应辛普森法

普通的方法为保证精度在时间方面无疑会受到 $n$ 的限制，我们应该找一种更加合适的方法。

现在唯一的问题就是如何进行分段。如果段数少了计算误差就大，段数多了时间效率又会低。我们需要找到一个准确度和效率的平衡点。

我们这样考虑：假如有一段图像已经很接近二次函数的话，直接带入公式求积分，得到的值精度就很高了，不需要再继续分割这一段了。