未知関数 $x(t)$ とその導関数から構成される方程式を、常微分方程式（ODE, Ordinary Differential Equation）と言います。
多くの物理現象やエンジニアリングの問題は常微分方程式で記述されますが、解析的に解ける方程式は限られています。そこで数値的に解を求める手法が重要になります。その中で最も基本的な手法がオイラー法です。
本記事の内容
- 常微分方程式の初期値問題
- オイラー法のアルゴリズムと導出
- 誤差解析（局所誤差と大域誤差）
- 改良オイラー法（ホイン法）
- Pythonでの実装と比較
常微分方程式の初期値問題
1階常微分方程式の初期値問題は次の形で定式化されます。
$$ \frac{dx}{dt} = f(t, x), \quad x(t_0) = x_0 $$
ここで $f(t, x)$ は右辺関数、$x_0$ は初期条件です。
例: 指数減衰
$$ \frac{dx}{dt} = -\alpha x, \quad x(0) = x_0 $$
解析解は $x(t) = x_0 e^{-\alpha t}$ です。
オイラー法の導出
テイラー展開を $t_n$ の周りで行うと、
$$ x(t_{n+1}) = x(t_n) + h \cdot x'(t_n) + \frac{h^2}{2} x”(\xi_n) $$
ここで $h = t_{n+1} – t_n$ は時間刻み幅、$\xi_n \in [t_n, t_{n+1}]$ です。
2次以降の項を無視すると、オイラー法の漸化式が得られます。
$$ x_{n+1} = x_n + h \cdot f(t_n, x_n) $$
幾何学的解釈
オイラー法は、各時刻 $t_n$ における接線の傾き $f(t_n, x_n)$ を使って、次の時刻の値を線形外挿する方法です。
誤差解析
局所切断誤差（Local Truncation Error）
1ステップあたりの誤差は、テイラー展開の打ち切り項に対応します。
$$ \tau_n = \frac{h^2}{2} x”(\xi_n) = O(h^2) $$
大域誤差（Global Error）
$N = (T – t_0)/h$ ステップ後の累積誤差は、
$$ E_N = O(h) $$
つまり、オイラー法は1次精度の手法です。刻み幅を半分にすると、大域誤差もおよそ半分になります。
改良オイラー法（ホイン法）
オイラー法の精度を向上させたのが改良オイラー法（Heun’s method）です。
$$ \begin{align} \tilde{x}_{n+1} &= x_n + h \cdot f(t_n, x_n) \quad (\text{予測子}) \\ x_{n+1} &= x_n + \frac{h}{2}[f(t_n, x_n) + f(t_{n+1}, \tilde{x}_{n+1})] \quad (\text{修正子}) \end{align} $$
改良オイラー法は2次精度（$E_N = O(h^2)$）であり、同じ刻み幅でオイラー法より高精度です。
ルンゲ・クッタ法（4次）
さらに高精度な手法として、4次のルンゲ・クッタ法（RK4）があります。
$$ \begin{align} k_1 &= f(t_n, x_n) \\ k_2 &= f(t_n + h/2, x_n + hk_1/2) \\ k_3 &= f(t_n + h/2, x_n + hk_2/2) \\ k_4 &= f(t_n + h, x_n + hk_3) \\ x_{n+1} &= x_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{align} $$
RK4は4次精度（$E_N = O(h^4)$）で、多くの実用問題で標準的に使われています。
Pythonでの実装と比較
import numpy as np
import matplotlib.pyplot as plt
def euler_method(f, t0, x0, t_end, h):
"""オイラー法"""
t_values = [t0]
x_values = [x0]
t, x = t0, x0
while t < t_end - h/2:
x = x + h * f(t, x)
t = t + h
t_values.append(t)
x_values.append(x)
return np.array(t_values), np.array(x_values)
def improved_euler(f, t0, x0, t_end, h):
"""改良オイラー法（ホイン法）"""
t_values = [t0]
x_values = [x0]
t, x = t0, x0
while t < t_end - h/2:
k1 = f(t, x)
x_pred = x + h * k1
k2 = f(t + h, x_pred)
x = x + h / 2 * (k1 + k2)
t = t + h
t_values.append(t)
x_values.append(x)
return np.array(t_values), np.array(x_values)
def rk4_method(f, t0, x0, t_end, h):
"""4次ルンゲ・クッタ法"""
t_values = [t0]
x_values = [x0]
t, x = t0, x0
while t < t_end - h/2:
k1 = f(t, x)
k2 = f(t + h/2, x + h*k1/2)
k3 = f(t + h/2, x + h*k2/2)
k4 = f(t + h, x + h*k3)
x = x + h/6 * (k1 + 2*k2 + 2*k3 + k4)
t = t + h
t_values.append(t)
x_values.append(x)
return np.array(t_values), np.array(x_values)
# 問題1: 指数減衰 dx/dt = -x
f1 = lambda t, x: -x
exact1 = lambda t: np.exp(-t)
# 問題2: 振動 dx/dt = cos(t)
f2 = lambda t, x: np.cos(t)
exact2 = lambda t: np.sin(t)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# 問題1: 各手法の比較
h = 0.5
t_exact = np.linspace(0, 5, 500)
t_e, x_e = euler_method(f1, 0, 1, 5, h)
t_ie, x_ie = improved_euler(f1, 0, 1, 5, h)
t_rk, x_rk = rk4_method(f1, 0, 1, 5, h)
axes[0, 0].plot(t_exact, exact1(t_exact), 'k-', linewidth=2, label='Exact')
axes[0, 0].plot(t_e, x_e, 'ro-', label=f'Euler (h={h})', markersize=4)
axes[0, 0].plot(t_ie, x_ie, 'bs-', label=f'Improved Euler', markersize=4)
axes[0, 0].plot(t_rk, x_rk, 'g^-', label='RK4', markersize=4)
axes[0, 0].set_xlabel('t')
axes[0, 0].set_ylabel('x(t)')
axes[0, 0].set_title("dx/dt = -x, x(0) = 1")
axes[0, 0].legend()
axes[0, 0].grid(True)
# 問題1: 誤差の比較
axes[0, 1].semilogy(t_e, np.abs(x_e - exact1(t_e)), 'ro-', label='Euler', markersize=4)
axes[0, 1].semilogy(t_ie, np.abs(x_ie - exact1(t_ie)), 'bs-', label='Improved Euler', markersize=4)
axes[0, 1].semilogy(t_rk, np.abs(x_rk - exact1(t_rk)) + 1e-16, 'g^-', label='RK4', markersize=4)
axes[0, 1].set_xlabel('t')
axes[0, 1].set_ylabel('|Error|')
axes[0, 1].set_title('Absolute Error Comparison')
axes[0, 1].legend()
axes[0, 1].grid(True)
# 刻み幅と大域誤差の関係
h_values = np.logspace(-3, 0, 20)
errors_euler = []
errors_ie = []
errors_rk4 = []
t_end = 2.0
for h_val in h_values:
t_e, x_e = euler_method(f1, 0, 1, t_end, h_val)
t_ie, x_ie = improved_euler(f1, 0, 1, t_end, h_val)
t_rk, x_rk = rk4_method(f1, 0, 1, t_end, h_val)
errors_euler.append(np.abs(x_e[-1] - exact1(t_e[-1])))
errors_ie.append(np.abs(x_ie[-1] - exact1(t_ie[-1])))
errors_rk4.append(np.abs(x_rk[-1] - exact1(t_rk[-1])) + 1e-16)
axes[1, 0].loglog(h_values, errors_euler, 'ro-', label='Euler (O(h))', markersize=4)
axes[1, 0].loglog(h_values, errors_ie, 'bs-', label='Improved Euler (O(h^2))', markersize=4)
axes[1, 0].loglog(h_values, errors_rk4, 'g^-', label='RK4 (O(h^4))', markersize=4)
# 参照線
axes[1, 0].loglog(h_values, h_values, 'k--', alpha=0.3, label='O(h)')
axes[1, 0].loglog(h_values, h_values**2, 'k-.', alpha=0.3, label='O(h^2)')
axes[1, 0].loglog(h_values, h_values**4, 'k:', alpha=0.3, label='O(h^4)')
axes[1, 0].set_xlabel('Step Size h')
axes[1, 0].set_ylabel('Global Error at t=2')
axes[1, 0].set_title('Convergence Order')
axes[1, 0].legend(fontsize=8)
axes[1, 0].grid(True)
# 2体問題（振り子）への応用
def pendulum(t, state):
theta, omega = state
return np.array([omega, -9.81 * np.sin(theta)])
h = 0.01
t_end_pend = 10.0
n_steps = int(t_end_pend / h)
t_pend = np.zeros(n_steps + 1)
x_euler = np.zeros((n_steps + 1, 2))
x_rk4_pend = np.zeros((n_steps + 1, 2))
x_euler[0] = [np.pi / 4, 0] # 初期角度45度, 初期角速度0
x_rk4_pend[0] = [np.pi / 4, 0]
for i in range(n_steps):
t_pend[i + 1] = t_pend[i] + h
# オイラー法
x_euler[i + 1] = x_euler[i] + h * pendulum(t_pend[i], x_euler[i])
# RK4
k1 = pendulum(t_pend[i], x_rk4_pend[i])
k2 = pendulum(t_pend[i] + h/2, x_rk4_pend[i] + h*k1/2)
k3 = pendulum(t_pend[i] + h/2, x_rk4_pend[i] + h*k2/2)
k4 = pendulum(t_pend[i] + h, x_rk4_pend[i] + h*k3)
x_rk4_pend[i + 1] = x_rk4_pend[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
axes[1, 1].plot(t_pend, np.degrees(x_euler[:, 0]), 'r-', label='Euler', alpha=0.7)
axes[1, 1].plot(t_pend, np.degrees(x_rk4_pend[:, 0]), 'b-', label='RK4', linewidth=2)
axes[1, 1].set_xlabel('Time [s]')
axes[1, 1].set_ylabel('Angle [deg]')
axes[1, 1].set_title('Pendulum: Euler vs RK4')
axes[1, 1].legend()
axes[1, 1].grid(True)
plt.tight_layout()
plt.show()
まとめ
本記事では、常微分方程式をオイラー法で解く方法について解説しました。
- オイラー法はテイラー展開の1次近似に基づく最も基本的な数値解法である
- 漸化式 $x_{n+1} = x_n + h f(t_n, x_n)$ で計算でき、実装が非常に簡単
- 大域誤差は $O(h)$（1次精度）で、刻み幅を半分にすると誤差も半分になる
- 改良オイラー法は2次精度、ルンゲ・クッタ法（RK4）は4次精度で高精度
- 実用的には、精度と計算コストのバランスからRK4が標準的に使われる
次のステップとして、以下の記事も参考にしてください。