跳转至

數值積分

定積分的定義

簡單來説，函數 \(f(x)\) 在區間 \([l,r]\) 上的定積分 \(\int_{l}^{r}f(x)\mathrm{d}x\) 指的是 \(f(x)\) 在區間 \([l,r]\) 中與 \(x\) 軸圍成的區域的面積（其中 \(x\) 軸上方的部分為正值，\(x\) 軸下方的部分為負值）。

很多情況下，我們需要高效，準確地求出一個積分的近似值。下面介紹的 辛普森法，就是這樣一種求數值積分的方法。

辛普森法

這個方法的思想是將被積區間分為若干小段，每段套用二次函數的積分公式進行計算。

二次函數積分公式（辛普森公式）

對於一個二次函數 \(f(x)=ax^2+bx+c\)，有：

\[ \int_l^r f(x) {\mathrm d}x = \frac{(r-l)(f(l)+f(r)+4 f(\frac{l+r}{2}))}{6} \]

推導過程： 對於一個二次函數 \(f(x)=ax^2+bx+c\)； 求積分可得 \(F(x)=\int_0^x f(x) {\mathrm d}x = \frac{a}{3}x^3+\frac{b}{2}x^2+cx+D\) 在這裏 D 是一個常數，那麼

\[ \begin{aligned} \int_l^r f(x) {\mathrm d}x &= F(r)-F(l) \\ &= \frac{a}{3}(r^3-l^3)+\frac{b}{2}(r^2-l^2)+c(r-l) \\ &=(r-l)(\frac{a}{3}(l^2+r^2+lr)+\frac{b}{2}(l+r)+c) \\ &=\frac{r-l}{6}(2al^2+2ar^2+2alr+3bl+3br+6c)\\ &=\frac{r-l}{6}((al^2+bl+c)+(ar^2+br+c)+4(a(\frac{l+r}{2})^2+b(\frac{l+r}{2})+c)) \\ &=\frac{r-l}{6}(f(l)+f(r)+4f(\frac{l+r}{2})) \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]\) 的某個值。

實現

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
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;
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
N = 1000 * 1000

def simpson_integration(a, b):
    h = (b - a) / N
    s = f(a) + f(b)
    for i in range(1, N):
        x = a + h * i
        if i & 1:
            s = s + f(x) * 4
        else:
            s = s + f(x) * 2
    s = s * (h / 3)
    return s

自適應辛普森法

普通的方法為保證精度在時間方面無疑會受到 \(n\) 的限制，我們應該找一種更加合適的方法。

現在唯一的問題就是如何進行分段。如果段數少了計算誤差就大，段數多了時間效率又會低。我們需要找到一個準確度和效率的平衡點。

我們這樣考慮：假如有一段圖像已經很接近二次函數的話，直接帶入公式求積分，得到的值精度就很高了，不需要再繼續分割這一段了。

於是我們有了這樣一種分割方法：每次判斷當前段和二次函數的相似程度，如果足夠相似的話就直接代入公式計算，否則將當前段分割成左右兩段遞歸求解。

現在就剩下一個問題了：如果判斷每一段和二次函數是否相似？

我們把當前段直接代入公式求積分，再將當前段從中點分割成兩段，把這兩段再直接代入公式求積分。如果當前段的積分和分割成兩段後的積分之和相差很小的話，就可以認為當前段和二次函數很相似了，不用再遞歸分割了。

上面就是自適應辛普森法的思想。在分治判斷的時候，除了判斷精度是否正確，一般還要強制執行最少的迭代次數。

參考代碼如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
double simpson(double l, double r) {
  double mid = (l + r) / 2;
  return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6;  // 辛普森公式
}

double asr(double l, double r, double eps, double ans, int step) {
  double mid = (l + r) / 2;
  double fl = simpson(l, mid), fr = simpson(mid, r);
  if (abs(fl + fr - ans) <= 15 * eps && step < 0)
    return fl + fr + (fl + fr - ans) / 15;  // 足夠相似的話就直接返回
  return asr(l, mid, eps / 2, fl, step - 1) +
         asr(mid, r, eps / 2, fr, step - 1);  // 否則分割成兩段遞歸求解
}

double calc(double l, double r, double eps) {
  return asr(l, r, eps, simpson(l, r), 12);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
def simpson(l, r):
    mid = (l + r) / 2
    return (r - l) * (f(l) + 4 * f(mid) + f(r)) / 6 # 辛普森公式
def asr(l, r, eps, ans, step):
    mid = (l + r) / 2
    fl = simpson(l, mid); fr = simpson(mid, r)
    if abs(fl + fr - ans) <= 15 * eps and step < 0:
        return fl + fr + (fl + fr - ans) / 15 # 足夠相似的話就直接返回
    return asr(l, mid, eps / 2, fl, step - 1) + \
          asr(mid, r, eps / 2, fr, step - 1) # 否則分割成兩段遞歸求解
def calc(l, r, eps):
    return asr(l, r, eps, simpson(l, r), 12)

習題

參考資料

https://doi.org/10.1145/321526.321537：該文章討論了自適應 Simpson 法的改進方案，其中詳細論述了上文代碼中的常數 15 的由來與優勢。

本页面最近更新：，更新历史
发现错误？想一起完善？ 在 GitHub 上编辑此页！
本页面贡献者：H-J-Granger, Chrogeek, countercurrent-time, Enter-tainer, Great-designer, iamtwz, Ir1d, ksyx, mao1t, Menci, NachtgeistW, Nanarikom, ShaoChenHeng, StudyingFather, SukkaW, Tiphereth-A, zyj-111
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用