こんにちは、Rafkaです。
前回から間が書いてしまいましたが、今回は常微分方程式の数値解法の１つで、最も単純なオイラー法の導出と実装を行っていきます。
前回 → 数値シミュレーション －微分方程式－
実装をすると言っても、前回の時に実装は済ませてしまっているので、導出した後に改めてソースを見直して、アルゴリズムとの照らし合わせをする感じになりそうです。
数値解法の共通の考え方
第１回の時に、微分方程式の数値解法は、未知関数 \(\boldsymbol{x}(t)\) の変化量を定義している右辺の式 \(\boldsymbol{f}\bigl( t, \boldsymbol{x}(t) \bigr)\) の値を数値的に積分することと書きました。もっと単純に書くと、右辺に定義される式から計算できる微小時間 \(\Delta t\) あたりの変化量 \(\Delta \boldsymbol{x}(t)\) を計算し、その値を用いて \(\boldsymbol{x}(t)\) の値を更新していきます。
$$\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \Delta \boldsymbol{x}(t)$$
微分方程式は限りなく小さい微小時間についての定義であることや、上に示した更新式からも分かる通り、\(\Delta t\) をなるべく小さくすることで、より精緻なシミュレーションを行うことができますが、そもそもコンピュータを用いて計算を行うので、取れる小さい値には限度があることや、細かい \(\Delta t\) を用いると計算ステップ数が増大してコストが大きくなってしまう問題があります。ですので計算誤差と計算コストのバランスを考えて \(\Delta t\)の値を設定する必要があるのですが、計算誤差は \(\Delta \boldsymbol{x}(t)\) の計算方法を工夫することである程度改善することができます。その計算方法の種類に、今回紹介するオイラー法やルンゲクッタ法などが存在します。
オイラー法の導出
いくつかあるのですが、他の方法よりもいくらか直感的に理解できる方法から解説したいと思います。
● 直感的な導出
微分方程式の最初の式に手を加えていきます。
$$\frac{d\boldsymbol{x}(t)}{dt} = \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)$$
左辺の微分の定義は以下の式です。
$$\frac{d\boldsymbol{x}(t)}{dt} = \lim_{\Delta t \to 0} \frac{\boldsymbol{x}(t + \Delta t) \;-\; \boldsymbol{x}(t)}{\Delta t}$$
\( \lim_{\Delta \to 0} \) を取り除いてしまって、前の式に代入してしまいましょう。すると、以下の式になります。
$$\frac{\boldsymbol{x}(t + \Delta t) \;-\; \boldsymbol{x}(t)}{\Delta t} = \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)$$
ここから式変形して先に示した更新式に近い形にしていきます。
$$\begin{align} \boldsymbol{x}(t + \Delta t) \;-\; \boldsymbol{x}(t) &= \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\Delta t\\ \boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\Delta t \end{align}$$
これで更新式の変化量を \(\Delta \boldsymbol{x}(t) = \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\Delta t\) とすれば良いことがわかりました。
以上で更新式の導出自体はできるのですが、誤差の議論とかは全くできないので、一応テイラー展開を用いた導出も書いておきます。
● テイラー展開を用いた導出
次の時間ステップの値 \(\boldsymbol{x}(t + \Delta t)\) を \(t\) の周りでテイラー展開します。
$$\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \frac{d\boldsymbol{x}(t)}{dt}\Delta t + O(\Delta t^2)$$
最初に示した微分方程式の定義より、\(d\boldsymbol{x}(t) / dt = \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\)ですのでこれを代入します。
$$\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\Delta t + O(\Delta t^2)$$
これで先程示した式と同じ式が導出できました。\(O(\Delta t^2)\) から、オイラー法が計算精度１次の手法であることも分かります。
実装
以上でオイラー法のアルゴリズムの導出が終わったので、実装の方に移ろうと思います。
ちなみに今回示したオイラー法は、厳密には陽的オイラー法と呼ばれるものです。陰的オイラー法と呼ばれるものもあるので、興味が湧いた方はぜひ調べてみてください。
さて、前回の記事で示したように、オイラー法は以下のように実装しました。
２行目のDeltaが \(\Delta t\) に相当します。７～９行目でシミュレーションの開始の時刻と終了時刻、微分方程式の初期値の読み込みを行っています。１２～１６行目のwhileが、実際にオイラー法で計算を行っています。１３行目でxの更新に使っているEulerUpdateメソッドが\(\Delta \boldsymbol{x}(t)\) に相当し、このメソッドの実装は１８行目にあり、導出した \(\boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\Delta t\) を返すようになっています。
まとめ
今回は微分方程式の数値解法の１つとして、オイラー法の紹介を行いました。
オイラー法での更新式は以下になります。
$$\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \boldsymbol{f} \bigl( t, \boldsymbol{x}(t) \bigr)\Delta t$$
次回、次次回ではオイラー法よりも誤差を小さくできるホイン法とルンゲクッタ法の紹介を行います。全て紹介し終わった後、斜方投射を題材に、各手法でどの程度誤差が小さくなるのかを示します。
次回 → 数値シミュレーション －ホイン法－