2階常微分方程式の数値解法

概要

2階常微分方程式の数値解法、Euler法とRunge-Kutta法をPythonのモジュールとしてまとめた。

1階常微分方程式の解法についてはこちら

基本形

以下のような方程式を考える。

(1)    \begin{equation*} \dfrac{d^2 y}{dx^2} = f\left ( x, y, \dfrac{dy}{dx}, ... \right ) \end{equation*}

(2)    \begin{equation*} y'' = f(x, y, y', ... ) \end{equation*}

y‘(x)をv(x)と置いて以下のように変形する。

(3)    \begin{align*} \frac{dy}{dx} &=y' = v \\ \frac{dv}{dx} &= v' = f(x, y, y', ...) \end{align*}

Euler法

Euler法では、iステップ目のΔv→Δyと計算し、それらを使ってi+1ステップ目のyvを計算する。

(4)    \begin{align*} \Delta v_i &= \Delata x \cdot f(x_i, y_i, y'_i, ...) \\ \Delta y_i &= \Delta x \cdot v_i \\ y_{i+1} &= y_i + \Delta y_i \\ v_{i+1} &= v_i + \Delta v_i \end{align*}

Runge-Kutta法

Runge-Kutta法も同様に、iステップ目のΔv1→Δy1を計算し、それらを使ってΔvi、Δyi (i = 1, 2, 3, 4)を計算していく。

(5)    \begin{align*} \Delta v_{1, i} &= \Delta x \cdot f(x_i, y_i, v_i, ...) \\ \Delta y_{1, i} &= \Delta x \cdot v_i \\ \Delta v_{2, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{1, i}}{2}, v_i + \dfrac{\Delta v_{1, i}}{2}, ... \right ) \\ \Delta y_{2, i} &= \Delta x \cdot \left (v_i + \dfrac{\Delta v_{1, i}}{2} \right ) \\ \Delta v_{3, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{2, i}}{2}, v_i + \dfrac{\Delta v_{2, i}}{2}, ... \right ) \\ \Delta y_{3, i} &= \Delta x \cdot \left (v_i + \dfrac{\Delta v_{2, i}}{2} \right ) \\ \Delta v_{4, i} &= \Delta x \cdot (v_i + \Delta v_{3, i}) \\ \Delta y_{4, i} &= \Delta x \cdot f(x_i + \Delta x, y_i + \Delta y_{3, i}, v_i + \Delta v_{3, i}, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \frac{\Delta y_{1, i} + 2 \Delta y_{2, i} + 2 \Delta y_{3, i} + \Delta y_{4, i}}{6} \\ v_{i+1} &= v_i + \frac{\Delta v_{1, i} + 2 \Delta v_{2, i} + 2 \Delta v_{3, i} + \Delta v_{4, i}}{6} \end{align*}

数値解法モジュール

コード

Euler法とRunge-Kutta法で2階常微分方程式を解くモジュールを、Pythonで以下の様に書いた。

関数euler2rk4_2の使い方はいずれも同じで、以下の通り。

  • x, y, v:現在の計算ステップでの独立変数x、xにおけるyの値、xにおけるyの1次微分値
  • dx:計算の分割幅
  • func:微分方程式の右辺の関数
  • params:関数funcに与えるパラメーター

パラメーターは、呼び出し側の関数定義で設定したパラメーター名をキーとした辞書で与える。その前提で、モジュール内で**paramsとしてアンパックしている。

戻り値は実数型の数値。

使い方

モジュールの使い方の概要は以下のとおり。

  1. 右辺の関数を定義する
  2. 関数に渡すパラメーターを設定する
  3. 計算開始点と終了点のxの値、分割幅dx, yvの初期値を設定する
  4. x, y, vの計算用リストを確保する
  5. 計算開始点のyvのリストに初期値をセットする
  6. 計算実行

単純な放物線の微分方程式を例にしてモジュールの使い方を示す。まず、解の方程式として以下を考える。

(6)    \begin{equation*} y = a x^2 \quad (-1 \le x \le 1) \end{equation*}

この式を2階微分した、以下の常微分方程式を解いていく。

(7)    \begin{equation*} y'' = 2a \quad \left( x = -1 \Rightarrow \left\{ \begin{array}{l} y = a\\ y' =-2a \end{array} \right) \end{equation*}

まず、必要なパッケージをインポート。diffeqは上で示したモジュールを収めている。

次に、微分方程式右辺の関数を定義する(以下のderivative)。この関数はパラメーターとして係数aを1つ持っている。また、パラメーターの値を辞書で定義する。

数値計算に計算開始点x_startと終了点x_end、分割幅dxyvの計算初期値を設定。

次に、予めxの計算点のリストを用意して、これに対応するy, vの値を計算・保存していく。計算点数は植木算で分割数+1とし、ループ計算に用いるリストは計算終了点を除いている。結果を保存しないなら、リストではなく逐次各ステップの計算値を次のステップの入力値にする。

また、計算開始点におけるy, vの初期値をセットしている。

計算自体はシンプルで、定義域リストの値を1つずつ進めながら、現在の値や関数への参照、パラメーターを与えてモジュールを呼び出す。

実行例

元の解析解と併せて表示した実行結果は以下の通りで、1回微分方程式の場合と同じく、Euler法の解は計算が進むごとに乖離が大きくなっている。Runge-Kutta法の解はかなり正確に解析解と一致している。

以下は、解析解も含めてグラフ表示するための、全体の実行コード。

 

1階常微分方程式の数値解法

概要

1階常微分方程式の数値解法、Euler法とRunge-Kutta法をPythonのモジュールとしてまとめた。

2階常微分方程式の解法についてはこちら

基本形

以下のような方程式を考える。

(1)    \begin{equation*} \dfrac{dy}{dx} = y' = f(x, y, ...) \end{equation*}

Euler法

Euler法はシンプルな方法で、ある計算ステップiの時の関数上の点(xi, yi)とその点における微分係数yiを使って、計算幅dxだけ先にある点xi+1 = xi + dxのに対応するyiの値を近似的に計算する。

(2)    \begin{align*} \Delta y_i &= \Delta x \cdot f(x_i, y_i, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \Delta y_i \end{align*}

概念図からもわかるように、たとえば下に凸の関数の場合、微分係数で線形的に予測した値は常に真値より小さくなるため、誤差が拡大していくことが予想される。

Runge-Kutta法

Runge-Kutta法には1次、2次など異なる次数の解法があるが、4次の方法がよく紹介されている。

(3)    \begin{align*} \Delta y_{1, i} &= \Delta x \cdot f(x_i, y_i, ...) \\ \Delta y_{2, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{1, i}}{2}, ... \right ) \\ \Delta y_{3, i} &= \Delta x \cdot f \left ( x_i + \dfrac{\Delta x}{2}, y_i + \dfrac{\Delta y_{2, i}}{2}, ... \right ) \\ \Delta y_{4, i} &= \Delta x \cdot f(x_i + \Delta x, y_i + \Delta y_{3, i}, ...) \\ x_{i+1} &= x_i + \Delta x \\ y_{i+1} &= y_i + \frac{\Delta y_{1, i} + 2 \Delta y_{2, i} + 2 \Delta y_{3, i} + \Delta y_{4, i}}{6} \end{align*}

これを説明するのに、xy平面上での関数を描き、その微分係数を複数描いているものが多い。

本来は、右辺のz = f(…)上のある曲線があって、それをxy平面に投影した曲線のある点における微分係数と、その点におけるzの値が等しくなるような曲線を求めるという意味となる。実際にzの値と投影曲線の微分係数が等しいというのを直感的に示すのはやっかいだなと思っていた。

いろいろ探してみると、どうやらこの式形は元の関数をTaylor展開して、次数を上げて精度を高めるというものらしい。文献を見てみると、その導出過程でかなり長々とした式が出ていた。時間がある時に試してみたい。

数値解法モジュール

コード

Euler法とRunge-Kutta法で1階常微分方程式を解くモジュールを、Pythonで以下の様に書いた。

関数euler1rk4_1の使い方はいずれも同じで、以下の通り。

  • x, y:現在の計算ステップでの独立変数xと関数値yの値
  • dx:計算の分割幅
  • func:微分方程式の右辺の関数
  • params:関数funcに与えるパラメーター

パラメーターは、呼び出し側の関数定義で設定したパラメーター名をキーとした辞書で与える。その前提で、モジュール内で**paramsとしてアンパックしている。

戻り値は実数型の数値。

使い方

モジュールの使い方の概要は以下のとおり。

  1. 右辺の関数を定義する
  2. 関数に渡すパラメーターを設定する
  3. 計算開始点と終了点のxの値、分割幅dx, yの初期値を設定する
  4. x, yの計算用リストを確保する
  5. 計算開始点のyのリストに初期値をセットする
  6. 計算実行

単純な放物線の微分方程式を例にしてモジュールの使い方を示す。まず、解の方程式として以下を考える。

(4)    \begin{equation*} y = a x^2 \quad (-1 \le x \le 1) \end{equation*}

この式を1階微分した、以下の常微分方程式を解いていく。

(5)    \begin{equation*} y'= 2ax \quad (x = -1 \Rightarrow y = 1) \end{equation*}

まず、必要なパッケージをインポート。diffeqは上で示したモジュールを収めている。

次に、微分方程式右辺の関数を定義する(以下のderivative)。この関数はパラメーターとして係数aを1つ持っている。また、パラメーターの値を辞書で定義する。

数値計算に計算開始点x_startと終了点x_end、分割幅dxyの計算初期値を設定。

次に、予めxの計算点のリストを用意して、これに対応するyの値を計算・保存していく。計算点数は植木算で分割数+1とし、ループ計算に用いるリストは計算終了点を除いている。結果を保存しないなら、リストではなく逐次各ステップの計算値を次のステップの入力値にする。

また、計算開始点におけるyの初期値をセットしている。

計算自体はシンプルで、定義域リストの値を1つずつ進めながら、現在の値や関数への参照、パラメーターを与えてモジュールを呼び出す。

実行例

元の解析解と併せて表示した実行結果は以下の通りで、予想通りEuler法の解は計算が進むごとに乖離が大きくなっている。Runge-Kutta法の解はかなり正確に解析解と一致している。

以下は、解析解も含めてグラフ表示するための、全体の実行コード。