https://pydocument.hatenablog.com/entry/2023/05/08/232915

この記事では、オイラー法を用いて微分方程式を数値的に解く方法を解説します。主に、Pythonを使用した実装例を通じて、その手順と注意点について説明します。

オイラー法とは

オイラー法は、微分方程式の数値解法の中で最も基本的な手法の一つです。この方法では、微小な時間ステップ $(\Delta t)$ を用いて、微分方程式を差分方程式に近似し、未来の値を逐次的に計算します。

具体例

以下の初期値問題の常微分方程式を考えます。

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

ここで、 $t$ は時間、 $y$ は求めたい関数、 $f(t, y)$ は既知の関数、 $y_0$ は初期値です。

オイラー法では、微小時間 $\Delta t$ 後の $y$ の値 $(y_{i+1})$ を以下の式で近似します。

$y_{i+1} = y_i + f(t_i, y_i) \Delta t$

ここで、 $y_i$ は時刻 $t_i$ での $y$ の値、 $t_{i+1} = t_i + \Delta t$ です。

Pythonによる実装例

例題

次の微分方程式をオイラー法で解きます。

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

この方程式の厳密解は

$y = t - 1 + 2e^{-t}$

です。

1. ライブラリのインポート

まず、数値計算に必要なNumPyと、結果の描画に使うMatplotlibをインポートします。

import numpy as np
import matplotlib.pyplot as plt

2. 微分方程式と初期値、時間ステップの設定

次に、微分方程式の右辺を表す関数 f、初期値 y0、計算する時間の範囲 t_start から t_end、時間刻み幅 dt を設定します。

# 微分方程式の右辺を定義
def f(t, y):
    return t - y

# 初期値
y0 = 1

# 時間の範囲と刻み幅
t_start = 0
t_end = 5
dt = 0.1  # 刻み幅を小さくすると精度が上がるが、計算時間が増加する

3. オイラー法の基本的な実装

時間ステップごとの y の値を格納するための配列を用意し、オイラー法の式に従って値を更新していきます。 np.arange(t_start, t_end + dt, dt)によりt_start から t_end まで dt 刻みで時間 t の配列を作成します。t_end を計算に含めるため、t_end + dt としています。

# 時間の配列を作成
t = np.arange(t_start, t_end + dt, dt) # t_end も含むように修正
# yの値を格納する配列を初期化
y = np.zeros(len(t))
y[0] = y0  # 初期値を設定

# オイラー法の計算
for i in range(len(t) - 1):
    y[i+1] = y[i] + f(t[i], y[i]) * dt

4. オイラー法の関数化

オイラー法の処理を関数 euler としてまとめ、再利用性を高めます。

def euler(f, y0, t):
    """
    オイラー法で常微分方程式を解く関数

    Args:
        f: 微分方程式の右辺を表す関数 (t, yを引数とする)
        y0: 初期値
        t: 時間の配列

    Returns:
        y: 数値解の配列
    """
    y = np.zeros(len(t))
    y[0] = y0
    dt = t[1] - t[0]  # 時間刻み幅
    for i in range(len(t) - 1):
        y[i+1] = y[i] + f(t[i], y[i]) * dt
    return y

5. 厳密解との比較 (グラフ化)

数値解と厳密解を比較するために、結果をグラフにプロットします。

# オイラー法で計算
t = np.arange(t_start, t_end + dt, dt) # 上で定義した dt を利用
y_euler = euler(f, y0, t)

# 厳密解を計算
y_exact = t - 1 + 2 * np.exp(-t)

# グラフを描画
plt.plot(t, y_euler, label='Euler Method', marker='o', linestyle='-')
plt.plot(t, y_exact, label='Exact Solution', linestyle='--')
plt.xlabel('t')
plt.ylabel('y')
plt.title('Comparison of Euler Method and Exact Solution')
plt.legend()
plt.grid(True)
plt.show()

6. 刻み幅の影響 (発展)

オイラー法の精度の確認のために、異なる刻み幅(dt)で計算を行い、数値解が厳密解にどのように近づくか（または遠ざかるか）を確認します。

# 異なる刻み幅で計算
dt_values = [0.5, 0.1, 0.01]
for dt in dt_values:
    t = np.arange(t_start, t_end + dt, dt)
    y_euler = euler(f, y0, t)
    plt.plot(t, y_euler, label=f'Euler (dt={dt})')

# 厳密解
y_exact = t - 1 + 2 * np.exp(-t)
t_fine = np.linspace(t_start,t_end,1000) # 厳密解を滑らかに描画するため
y_exact = t_fine - 1 + 2 * np.exp(-t_fine)
plt.plot(t_fine, y_exact, label='Exact Solution', linestyle='--')

plt.xlabel('t')
plt.ylabel('y')
plt.title('Effect of Step Size (dt) on Euler Method')
plt.legend()
plt.grid(True)
plt.show()

まとめと注意点、発展

オイラー法は非常にシンプルで実装が容易です。ただし、1次の精度しかないため、高い精度を得るには時間刻み幅 $(\Delta t)$ を非常に小さくする必要があります。これにより計算コストが増加します。オイラー法は、各ステップで局所的な誤差が生じ、それが蓄積することで全体の誤差（大域的誤差）が大きくなる可能性があります。より高精度な数値解法として、ルンゲ・クッタ法などがあります。

表形式での比較

手法	特徴	精度	計算コスト
オイラー法	最も基本的な手法。1ステップずつ未来の値を計算。実装が簡単。	1次精度	低い
ルンゲ・クッタ法	オイラー法よりも高精度な解法。中間のステップでの傾きを複数回計算することで、より正確な近似を行う。4次のルンゲ・クッタ法（RK4）がよく使われる。	4次精度など	高い
陰的解法	微分方程式の未来の値を使って現在の値を計算する方法。特に硬い微分方程式（変化の激しい方程式）に対して安定して解ける。反復計算が必要になる場合があるため、計算コストは高くなる傾向にある。	様々	高い

[PR]

click.linksynergy.com

問題解決のための「アルゴリズム×数学」が基礎からしっかり身につく本 [ 米田優峻 ]