シンプソン則とは

数値積分の解法の一つ
関数$f(x)$において、微小区間の関数値を二次方程式で近似
微小区間$[x_0,x_2]$の端点とその中点$x_1$を用いて、二次方程式を求める

算法

三点の$y$座標をそれぞれ$f(x_0)=y_0,f(x_1)=y_1,f(x_2)=y_2$とする。
近似に用いる二次方程式$y=ax^2+bx+c$は、$(x_0,y_0),(x_1,y_1),(x_2,y_2)$を通る、また、$x_2-x_1=h,x_1-x_0=h$とすると、
$$ a=\frac{y_0-2y_1+y_2}{2h^2}$$$$b=\frac{-(x_1+x_2)y_0+2(x_2+x_0)y_1-(x_0+x_1)y_2}{2h^2}$$$$c=\frac{x_1x_2y_0-2x_2x_0y_1+x_0x_1y_2}{2h^2} $$
この二次方程式$y=ax^2+bx+c$を$x_0$から$x_2$まで積分すると、
$$ \int_{x_0}^{x_2}(ax^2+bx+c)dx=\frac{a}{3}(x_2^3-x_0^3)+\frac{b}{2}(x_2^2-x_0^2)+c(x_2-x_0) $$$$=\frac{h}{3}(y_0+4y_1+y_2)$$

区間$[a,b]$が微小距離$h$で$2n$等分されているとし、各$x_i\ (i=0,1,\cdots,2n)$の$y$座標を$f(x_i)=y_i$とすると、区間$[a,b]$内の積分は次のようになる。
$$\int_a^bf(x)dx=\frac{h}{3}{(y_0+4y_1+y_2)+(y_2+4y_3+y_4)+\cdots+(y_{2n-2}+4y_{2n-1}+y_{2n})}$$$$=\frac{h}{3}(y_0+4y_1+2y_2+4y_3+2y_4+\cdots+2y_{2n-2}+4y_{2n-1}+y_{2n})$$
ただし、$h=(b-a)/(2n)$

サンプルコード

$f(x)=\sqrt{1-x^2}$において、区間$[0,1]$の定積分の値を求めるプログラム。
分割数は4。
解析解は$\pi/4$です。

simpsons_rule.c

#include<stdio.h>
#include<math.h>

double f (double x) {
  return sqrt(1-x*x);
}

/* シンプソン則（区間[a,b]をn分割） */
double simpsons_rule (double a, double b, int n) {
  double h;
  int i;
  double value=0;

  h = (b - a) / (2*n);
  for (i = 0; i <= 2*n; i++) {
    if (i == 0 || i == 2*n) value += f(a + i*h);   // 0か2nのとき
    else if (i % 2 == 1) value += 4 * f(a + i*h);   // 奇数
    else value += 2 * f(a + i*h);   // 偶数
  }
  value = value*h/3;
  return value;
}

int main (void) {
  printf("Analytical solution: %f\n", M_PI/4);
  printf("Numerical solution : %f\n", simpsons_rule(0, 1, 4));
  return 0;
}

実行結果

Analytical solution: 0.785398
Numerical solution : 0.780297

特徴

一区間の誤差評価：$O(h^5/90)$
全区間合計の誤差評価：$O(h^4)$