数値シミュレーション －オイラー法－

こんにちは、Rafka です。

前回から間が書いてしまいましたが、今回は常微分方程式の数値解法の１つで、最も単純なオイラー法の導出と実装を行っていきます。

実装をすると言っても、前回の時に実装は済ませてしまっているので、導出した後に改めてソースを見直して、アルゴリズムとの照らし合わせをする感じになりそうです。

第 1 回の時に、 微分方程式の数値解法は、 未知関数 $\boldsymbol{x}(t)$ の変化量を定義している 右辺の式 $\boldsymbol{f}\bigl( \boldsymbol{x}(t), t \bigr)$ の値を 数値的に積分することと書きました。 もっと単純に書くと、 右辺に定義される式から計算できる微小時間 $\Delta t$ あたりの 変化量 $\Delta \boldsymbol{x}(t)$ を計算し、 その値を用いて $\boldsymbol{x}(t)$ の値を更新していきます。

微分方程式は限りなく小さい微小時間についての定義であることや、 上に示した更新式からも分かる通り、 $\Delta t$ をなるべく小さくすることで、 より精緻なシミュレーションを行うことができますが、 そもそもコンピュータを用いて計算を行うので、 取れる小さい値には限度があることや、 細かい $\Delta t$ を用いると計算ステップ数が増大してコストが大きくなってしまう問題があります。 ですので、計算誤差と計算コストのバランスを考えて $\Delta t$ の値を設定する必要があるのですが、 計算誤差は $\Delta \boldsymbol{x}(t)$ の計算方法を工夫することである程度改善することができます。 その計算方法の種類に、今回紹介するオイラー法やルンゲクッタ法などが存在します。

いくつかあるのですが、他の方法よりもいくらか直感的に理解できる方法から解説したいと思います。

微分方程式の最初の式に手を加えていきます。

左辺の微分の定義は以下の式です。

$\lim_{\Delta \to 0}$ を取り除いてしまって、 前の式に代入してしまいましょう。 すると、以下の式になります。

ここから式変形して先に示した更新式に近い形にしていきます。

これで更新式の 変化量を $\Delta \boldsymbol{x}(t) = \boldsymbol{f} \bigl( \boldsymbol{x}(t), t \bigr)\Delta t$ と すれば良いことがわかりました。

以上で更新式の導出自体はできるのですが、誤差の議論とかは全くできないので、一応テイラー展開を用いた導出も書いておきます。

次の時間ステップの値 $\boldsymbol{x}(t + \Delta t)$ を $t$ の周りでテイラー展開します。

最初に示した微分方程式の定義より、 $d\boldsymbol{x}(t) / dt = \boldsymbol{f} \bigl( \boldsymbol{x}(t), t \bigr)$ ですのでこれを代入します。

これで先程示した式と同じ式が導出できました。 $O(\Delta t^2)$ から、オイラー法の式の範囲では 2 次以上の精度の誤差が無視されていることが分かり、 オイラー法が計算精度 1 次の手法であることも分かります。

計算精度が 1 次ということは、刻み幅 $\Delta t$ を小さくしていくと、誤差がおおむね $\Delta t$ に比例して小さくなるということです。 簡単な微分方程式 $dx/dt = x$（初期値 $x(0) = 1$、厳密解 $x(t) = e^t$）で確かめてみましょう。

$\Delta t$ を半分にするたびに、$t = 2$ での誤差は 3.389 → 2.327 → 1.429 → 0.806 → 0.430 と小さくなり、 1 回あたりの縮み方は半分に近づいていきます。

近似解がいつも厳密解より下にあるのは、各ステップの始点での傾きだけを使って進めているからです。 $e^t$ のように傾きが増えていく関数では、始点の傾きが区間の中で一番小さいので、近似解は厳密解に追いつけません。

以上でオイラー法のアルゴリズムの導出が終わったので、実装の方に移ろうと思います。

ちなみに今回示したオイラー法は、厳密には陽的オイラー法と呼ばれるものです。陰的オイラー法と呼ばれるものもあるので、興味が湧いた方はぜひ調べてみてください。

さて、前回の記事で示したように、オイラー法は以下のように実装しました。

void Euler(ModelFunction f, Span<double> x, double t, double dt) {
  var n = x.Length;
  Span<double> dx = stackalloc double[n];
 
  f(x, t, dx);
  ScaleAdd(x, dx, dt, x);
}

5 行目の f(x, t, dx) が $\boldsymbol{f} \bigl( \boldsymbol{x}(t), t \bigr)$ に相当し、 6 行目の処理で、その値に $\Delta t$ を掛けて、 前の時間ステップの値 x に足し込むことで、 $\boldsymbol{x}(t)$ の更新を行っています。

今回は微分方程式の数値解法の１つとして、オイラー法の紹介を行いました。

次回、次次回ではオイラー法よりも誤差を小さくできるホイン法とルンゲクッタ法の紹介を行います。 全て紹介し終わった後、斜方投射を題材に、各手法でどの程度誤差が小さくなるのかを示します。