オイラー法と修正オイラー法の説明と実装例[C言語]

学科的に流体シミュレーションとか多分そのうちやると思うのですが，ルンゲクッタ法を理解していなかったのでとりあえず常微分方程式の数値計算の基本であるオイラー法と修正オイラー法について記事を書いてみます．

1.オイラー法

オイラー法は精度は悪いですが，非常に簡単で計算が高速です．研究などでは恐らく使われることはありませんが，ゲームプログラミング等をしているとシューティングの弾の運動などは大抵オイラー法が使われています．
今回は説明のための例として，次のような常微分方程式があったとします．今回は初期値f(x_0)がわかっているものとします． $$f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x))$$ ここで　$x_k = x_0 + k \Delta x$　と置くと　$f(x_{k})$　の値は次のようにして求められることがわかります． $$\begin{eqnarray*} f(x_{k}) &=& f(x_0) + \int_{x_0}^{x_k} g(x,f(x)) dx \\ &=& f(x_0) + \sum_{i=0}^{k-1} \int_{x_{i}}^{x_{i+1}} g(x,f(x)) dx \end{eqnarray*}$$ ここで　$f(x_{k+1})$　と　$f(x_{k})$　の漸化式を考えます． $$\begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*}$$ この $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分がキモです．この部分をどうするかで精度が大きく変わってきます．

オイラー法ではg(x,f(x))が区間$[x_k,x_{k+1}]$でほとんど変化しないものとして，次のように近似しています． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx g(x_k,f(x_k)) \Delta x$$ これを漸化式に代入したものがオイラー法です． $$f(x_{k+1}) \approx f(x_k) + g(x_k,f(x_k)) \Delta x$$

2.修正オイラー法

オイラー法は非常に簡単なのですが精度がよろしくないです．なのでもう少し精度を上げる必要があります．

修正オイラー法はその改良版の一つです．もっとも，これもあまり制度が良いというわけではないですが．

先ほど， $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分がキモと言いました．オイラー法では関数g(x,f(x))があまり変化しないものとして計算していました．イメージとしてはこんな感じです．

矩形近似とか長方形近似とか呼ばれているタイプの積分方法です．見るからに精度がよろしくなさそうなのがわかります．

そこで台形公式というものを使って修正オイラー法というものが提案されました．イメージとしてはこんな感じです．

名前の通り台形で近似するというものです．数式で表すと次のようになります．台形の面積は((高さ)×（上辺）+(下辺))÷2という小学生の時に習った公式を用いて，

$$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2}$$ という近似式で計算します．しかしこれを漸化式に代入すると $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2}$$ となり，$f(x_{k+1})$ を求めるために $f(x_{k+1})$ を使っているのでこれでは使えません．なので右辺の $f(x_{k+1})$ だけは，オイラー法で使った矩形近似を用いて次のように求めておきます． $$\widehat{f(x_{k+1})} = f(x_k) + g(x_k,f(x_k)) \Delta x$$ この　$\widehat{f(x_{k+1})}$　を用いて漸化式に代入したものが修正オイラー法です． $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2}$$

3.利用例

ここまで抽象的すぎたのでここで一度，次の1階微分方程式を数値的に解いてみたいと思います． $$\frac{df(x)}{dx} - 2f(x) = x$$ まぁ適当に作った式なので意味は無いです．ちなみに厳密解はCを定数として $$f(x) = -\frac{1}{4}(2x+1) + Ce^{2x}$$ と表されます．ちなみに定数Cは関数f(x)の初期値f(0) = aとした時，$C = a + \frac{1}{4}$になります．
今回は初期値 f(0) = a = 1 として数値的に求めてみます．

上記の微分方程式を変形してみます． $$\begin{eqnarray*} \frac{df(x)}{dx} &=& x + 2f(x) \\ &=& g(x,f(x)) \end{eqnarray*}$$ これでg(x,f(x))が求まりましたので，先ほど示した漸化式に代入してみれば，オイラー法では $$f(x_{k+1}) \approx f(x_k) + (x_k + 2f(x_k)) \Delta x$$ となります.

修正オイラー法では，まず　$\widehat{f(x_{k+1})}$　が $$\begin{eqnarray*} \widehat{f(x_{k+1})} &=& f(x_k) + g(x_k,f(x_k)) \Delta x \\ &=& f(x_k) + (x_k + 2f(x_k)) \Delta x \\ \end{eqnarray*}$$ となることから，次の漸化式に代入してやれば良いです． $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2} \\ $$ 全部数式を代入すれば,x_kとf(x_k)とΔxだけの式にできますが，数式が長くなるので省略．

もっとも，関数をプログラム内で定義してやればそんなに気にならないと思います．

4.精度の評価

ソースコードは長くなりますので先にオイラー法，修正オイラー法の精度を示しときます．

もっとも，評価と言ってもこの例での真値との差を見るだけですが．以下数値的に求めたものと，真値との差です．
刻み幅Δx=0.001，反復回数1000回です．euler's method 2の方が修正オイラー法です．適切な英語が見つからなかったので……．

1000回繰り返して差が10^-7オーダーですね．とりあえずそれらしい値は取れているのがわかります．

これを見る限り修正オイラー法はオイラー法と比べたら優秀そうではあります．
ただ傾きが指数関数的に増えていきそうなので後日ルンゲクッタ法を実装してそれと比較したいです．

5.ソースコード

C言語による実装．オイラー法はこんな感じです．

	/* Euler法による微分方程式 f' - 2f = x の計算 */
	#include <stdio.h>
	#include <math.h>
	double g(double x,double fx){
	return x + 2.0 * fx;
	}
	int main(void){
	int k; // ループカウンタ
	int N = 10000; // 繰り返し数
	double dx = 1e-3; // 刻み幅Δx
	double f = 1.0; // 関数f(x)の値(初期値)
	double x_0 = 0.0; // 初期値の時のxの値
	double x_k; // x_kの値
	for( k = 1 ; k <= N ; k++ ){
	x_k = x_0 + k * dx; // x_kの値を求める
	f = f + g(x_k,f) * dx; // f_k + g(x_k,f(x_k))をfに代入
	// 値の表示
	printf("%lf,%e\n",x_k,f);
	}
	return 0;
	}

view raw EulersMethod.c hosted with ❤ by GitHub

修正オイラー法はこんな感じです．

	/* 修正Euler法による微分方程式 f' - 2f = x の計算 */
	#include <stdio.h>
	#include <math.h>
	double g(double x,double fx){
	return x + 2.0 * fx;
	}
	int main(void){
	int k; // ループカウンタ
	int N = 10000; // 繰り返し数
	double dx = 1e-3; // 刻み幅Δx
	double f = 1.0; // 関数f(x)の値(初期値)
	double f_bar; // f(x_k+1)を矩形積分でとりあえず求めた値
	double x_0 = 0.0; // 初期値の時のxの値
	double x_k; // x_kの値
	for( k = 1 ; k <= N ; k++ ){
	x_k = x_0 + k * dx; // x_kの値を求める
	// まず矩形積分でf_k+1を見積もって
	f_bar = f + g(x_k,f) * dx;
	// 台形公式で積分する
	f = f + (g(x_k,f) + g(x_k + dx, f_bar)) * dx /2.0;
	// 値の表示
	printf("%lf,%e\n",x_k,f);
	}
	return 0;
	}

view raw EulersMethod2.c hosted with ❤ by GitHub

今日はここまで．