為防止連結失效，轉載c++實現數值積分的方法（辛普森法）。

數值積分

普通辛普森法

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

// C++ Version
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;
}
 

# Python Version
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的限制，我們應該找一種更加合適的方法。

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

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

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

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

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

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

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

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

文章轉載自https://oi-wiki.org/math/integral/，或點選這裡跳轉。