ニュートン・コーツの公式とシンプソンの公式の導出と実装[C言語]

ルンゲクッタ法はどうやって導かれたのか謎だったのですが，『情報処理入門コース　数値計算』によれば，シンプソンの公式から来ているとのこと．
なのでシンプソンの公式がどういう経緯でできたのかメモしていきます．

※数式がうまく表示されない時はPCのブラウザから見ることを推奨します

1.ニュートン・コーツの積分公式

まずラグランジュ補間の公式が与えられているものとして使います．f(x)が近似したい式，p(x)が近似式．x_0,...,x_nが関数f(x)上のわかっている点． $$p(x)=\sum_{i=0}^{n}f(x_i)N(i,x) \\ N(i,x) = \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j}$$ ラグランジュ補間の出処はよくわかんないけど関数N(x)がクロネッカーのデルタ関数のようになっていることから意味をお察しください．とりあえずラグランジュ補間については前書いたのでそっちも見てくださいな．

次に，このラグランジュ補間で求めた近似式を用いて関数f(x)を区間[a,b]で積分する事を考えます． $$a = x_0,x_1,x_2,...,x_n = b$$ という様に区間をn分割し，上記のラグランジュ補間を適用すると次のようになります． $$\begin{eqnarray*} S &=& \int_{a}^{b}f(x)dx \\ & \approx & \int_{a}^{b}p(x)dx \\ &=& \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ \end{eqnarray*}$$ ここで少しわかりづらいのですが，x_iは定数なのでf(x_i)はxの関数では無いことから， $$\begin{eqnarray*} S & \approx & \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \end{eqnarray*}$$ という式が得られます．ここで関数N(i,x)はx_0,...,x_nが[a,b]をn等分した値，すなわち $$x_k = a + k\Delta x \\ \Delta x = \frac{b-a}{n}$$ と表せることを用いて，次の変数変換をすることを考えます． $$t = \frac{x-a}{\Delta x}\\ dt = \Delta x \, dx$$ ここで，N(i,x)は上記のtを代入して次のように変形したものをM(i,t)とする． $$\begin{eqnarray*} N(i,x) &=& \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} \\ &=& \prod_{j=0,j \neq i}^{n} \frac{(a+ t \Delta x) - (a + j \Delta x)}{(a+i \Delta x) - (a + j \Delta x ) } \\ &=& \prod_{j=0,j \neq i}^{n} \frac{t-j}{i-j} \\ &=& M(i,t) \end{eqnarray*}$$ このN(i,x)を変形した関数M(i,t)を用いると先ほどの積分はこうなる． $$\begin{eqnarray*} S & \approx & \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \Delta x \int_{0}^{n}M(i,t) dt \end{eqnarray*}$$ M(i,t)は変数xや，区間[a,b]などに依存しないから， $$C_i = \int_{0}^{n}M(i,t) dt$$ と置いて，先に手計算で計算しておけば使い回しが効いて，最終的に次のようになります． $$S \approx \sum_{i=0}^{n}f(x_i) C_i \Delta x$$ この公式をニュートン・コーツの積分公式と呼ぶとのこと．

2.シンプソンの公式

うえで求めたニュートン・コーツの積分公式にn=2を代入したものが世間一般にシンプソンの公式と呼ばれているものとのこと．C_iの値は上式に突っ込んで $$\begin{eqnarray*} C_0 &=& \int_0^2 \frac{t-1}{-1} \frac{t-2}{-2} dt &=& \frac{1}{3} \\ C_1 &=& \int_0^2 \frac{t}{1} \frac{t-2}{-1} dt &=& \frac{4}{3} \\ C_2 &=& \int_0^2 \frac{t}{2} \frac{t-1}{1} dt &=& \frac{1}{3} \\ \end{eqnarray*}$$ となる．この公式のすごいところはこの係数をあらかじめ求めておけば面倒な計算を端折ることができるところだと思います．
この値をニュートン・コーツの積分公式に代入すると， $$\begin{eqnarray*} \int_{a}^{b} f(x) dx &=& \sum_{i=0}^{2}f(x_i) C_i \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) {\frac{b-a}{2}}\\ &=& \frac{1}{6}(b-a)(f(x_0)+4f(x_1) +f(x_2)) \end{eqnarray*}$$ となります．

シンプソンの公式を図で表すとこんな感じ．青色の線が二次式で近似した曲線．その曲線の積分を求めてる．

ちなみにn=3の場合などはニュートン・コーツの公式のWikipediaに載ってたりするのでそちらを参照すると良いかと思います．

3.プログラム

抽象的な話はこの辺にして，実際にプログラムを書いてみます．試しに$f(x) = \frac{4}{x^2+1}$を[0,1]区間で積分してみます．

ちなみに$\int_{0}^{1}f(x)dx$の厳密解はπです．

まず $\Delta x = \frac{1}{n}$ とし，$x_k = k \Delta x$ と置くと， $$\begin{eqnarray*} \int_{0}^{1}f(x)dx = \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}}f(x) dx \end{eqnarray*}$$ という微小区間 $[x_k,x_{k+1}]$ の積分となります．この微小な区間に対してシンプソンの公式を適用すると次のようになります． $$\int_{x_k}^{x_{k+1}}f(x) dx \approx \frac{1}{6}(f(x_k)+4f(x_k+\frac{\Delta x}{2}) + f(x_{k+1}) )\Delta x$$ 実際はf(x_k)の足し算がかぶっているのでもう少しスマートにできるのですが，今回は愚直にそのままこれを実装します．

	/* シンプソンの公式による1/(x^2+1)の計算 */
	#include <stdio.h>
	// 被積分関数
	double f(double x){
	return 4.0/(x*x+1.0);
	}
	int main(void){
	int k; // ループカウンタ
	int N = 1000; // 分割数
	double dx = 1.0/N; // 刻み幅Δx
	double S = 0.0; // 和を格納する変数
	double a; // x_kの値を格納する
	double b; // x_{k+1}の値を格納する
	double m; // a,bの中間の値を格納する
	for( k = 0 ; k < N ; k++ ){
	// 近似を行う3点を計算
	a = k * dx;
	m = a + dx / 2.0;
	b = (k+1) * dx;
	// f_k + g(x_k,f(x_k))をfに代入
	S += (1.0/6.0)*(f(a)+4.0*f(m)+f(b))*dx;
	// 値の表示
	printf("%lf,%e\n",b,S);
	}
	return 0;
	}

view raw SimpsonsRule.c hosted with ❤ by GitHub

実行結果

結果が3.141593なので円周率が概ねちゃんと求まっているので大丈夫そうです．