はじめに

　本記事では、常微分方程式を解く解法の一つであるRKG法をPythonのライブラリNumpyを用いて解説していきます。

RKG法とは

常微分方程式（ODE: Ordinary Differential Equation）は、自然科学や工学の様々な分野で現れる基本的な数理モデルです。これらを数値的に解くための手法として、「Runge-Kutta法」が広く知られています。本稿では、その中でもGillによって提案された「Runge-Kutta-Gill法（RKG法）」について、わかりやすく紹介します。

常微分方程式の基本形

まず、1階の常微分方程式の一般的な形は次のようになります。

$$ \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0 $$

この方程式に対して、時間 $t$ をステップ幅 $h$ で刻みながら、$y$ の値を逐次求めていくのが数値解法の基本的な考え方です。

Runge-Kutta-Gill法の概要

RKG法は、4次のRunge-Kutta法の一種であり、標準的なRK4法に比べて丸め誤差の影響を抑える工夫が施されています。

Gill法のステップは以下のように定義されます：

$$ k_1 = h f(t_n, y_n) $$

$$ k_2 = h f\left(t_n + \frac{h}{2},\ y_n + \frac{1}{2}k_1\right) $$

$$ k_3 = h f\left(t_n + \frac{h}{2},\ y_n + \frac{1}{2}\left(\frac{\sqrt{2}-1}{2}k_1 + \frac{2-\sqrt{2}}{2}k_2\right)\right) $$

$$ k_4 = h f\left(t_n + h,\ y_n - \frac{\sqrt{2}}{2}k_2 + \frac{2+\sqrt{2}}{2}k_3\right) $$

最終的な更新式は：

$$ y_{n+1} = y_n + \frac{1}{6}(k_1 + (2 - \sqrt{2})k_2 + (2 + \sqrt{2})k_3 + k_4) $$

この手法では、$k_2$ および $k_3$ の係数に $\sqrt{2}$ が登場することで、浮動小数点演算における誤差の伝播を低減しています。

Runge-Kutta-Gill法は、標準的な4次RK法と同程度の計算コストでありながら、浮動小数点誤差に強いという特長を持っています。特に長時間にわたる数値積分や、厳密な精度が求められる場面で有効な選択肢となるでしょう。

PythonでRKG法を実装してみる

ここでは、PythonでRKG法を実装した関数を紹介します。 NumPyを使ってベクトル・スカラーどちらにも対応できるようにしています。

import numpy as np

def rkg(
    f: callable,
    y0: np.ndarray,
    t: np.ndarray
) -> np.ndarray:
    """
    RKG法で常微分方程式を数値的に解く

    Parameters
    ----------
    f : function
        dy/dt = f(t, y) を返す関数
    y0 : float or np.ndarray
        初期値
    t : np.ndarray
        時刻配列（等間隔）

    Returns
    -------
    y : np.ndarray
        各時刻における数値解
    """
    y0 = np.atleast_1d(y0)
    y = np.zeros((len(t), y0.size))
    y[0] = y0
    h = t[1] - t[0]
    sqrt2 = np.sqrt(2)
    for n in range(len(t)-1):
        k1 = h * f(t[n], y[n])
        k2 = h * f(t[n] + h/2, y[n] + k1/2)
        k3 = h * f(t[n] + h/2, y[n] + (sqrt2-1)*k1/2 + (2-sqrt2)*k2/2)
        k4 = h * f(t[n] + h, y[n] - sqrt2*k2/2 + (1+sqrt2/2)*k3)
        y[n+1] = y[n] + (k1 + (2-sqrt2)*k2 + (2+sqrt2)*k3 + k4)/6
    # y0がスカラーなら1次元配列で返す
    return y.squeeze()

実際に常微分方程式をPythonで解く

それでは、RKG法を使って実際に常微分方程式

$$ \frac{dy}{dt} = -y, \quad y(0) = 1 $$

を解いてみます。解析解（$y = e^{-t}$）と比較し、グラフで可視化します。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用（なくてもOK）

# RKG法の関数（上記参照）

# 微分方程式 dy/dt = -y
def f(t, y):
    return -y

# 初期値
y0 = 1.0

# 時刻配列
t = np.linspace(0, 5, 101)

# RKG法で数値解を求める
y = rkg(f, y0, t)

# 解析解
y_true = np.exp(-t)

# 結果のプロット
plt.plot(t, y, label='RKG法')
plt.plot(t, y_true, '--', label='解析解')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.title('dy/dt = -y のRKG法による数値解')
plt.show()