数値計算 理論

数値シミュレーション -ホイン法-

こんにちは、Rafkaです。

前回 → 数値シミュレーション -オイラー法-

今回は、前回紹介したオイラー法よりも誤差を抑えられる、ホイン法(2次のルンゲクッタ法とも呼ばれます)の導出と実装を行っていきます。

 

アルゴリズム


導出する前に、ホイン法のアルゴリズムを紹介します。ちなみにホインの綴りは少し特殊で、Heunと書きます。

ホイン法では以下の更新式で \(\boldsymbol{x}(t)\) の値を更新していきます。

$$\begin{align}\boldsymbol{k}_1 &= \boldsymbol{f} \bigl(t, \boldsymbol{x}(t)\bigr)\\ \boldsymbol{k}_2 &= \boldsymbol{f} \bigl(t + \Delta t, \boldsymbol{x}(t) + \boldsymbol{k}_1\Delta t\bigr)\\ \boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \frac{\Delta t}{2}(\boldsymbol{k}_1 + \boldsymbol{k}_2)\end{align}$$

上の式の感覚的な解釈としては、時間 \(t\)における傾き \(\boldsymbol{k}_1\) を計算し、次の時間ステップまでをその傾きの1次関数として近似した時の次のステップでの傾き \(\boldsymbol{k}_2\)を計算し、それらの平均を次の時間ステップへの値の更新に使用するといった感じです。

それでは導出の方に進みましょう。

 

導出


とりあえず \(\boldsymbol{x}(t + \Delta t)\) をテイラー展開します。前回のオイラー法よりも精度を上げたいので、今回は2次の項まで展開します。

$$\boldsymbol{x}(t + \Delta t) = \boldsymbol{x}(t) + \frac{d\boldsymbol{x}(t)}{dt}\Delta t + \frac{1}{2}\frac{d^2\boldsymbol{x}(t)}{dt^2}\Delta t^2 + O(\Delta t^3)$$

ここで厄介なのは2回微分の項です。2回微分の項は以下のように一度定義の式で \(\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\) に置き換えた後 \(\boldsymbol{f}\) について全微分を行うことで計算を進めることができます。

$$\begin{align}\frac{d^2\boldsymbol{x}(t)}{dt^2} &= \frac{\text{d}\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\text{d}t}\\ &= \frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial t} + \frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial\boldsymbol{x}}\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\end{align}$$

これをテイラー展開した式の2回微分の項に代入すると、以下のようになります。

$$\begin{align}\boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\Delta t\\ &+ \frac{1}{2}\Biggl(\frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial t} + \frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial\boldsymbol{x}}\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\Biggr)\Delta t^2 + O(\Delta t^3)\end{align}$$

偏微分の項が出てきてしまったので、未知定数 \(a, b\) を用いて表現された以下の式を組み合わせて上と同じ形に変形し、それぞれを比較して定数を決定して最終的な更新式の導出を行いたいと思います。

$$\begin{align}\boldsymbol{k}_1 &= \boldsymbol{f} \bigl(t, \boldsymbol{x}(t)\bigr)\\ \boldsymbol{k}_2 &= \boldsymbol{f} \bigl(t + \Delta t, \boldsymbol{x}(t) + \boldsymbol{k}_1\Delta t\bigr)\\ \boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \Delta t(a\boldsymbol{k}_1 + b\boldsymbol{k}_2)\end{align}$$

後で使うので、先に \(\boldsymbol{k}_2\) を1次の項までテイラー展開しておきます。

$$\begin{align}\boldsymbol{k}_2 &= \boldsymbol{f}\Bigl(t + \Delta t, \boldsymbol{x}(t) + \boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\Delta t\bigr)\Bigr)\\ &= \boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\Delta t\bigr) + \frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial t}\Delta t + \frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial\boldsymbol{x}}\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\Delta t\end{align}$$

\(\boldsymbol{k}_1\) とテイラー展開した \(\boldsymbol{k}_2\) を3番目更新式に代入して変形すると、以下の式になります。

$$\begin{align}\boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + (a + b)\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\Delta t\\ &+ b\Biggl(\frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial t} + \frac{\partial\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)}{\partial\boldsymbol{x}}\boldsymbol{f}\bigl(t, \boldsymbol{x}(t)\bigr)\Biggr)\Delta t^2\end{align}$$

同じ様な形になりましたね。未知定数の部分の比較をすると、\(a + b = 1 / 2\) 、\(b = 1 / 2\) であることが分かるので、

$$a = \frac{1}{2} \quad b = \frac{1}{2}$$

とすれば良いことが分かります。これらの未知定数を更新式に代入して整理すると、アルゴリズムの部分で示した以下の更新式が導かれます。

$$\begin{align}\boldsymbol{k}_1 &= \boldsymbol{f} \bigl(t, \boldsymbol{x}(t)\bigr)\\ \boldsymbol{k}_2 &= \boldsymbol{f} \bigl(t + \Delta t, \boldsymbol{x}(t) + \boldsymbol{k}_1\Delta t\bigr)\\ \boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \frac{\Delta t}{2}(\boldsymbol{k}_1 + \boldsymbol{k}_2)\end{align}$$

途中示したテイラー展開の式から、ホイン法が計算精度2次の手法であることが分かります。

実装


オイラー法のメソッドと同様な形になるように、以下のように実装しました。Solverクラスの中からホイン法の部分だけを抜き出してきています。

オイラー法と異なるのは更新式の形だけなので、メソッドの違いも10行目で使用している更新値の計算メソッドだけです。今回使用したHeunUpdateメソッドの実装は15行目から19行目にあり、示した更新式のように \(\boldsymbol{k}_1\) と \(\boldsymbol{k}_2\) を計算して、更新値を返すようにしています。

 

まとめ


今回はオイラー法よりも計算精度の良いホイン法の導出と実装を行いました。

ホイン法の更新式は以下になります。

$$\begin{align}\boldsymbol{k}_1 &= \boldsymbol{f} \bigl(t, \boldsymbol{x}(t)\bigr)\\ \boldsymbol{k}_2 &= \boldsymbol{f} \bigl(t + \Delta t, \boldsymbol{x}(t) + \boldsymbol{k}_1\Delta t\bigr)\\ \boldsymbol{x}(t + \Delta t) &= \boldsymbol{x}(t) + \frac{\Delta t}{2}(\boldsymbol{k}_1 + \boldsymbol{k}_2)\end{align}$$

次回はルンゲクッタ法の紹介を行いますが、導出は書くと大変なので割愛させていただきます。導出を行っている他のサイトさんを見つけたのでそちらの紹介を行い、アルゴリズムの紹介とその実装をメインに説明したいと思います。

次回 → 数値シミュレーション -ルンゲ=クッタ法-