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

2023-06-24 / tau / コメントする

概要

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ステップ目のyとvを計算する。

(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ステップ目のΔv₁→Δy₁を計算し、それらを使ってΔv_i、Δy_i (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で以下の様に書いた。

# 2階常微分方程式
# - y'' = f(x, y, ...)の形の方程式が対象
# - ある計算ステップのx(i), y(i), y'(i)を与え、次のステップのy(i+1), y'(i+1)の値を返す
# - 呼び出し側で、右辺f=func(x, y, v, param1, param2, ...)のように定義されている必要がある
#   ただしv=y'
# 使用例
# - y'' = derivative(x, y, param1, param2)の場合
# - params = {'param1': 2.0, 'param2': 3.14}のようにパラメーターを定義
# - 以下のように呼び出して、ステップ計算を進める
# - (y_next, v_next) = euler2(x=x_now, y=y_now, v=v_now, derivative, params)あるいは
# - (y_next, v_next) = rk4_2(x=x_now, y=y_now, v=v_now, derivative, params)
# funcに与える関数がパラメータを持たない場合、上の引数でparams=Noneを指定

# Euler法により2階常微分微分方程式を解く
def euler2(x, y, v, dx, func, params):
    if params is None:
        dv = dx * func(x=x, y=y, v=v)
        dy = dx * v
    else:
        dv = dx * func(x=x, y=y, v=v, **params)
        dy = dx * v
    return (y + dy, v + dv)

# 4次のRunge-Kutta法により2階常微分方程式を解く
def rk4_2(x, y, v, dx, func, params):
    if params is None:
        dv1 = dx * func(x=x, y=y, v=v)
        dy1 = dx * v
        dv2 = dx * func(x=x+dx/2, y=y+dy1/2, v=v+dv1/2)
        dy2 = dx * (v + dv1 / 2)
        dv3 = dx * func(x=x+dx/2, y=y+dy2/2, v=v+dv2/2)
        dy3 = dx * (v + dv2 / 2)
        dv4 = dx * func(x=x+dx, y=y+dy3, v=v+dv3)
        dy4 = dx * (v + dv3)
    else:
        dv1 = dx * func(x=x, y=y, v=v, **params)
        dy1 = dx * v
        dv2 = dx * func(x=x+dx/2, y=y+dy1/2, v=v+dv1/2, **params)
        dy2 = dx * (v + dv1 / 2)
        dv3 = dx * func(x=x+dx/2, y=y+dy2/2, v=v+dv2/2, **params)
        dy3 = dx * (v + dv2 / 2)
        dv4 = dx * func(x=x+dx, y=y+dy3, v=v+dv3, **params)
        dy4 = dx * (v + dv3)
    return (
        y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6,
        v + (dv1 + 2 * dv2 + 2 * dv3 + dv4) / 6
    )

# 2階常微分方程式

# - y'' = f(x, y, ...)の形の方程式が対象

# - ある計算ステップのx(i), y(i), y'(i)を与え、次のステップのy(i+1), y'(i+1)の値を返す

# - 呼び出し側で、右辺f=func(x, y, v, param1, param2, ...)のように定義されている必要がある

# ただしv=y'

# 使用例

# - y'' = derivative(x, y, param1, param2)の場合

# - params = {'param1': 2.0, 'param2': 3.14}のようにパラメーターを定義

# - 以下のように呼び出して、ステップ計算を進める

# - (y_next, v_next) = euler2(x=x_now, y=y_now, v=v_now, derivative, params)あるいは

# - (y_next, v_next) = rk4_2(x=x_now, y=y_now, v=v_now, derivative, params)

# funcに与える関数がパラメータを持たない場合、上の引数でparams=Noneを指定

# Euler法により2階常微分微分方程式を解く

def euler2(x, y, v, dx, func, params):

if params is None:

dv = dx * func(x=x, y=y, v=v)

dy = dx * v

else:

dv = dx * func(x=x, y=y, v=v, **params)

dy = dx * v

return (y + dy, v + dv)

# 4次のRunge-Kutta法により2階常微分方程式を解く

def rk4_2(x, y, v, dx, func, params):

if params is None:

dv1 = dx * func(x=x, y=y, v=v)

dy1 = dx * v

dv2 = dx * func(x=x+dx/2, y=y+dy1/2, v=v+dv1/2)

dy2 = dx * (v + dv1 / 2)

dv3 = dx * func(x=x+dx/2, y=y+dy2/2, v=v+dv2/2)

dy3 = dx * (v + dv2 / 2)

dv4 = dx * func(x=x+dx, y=y+dy3, v=v+dv3)

dy4 = dx * (v + dv3)

else:

dv1 = dx * func(x=x, y=y, v=v, **params)

dy1 = dx * v

dv2 = dx * func(x=x+dx/2, y=y+dy1/2, v=v+dv1/2, **params)

dy2 = dx * (v + dv1 / 2)

dv3 = dx * func(x=x+dx/2, y=y+dy2/2, v=v+dv2/2, **params)

dy3 = dx * (v + dv2 / 2)

dv4 = dx * func(x=x+dx, y=y+dy3, v=v+dv3, **params)

dy4 = dx * (v + dv3)

return (

y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6,

v + (dv1 + 2 * dv2 + 2 * dv3 + dv4) / 6

)

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

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

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

戻り値は実数型の数値。

使い方

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

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

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

(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は上で示したモジュールを収めている。

import numpy as np
import matplotlib.pyplot as plt
import diffeq

import numpy as np

import matplotlib.pyplot as plt

import diffeq

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

# 微分方程式の右辺の関数
def derivative(x, y, v, a):
    return 2 * a

# 関数に渡すパラメーター
params = {'a': 1.0}

# 微分方程式の右辺の関数

def derivative(x, y, v, a):

return 2 * a

# 関数に渡すパラメーター

params = {'a': 1.0}

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

# 計算に必要な値の設定
x_start = -1.0  # 計算開始点のxの値
x_end = 1.0     # 計算終了点のxの値
dx = 0.2        # 分割幅
y_initial = 1.0 # 計算開始点におけるyの初期値
v_initial = -2.0# 計算開始点におけるvの値

# 計算に必要な値の設定

x_start = -1.0 # 計算開始点のxの値

x_end = 1.0 # 計算終了点のxの値

dx = 0.2 # 分割幅

y_initial = 1.0 # 計算開始点におけるyの初期値

v_initial = -2.0# 計算開始点におけるvの値

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

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

# 以下、計算用リストの確保
# 計算点数
num_points = int((x_end - x_start) / dx) + 1
# for分で回すためのxのリスト
# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない
x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)
# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)
y_euler = np.linspace(x_start, x_end, num=num_points)
v_euler = [0.0] * y_euler.size
y_rk = [0.0] * y_euler.size
v_runge_kutta = [0.0] * y_euler.size

# 計算開始点に初期値をセット
y_euler[0] = y_initial
v_euler[0] = v_initial
y_rk[0] = y_initial
v_runge_kutta[0] = v_initial

# 以下、計算用リストの確保

# 計算点数

num_points = int((x_end - x_start) / dx) + 1

# for分で回すためのxのリスト

# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない

x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)

# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)

y_euler = np.linspace(x_start, x_end, num=num_points)

v_euler = [0.0] * y_euler.size

y_rk = [0.0] * y_euler.size

v_runge_kutta = [0.0] * y_euler.size

# 計算開始点に初期値をセット

y_euler[0] = y_initial

v_euler[0] = v_initial

y_rk[0] = y_initial

v_runge_kutta[0] = v_initial

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

# 計算開始
for i, x in enumerate(x_calc):
    # Euler法
    (y_euler[i + 1], v_euler[i + 1]) = \
        diffeq.euler2(x=x, y=y_euler[i], v=v_euler[i], dx=dx,
                      func=derivative, params=params)
    # Runge-Kutta法
    (y_rk[i + 1], v_runge_kutta[i + 1]) = \
        diffeq.rk4_2(x=x, y=y_rk[i], v=v_runge_kutta[i], dx=dx,
                            func=derivative, params=params)

# 計算開始

for i, x in enumerate(x_calc):

# Euler法

(y_euler[i + 1], v_euler[i + 1]) = \

diffeq.euler2(x=x, y=y_euler[i], v=v_euler[i], dx=dx,

func=derivative, params=params)

# Runge-Kutta法

(y_rk[i + 1], v_runge_kutta[i + 1]) = \

diffeq.rk4_2(x=x, y=y_rk[i], v=v_runge_kutta[i], dx=dx,

func=derivative, params=params)

実行例

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

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

# 放物線の1階常微分方程式

import numpy as np
import matplotlib.pyplot as plt
import diffeq

# 解析解を計算する関数
def y_calc(x, a):
    return a * x**2

# 微分方程式の右辺の関数
def derivative(x, y, v, a):
    return 2 * a

# 関数に渡すパラメーター
params = {'a': 1.0}

# 計算に必要な値の設定
x_start = -1.0  # 計算開始点のxの値
x_end = 1.0     # 計算終了点のxの値
dx = 0.1        # 分割幅
y_initial = 1.0 # 計算開始点におけるyの初期値
v_initial = -2.0# 計算開始点におけるvの値

# 以下、計算用リストの確保
# 計算点数
num_points = int((x_end - x_start) / dx) + 1
# for分で回すためのxのリスト
# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない
x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)
# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)
y_euler = np.linspace(x_start, x_end, num=num_points)
v_euler = [0.0] * y_euler.size
y_rk = [0.0] * y_euler.size
v_runge_kutta = [0.0] * y_euler.size

# 計算開始点に初期値をセット
y_euler[0] = y_initial
v_euler[0] = v_initial
y_rk[0] = y_initial
v_runge_kutta[0] = v_initial

# 比較用の解析解表示のための定義域と計算値
x_domain = np.linspace(x_start, x_end, num=num_points)
y = y_calc(x=x_domain, a=params['a'])

# 計算開始
for i, x in enumerate(x_calc):
    # Euler法
    (y_euler[i + 1], v_euler[i + 1]) = \
        diffeq.euler2(x=x, y=y_euler[i], v=v_euler[i], dx=dx,
                      func=derivative, params=params)
    # Runge-Kutta法
    (y_rk[i + 1], v_runge_kutta[i + 1]) = \
        diffeq.rk4_2(x=x, y=y_rk[i], v=v_runge_kutta[i], dx=dx,
                            func=derivative, params=params)

# グラフ描画
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(x_domain, y, label="Solution", color="blue", linestyle="dashed", linewidth=4)
ax.plot(x_domain, y_euler, label="Euler", color="green", linestyle="solid", linewidth=2)
ax.plot(x_domain, y_rk, label="Runge-Kutta", color="red", linestyle="solid", linewidth=2)
ax.set_aspect('equal')
ax.text(0, 1, 'dx=' + str(dx))
ax.legend(loc="upper right")
plt.show()

# 放物線の1階常微分方程式

import numpy as np

import matplotlib.pyplot as plt

import diffeq

# 解析解を計算する関数

def y_calc(x, a):

return a * x**2

# 微分方程式の右辺の関数

def derivative(x, y, v, a):

return 2 * a

# 関数に渡すパラメーター

params = {'a': 1.0}

# 計算に必要な値の設定

x_start = -1.0 # 計算開始点のxの値

x_end = 1.0 # 計算終了点のxの値

dx = 0.1 # 分割幅

y_initial = 1.0 # 計算開始点におけるyの初期値

v_initial = -2.0# 計算開始点におけるvの値

# 以下、計算用リストの確保

# 計算点数

num_points = int((x_end - x_start) / dx) + 1

# for分で回すためのxのリスト

# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない

x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)

# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)

y_euler = np.linspace(x_start, x_end, num=num_points)

v_euler = [0.0] * y_euler.size

y_rk = [0.0] * y_euler.size

v_runge_kutta = [0.0] * y_euler.size

# 計算開始点に初期値をセット

y_euler[0] = y_initial

v_euler[0] = v_initial

y_rk[0] = y_initial

v_runge_kutta[0] = v_initial

# 比較用の解析解表示のための定義域と計算値

x_domain = np.linspace(x_start, x_end, num=num_points)

y = y_calc(x=x_domain, a=params['a'])

# 計算開始

for i, x in enumerate(x_calc):

# Euler法

(y_euler[i + 1], v_euler[i + 1]) = \

diffeq.euler2(x=x, y=y_euler[i], v=v_euler[i], dx=dx,

func=derivative, params=params)

# Runge-Kutta法

(y_rk[i + 1], v_runge_kutta[i + 1]) = \

diffeq.rk4_2(x=x, y=y_rk[i], v=v_runge_kutta[i], dx=dx,

func=derivative, params=params)

# グラフ描画

fig, ax = plt.subplots(figsize=(12, 6))

ax.plot(x_domain, y, label="Solution", color="blue", linestyle="dashed", linewidth=4)

ax.plot(x_domain, y_euler, label="Euler", color="green", linestyle="solid", linewidth=2)

ax.plot(x_domain, y_rk, label="Runge-Kutta", color="red", linestyle="solid", linewidth=2)

ax.set_aspect('equal')

ax.text(0, 1, 'dx=' + str(dx))

ax.legend(loc="upper right")

plt.show()

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

2023-06-24 / tau / コメントする

概要

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

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

基本形

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

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

Euler法

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

(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*}$

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

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

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

数値解法モジュール

コード

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

# 1階常微分方程式
# - y' = f(x, y, ...)の形の方程式が対象
# - ある計算ステップのx(i), y(i)を与え、次のステップのy(i+1)の値を返す
# - 呼び出し側で、右辺f=func(x, y, param1, param2, ...)のように定義されている必要がある
# 使用例
# - dy/dt = derivative(x, y, param1, param2)の場合
# - params = {'param1': 2.0, 'param2': 3.14}のようにパラメーターを定義
# - 以下のように呼び出して、ステップ計算を進める
# - y_next = euler1(x=x_now, y=y_now, derivative, params)あるいは
# - y_next = rk4_1(x=x_now, y=y_now, derivative, params)
# funcに与える関数がパラメータを持たない場合、上の引数でparams=Noneを指定

# Euler法により1階常微分微分方程式を解く
def euler1(x, y, dx, func, params):
    if params is None:
        return y + dx * func(x=x, y=y)
    else:
        return y + dx * func(x=x, y=y, **params)

# 4次のRunge-Kutta法により1階常微分方程式を解く
def rk4_1(x, y, dx, func, params):
    if params is None:
        dy1 = dx * func(x=x, y=y)
        dy2 = dx * func(x=x+dx/2, y=y+dy1/2)
        dy3 = dx * func(x=x+dx/2, y=y+dy2/2)
        dy4 = dx * func(x+dx, y+dy3)
        return y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6
    else:
        dy1 = dx * func(x=x, y=y, **params)
        dy2 = dx * func(x=x+dx/2, y=y+dy1/2, **params)
        dy3 = dx * func(x=x+dx/2, y=y+dy2/2, **params)
        dy4 = dx * func(x+dx, y+dy3, **params)
        return y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6

# 1階常微分方程式

# - y' = f(x, y, ...)の形の方程式が対象

# - ある計算ステップのx(i), y(i)を与え、次のステップのy(i+1)の値を返す

# - 呼び出し側で、右辺f=func(x, y, param1, param2, ...)のように定義されている必要がある

# 使用例

# - dy/dt = derivative(x, y, param1, param2)の場合

# - params = {'param1': 2.0, 'param2': 3.14}のようにパラメーターを定義

# - 以下のように呼び出して、ステップ計算を進める

# - y_next = euler1(x=x_now, y=y_now, derivative, params)あるいは

# - y_next = rk4_1(x=x_now, y=y_now, derivative, params)

# funcに与える関数がパラメータを持たない場合、上の引数でparams=Noneを指定

# Euler法により1階常微分微分方程式を解く

def euler1(x, y, dx, func, params):

if params is None:

return y + dx * func(x=x, y=y)

else:

return y + dx * func(x=x, y=y, **params)

# 4次のRunge-Kutta法により1階常微分方程式を解く

def rk4_1(x, y, dx, func, params):

if params is None:

dy1 = dx * func(x=x, y=y)

dy2 = dx * func(x=x+dx/2, y=y+dy1/2)

dy3 = dx * func(x=x+dx/2, y=y+dy2/2)

dy4 = dx * func(x+dx, y+dy3)

return y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6

else:

dy1 = dx * func(x=x, y=y, **params)

dy2 = dx * func(x=x+dx/2, y=y+dy1/2, **params)

dy3 = dx * func(x=x+dx/2, y=y+dy2/2, **params)

dy4 = dx * func(x+dx, y+dy3, **params)

return y + (dy1 + 2 * dy2 + 2 * dy3 + dy4) / 6

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

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

戻り値は実数型の数値。

使い方

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

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

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

(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は上で示したモジュールを収めている。

import numpy as np
import matplotlib.pyplot as plt
import diffeq

import numpy as np

import matplotlib.pyplot as plt

import diffeq

# 微分方程式の右辺の関数
def derivative(x, y, a):
    return 2 * a * x

# 関数に渡すパラメーター
params = {'a': 1.0}

# 微分方程式の右辺の関数

def derivative(x, y, a):

return 2 * a * x

# 関数に渡すパラメーター

params = {'a': 1.0}

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

# 計算に必要な値の設定
x_start = -1.0  # 計算開始点のxの値
x_end = 1.0     # 計算終了点のxの値
dx = 0.2        # 分割幅
y_initial = 1.0 # 計算開始点におけるyの初期値

# 計算に必要な値の設定

x_start = -1.0 # 計算開始点のxの値

x_end = 1.0 # 計算終了点のxの値

dx = 0.2 # 分割幅

y_initial = 1.0 # 計算開始点におけるyの初期値

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

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

# 以下、計算用リストの確保
# 計算点数
num_points = int((x_end - x_start) / dx) + 1
# for分で回すためのxのリスト
# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない
x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)
# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)
y_euler = np.linspace(x_start, x_end, num=num_points)
y_rk = [0.0] * y_euler.size

# 計算開始点に初期値をセット
y_euler[0] = y_initial
y_rk[0] = y_initial

# 以下、計算用リストの確保

# 計算点数

num_points = int((x_end - x_start) / dx) + 1

# for分で回すためのxのリスト

# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない

x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)

# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)

y_euler = np.linspace(x_start, x_end, num=num_points)

y_rk = [0.0] * y_euler.size

# 計算開始点に初期値をセット

y_euler[0] = y_initial

y_rk[0] = y_initial

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

# 計算開始
for i, x in enumerate(x_calc):
    # Euler法
    y_euler[i + 1] = \
        diffeq.euler1(x=x, y=y_euler[i], dx=dx, func=derivative, params=params)
    # Runge-Kutta法
    y_rk[i + 1] = \
        diffeq.rk4_1(x=x, y=y_rk[i], dx=dx, func=derivative, params=params)

# 計算開始

for i, x in enumerate(x_calc):

# Euler法

y_euler[i + 1] = \

diffeq.euler1(x=x, y=y_euler[i], dx=dx, func=derivative, params=params)

# Runge-Kutta法

y_rk[i + 1] = \

diffeq.rk4_1(x=x, y=y_rk[i], dx=dx, func=derivative, params=params)

実行例

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

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

# 放物線の1階常微分方程式

import numpy as np
import matplotlib.pyplot as plt
import diffeq

# 解析解を計算する関数
def y_calc(x, a):
    return a * x**2

# 微分方程式の右辺の関数
def derivative(x, y, a):
    return 2 * a * x

# 関数に渡すパラメーター
params = {'a': 1.0}

# 計算に必要な値の設定
x_start = -1.0  # 計算開始点のxの値
x_end = 1.0     # 計算終了点のxの値
dx = 0.1        # 分割幅
y_initial = 1.0 # 計算開始点におけるyの初期値

# 以下、計算用リストの確保
# 計算点数
num_points = int((x_end - x_start) / dx) + 1
# for分で回すためのxのリスト
# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない
x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)
# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)
y_euler = np.linspace(x_start, x_end, num=num_points)
y_rk = [0.0] * y_euler.size

# 計算開始点に初期値をセット
y_euler[0] = y_initial
y_rk[0] = y_initial

# 比較用の解析解表示のための定義域と計算値
x_domain = np.linspace(x_start, x_end, num=num_points)
y = y_calc(x=x_domain, a=params['a'])

# 計算開始
for i, x in enumerate(x_calc):
    # Euler法
    y_euler[i + 1] = \
        diffeq.euler1(x=x, y=y_euler[i], dx=dx, func=derivative, params=params)
    # Runge-Kutta法
    y_rk[i + 1] = \
        diffeq.rk4_1(x=x, y=y_rk[i], dx=dx, func=derivative, params=params)

# グラフ描画
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(x_domain, y, label="Solution", color="blue", linestyle="dashed", linewidth=4)
ax.plot(x_domain, y_euler, label="Euler", color="green", linestyle="solid", linewidth=2)
ax.plot(x_domain, y_rk, label="Runge-Kutta", color="red", linestyle="solid", linewidth=2)
ax.set_aspect('equal')
ax.text(0, 1, 'dx=' + str(dx))
ax.legend(loc="upper right")
plt.show()

# 放物線の1階常微分方程式

import numpy as np

import matplotlib.pyplot as plt

import diffeq

# 解析解を計算する関数

def y_calc(x, a):

return a * x**2

# 微分方程式の右辺の関数

def derivative(x, y, a):

return 2 * a * x

# 関数に渡すパラメーター

params = {'a': 1.0}

# 計算に必要な値の設定

x_start = -1.0 # 計算開始点のxの値

x_end = 1.0 # 計算終了点のxの値

dx = 0.1 # 分割幅

y_initial = 1.0 # 計算開始点におけるyの初期値

# 以下、計算用リストの確保

# 計算点数

num_points = int((x_end - x_start) / dx) + 1

# for分で回すためのxのリスト

# n-1ステップの計算が最後でnステップ目のyを計算するため、xの終了点は含まない

x_calc = np.linspace(x_start, x_end, num=num_points-1, endpoint=False)

# Euler法・Runge-Kutta法の数値解yの保存領域(終了点を含む)

y_euler = np.linspace(x_start, x_end, num=num_points)

y_rk = [0.0] * y_euler.size

# 計算開始点に初期値をセット

y_euler[0] = y_initial

y_rk[0] = y_initial

# 比較用の解析解表示のための定義域と計算値

x_domain = np.linspace(x_start, x_end, num=num_points)

y = y_calc(x=x_domain, a=params['a'])

# 計算開始

for i, x in enumerate(x_calc):

# Euler法

y_euler[i + 1] = \

diffeq.euler1(x=x, y=y_euler[i], dx=dx, func=derivative, params=params)

# Runge-Kutta法

y_rk[i + 1] = \

diffeq.rk4_1(x=x, y=y_rk[i], dx=dx, func=derivative, params=params)

# グラフ描画

fig, ax = plt.subplots(figsize=(12, 6))

ax.plot(x_domain, y, label="Solution", color="blue", linestyle="dashed", linewidth=4)

ax.plot(x_domain, y_euler, label="Euler", color="green", linestyle="solid", linewidth=2)

ax.plot(x_domain, y_rk, label="Runge-Kutta", color="red", linestyle="solid", linewidth=2)

ax.set_aspect('equal')

ax.text(0, 1, 'dx=' + str(dx))

ax.legend(loc="upper right")

plt.show()

TauStation

タグ: differential equation

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

概要

基本形

Euler法

Runge-Kutta法

数値解法モジュール

コード

使い方

実行例

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

概要

基本形

Euler法

Runge-Kutta法

数値解法モジュール

コード

使い方

実行例