シンプソン則とは
- 数値積分の解法の一つ
- 関数$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$です。
#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)$