Pythonでのオイラー法の実装

オイラー法は微分方程式を数値的に解く手法の一つです。この手法は、微小な時間刻み幅で微分方程式を近似的に解き、数値解を求めます。Pythonを使ってオイラー法を実装することができます。以下に、Pythonでのオイラー法の実装方法を具体例と共に説明します。この例では、次の微分方程式を解くことを考えます。

dy/dt = t - y 初期値:y(0) = 1

この微分方程式の厳密解は y = t - 1 + 2exp(-t) です。オイラー法を使って、この微分方程式の数値解を求めてみましょう。

必要なライブラリのインポート

まず、必要なライブラリをインポートします。

import matplotlib.pyplot as plt
import numpy as np

オイラー法の実装

次に、オイラー法を実装します。オイラー法では、微分方程式を以下のように近似的に解きます。

y[i+1] = y[i] + f(t[i], y[i]) * dt

ここで、y[i]はi番目の時刻におけるyの値、t[i]はi番目の時刻、dtは時間刻み幅、f(t[i], y[i])は微分方程式の右辺です。この式に従って、オイラー法をPythonで実装すると以下のようになります。

def euler(f, y0, t):
    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

このコードでは、関数eulerには3つの引数があります。fは微分方程式の右辺を表す関数、y0は初期値、tは時間の配列です。この関数は、数値解を表すyの配列を返します。次に、上で定義した関数を使って微分方程式の数値解を求めます。具体的には、次のコードを実行します。

def f(t, y):
    return t - y

t = np.linspace(0, 5, 1000)
y0 = 1

y_euler = euler(f, y0, t)

このコードでは、f関数には上で定義した微分方程式の右辺を渡しています。tは0から5までの時間の範囲を、1000点に分割した配列です。y0は初期値を表します。最後に、数値解をグラフ化して、厳密解と比較してみましょう。以下のコードを実行して、グラフを描画します。

y_exact = t - 1 + 2*np.exp(-t)

plt.plot(t, y_euler, label='Euler Method')
plt.plot(t, y_exact, label='Exact Solution')
plt.legend()
plt.show()

以上が、Pythonでのオイラー法の実装方法についての具体的な説明です。数値解の精度を向上させるには、オイラー法よりも高次の数値積分法を使用することをお勧めします。

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