アルゴリズム

[Python] 常微分方程式のアルゴリズムを実装する方法

2025-01-18更新日: 2025-01-18

Pythonで常微分方程式(ODE)を解くためには、数値的な手法を用いることが一般的です。

代表的な方法としては、オイラー法、ルンゲ＝クッタ法(特に4次のRK4法)が挙げられます。

Pythonでは、scipyライブラリのsolve_ivp関数を使うことで、これらのアルゴリズムを簡単に実装できます。

solve_ivpは、ODEの初期値問題を解くための汎用的な関数で、様々な数値解法(RK45、BDFなど)を選択可能です。

目次から探す

常微分方程式(ODE)とは
Pythonで常微分方程式を解くための準備
オイラー法による常微分方程式の解法
ルンゲ＝クッタ法(RK4法)による解法
scipyのsolve_ivpを使った解法
応用例：物理シミュレーション
応用例：生物学的モデル
応用例：経済学モデル
まとめ

常微分方程式(ODE)とは

常微分方程式(Ordinary Differential Equation, ODE)は、1つまたは複数の変数に依存する関数の導関数を含む方程式です。

これらの方程式は、物理学、工学、生物学、経済学などの多くの分野で現れる現象をモデル化するために使用されます。

ODEは、時間や空間における変化を記述するための強力なツールであり、特に動的システムの挙動を理解するために重要です。

一般的に、ODEは初期条件や境界条件とともに解かれ、特定の状況下での関数の挙動を予測します。

常微分方程式の解法には、解析的手法と数値的手法があり、後者は特に複雑な方程式に対して有効です。

Pythonで常微分方程式を解くための準備

必要なライブラリのインストール

常微分方程式を解くためには、Pythonの数値計算ライブラリであるscipyが必要です。

以下のコマンドを使用して、scipyをインストールできます。

pip install scipy

scipyライブラリの概要

scipyは、科学技術計算のためのPythonライブラリで、数値計算、最適化、統計、信号処理など、さまざまな機能を提供します。

特に、常微分方程式を解くための関数が豊富に用意されており、数値的な解法を簡単に実装できます。

scipy.integrateモジュールには、ODEを解くための主要な関数が含まれています。

solve_ivp関数の基本的な使い方

solve_ivp関数は、初期値問題を解くための主要な関数です。

この関数は、さまざまな数値解法を使用して常微分方程式を解くことができます。

基本的な使い方は以下の通りです。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# 定義する微分方程式
def model(t, y):
    return -2 * y
# 初期条件
y0 = [1]
t_span = (0, 5)  # 時間の範囲
# ODEを解く
solution = solve_ivp(model, t_span, y0, t_eval=np.linspace(0, 5, 100))
# 結果のプロット
plt.plot(solution.t, solution.y[0])
plt.title('ODE Solution')
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.grid()
plt.show()

このコードでは、単純な微分方程式 $\frac{d y}{d t} = - 2 y$ を解いています。

solve_ivp関数を使用することで、初期条件に基づいて時間の範囲内での解を得ることができます。

出力結果は、時間に対する関数の挙動を示すグラフとして表示されます。

オイラー法による常微分方程式の解法

オイラー法の基本原理

オイラー法は、常微分方程式を数値的に解くための最も基本的な手法の一つです。

この方法は、微分方程式の解を近似するために、関数の接線を用いて次の点を推定します。

具体的には、次のように定義されます。

初期条件 $y (t_{0}) = y_{0}$ を設定します。
微小時間 $h$ を選びます。
次の点を次の式で計算します：

$y (t_{n + 1}) = y (t_{n}) + h \cdot f (t_{n}, y (t_{n}))$

ここで、 $f (t, y)$ は微分方程式の右辺です。

このプロセスを繰り返すことで、時間の経過に伴う関数の値を近似的に求めることができます。

Pythonでのオイラー法の実装

以下は、オイラー法を用いて常微分方程式を解くPythonのサンプルコードです。

import numpy as np
import matplotlib.pyplot as plt
# 定義する微分方程式
def f(t, y):
    return -2 * y
# オイラー法の実装
def euler_method(f, y0, t0, t_end, h):
    n = int((t_end - t0) / h)  # ステップ数
    t_values = np.linspace(t0, t_end, n + 1)
    y_values = np.zeros(n + 1)
    y_values[0] = y0
    for i in range(n):
        y_values[i + 1] = y_values[i] + h * f(t_values[i], y_values[i])
    return t_values, y_values
# 初期条件とパラメータ
y0 = 1
t0 = 0
t_end = 5
h = 0.1
# オイラー法を実行
t_values, y_values = euler_method(f, y0, t0, t_end, h)
# 結果のプロット
plt.plot(t_values, y_values, label='Euler Method', marker='o')
plt.title('Euler Method for ODE')
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.grid()
plt.legend()
plt.show()

このコードでは、微分方程式 $\frac{d y}{d t} = - 2 y$ をオイラー法で解いています。

初期条件と時間の範囲を設定し、オイラー法を実行して得られた結果をプロットしています。

オイラー法のメリットとデメリット

メリット	デメリット
実装が簡単で理解しやすい	精度が低く、大きなステップサイズでは不安定になることがある
計算が高速で、少ないリソースで実行可能	高次の精度が必要な場合には不適切
初期値問題に対して適用しやすい	ステップサイズを小さくしないと精度が保証されない

オイラー法は、シンプルで使いやすい手法ですが、精度や安定性に関しては限界があります。

特に、非線形の問題や急激な変化がある場合には、他の数値解法を検討する必要があります。

ルンゲ＝クッタ法(RK4法)による解法

ルンゲ＝クッタ法の概要

ルンゲ＝クッタ法(Runge-Kutta method)は、常微分方程式を数値的に解くための非常に一般的で強力な手法です。

特に、4次のルンゲ＝クッタ法(RK4法)は、精度と計算コストのバランスが良いため、広く使用されています。

RK4法は、微分方程式の解を求める際に、複数の中間点を計算し、それらを用いて次のステップの値を推定します。

この方法は、オイラー法よりも高い精度を持ち、特に非線形の問題に対しても効果的です。

4次ルンゲ＝クッタ法の数式

RK4法では、次のステップの値を以下の数式を用いて計算します。

初期条件 $y (t_{n})$ から次の点 $y (t_{n + 1})$ を求めるために、次の4つの中間値を計算します。

$k_{1} = h \cdot f (t_{n}, y_{n})$
$k_{2} = h \cdot f (t_{n} + \frac{h}{2}, y_{n} + \frac{k_{1}}{2})$
$k_{3} = h \cdot f (t_{n} + \frac{h}{2}, y_{n} + \frac{k_{2}}{2})$
$k_{4} = h \cdot f (t_{n} + h, y_{n} + k_{3})$

次に、次のステップの値は次のように計算されます：

$y_{n + 1} = y_{n} + \frac{1}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4})$

PythonでのRK4法の実装

以下は、RK4法を用いて常微分方程式を解くPythonのサンプルコードです。

import numpy as np
import matplotlib.pyplot as plt
# 定義する微分方程式
def f(t, y):
    return -2 * y
# ルンゲ＝クッタ法の実装
def rk4_method(f, y0, t0, t_end, h):
    n = int((t_end - t0) / h)  # ステップ数
    t_values = np.linspace(t0, t_end, n + 1)
    y_values = np.zeros(n + 1)
    y_values[0] = y0
    for i in range(n):
        k1 = h * f(t_values[i], y_values[i])
        k2 = h * f(t_values[i] + h / 2, y_values[i] + k1 / 2)
        k3 = h * f(t_values[i] + h / 2, y_values[i] + k2 / 2)
        k4 = h * f(t_values[i] + h, y_values[i] + k3)
        
        y_values[i + 1] = y_values[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6
    return t_values, y_values
# 初期条件とパラメータ
y0 = 1
t0 = 0
t_end = 5
h = 0.1
# RK4法を実行
t_values, y_values = rk4_method(f, y0, t0, t_end, h)
# 結果のプロット
plt.plot(t_values, y_values, label='RK4 Method', marker='o')
plt.title('RK4 Method for ODE')
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.grid()
plt.legend()
plt.show()

このコードでは、微分方程式 $\frac{d y}{d t} = - 2 y$ をRK4法で解いています。

初期条件と時間の範囲を設定し、RK4法を実行して得られた結果をプロットしています。

RK4法の精度と計算コスト

RK4法は、オイラー法や他の低次の手法に比べて高い精度を持っています。

具体的には、RK4法は4次の精度を持ち、ステップサイズを半分にすると誤差が約16分の1に減少します。

しかし、計算コストは高く、各ステップで4回の関数評価が必要です。

これにより、計算量が増加し、特に大規模な問題や長時間のシミュレーションでは計算時間が長くなる可能性があります。

手法	精度	計算コスト
オイラー法	1次	1回の関数評価
ルンゲ＝クッタ法(RK4)	4次	4回の関数評価

RK4法は、精度が求められる問題に対して非常に有効ですが、計算コストを考慮する必要があります。

scipyのsolve_ivpを使った解法

solve_ivpの基本的な使い方

scipy.integrateモジュールのsolve_ivp関数は、初期値問題を解くための強力なツールです。

基本的な使い方は以下の通りです。

微分方程式を定義します。
初期条件と時間の範囲を設定します。
solve_ivp関数を呼び出して解を求めます。

以下は、solve_ivpを使用した基本的な例です。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# 定義する微分方程式
def model(t, y):
    return -2 * y
# 初期条件
y0 = [1]
t_span = (0, 5)  # 時間の範囲
# ODEを解く
solution = solve_ivp(model, t_span, y0, t_eval=np.linspace(0, 5, 100))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='solve_ivp Solution')
plt.title('ODE Solution using solve_ivp')
plt.xlabel('Time')
plt.ylabel('y(t)')
plt.grid()
plt.legend()
plt.show()

このコードでは、微分方程式 $\frac{d y}{d t} = - 2 y$ を解いています。

solve_ivp関数を使用することで、初期条件に基づいて時間の範囲内での解を得ることができます。

solve_ivpでの解法選択(RK45, BDFなど)

solve_ivp関数では、さまざまな数値解法を選択することができます。

主な解法には以下のものがあります。

RK45: ルンゲ＝クッタ法の4次および5次の適応型手法。

一般的に使用される。

RK23: ルンゲ＝クッタ法の2次および3次の適応型手法。

計算コストが低い。

BDF: バックワード差分法。

剛性のある問題に適している。

LSODA: 自動的に非剛性と剛性の問題を切り替える手法。

解法を指定するには、method引数を使用します。

例えば、RK45を使用する場合は次のようにします。

solution = solve_ivp(model, t_span, y0, method='RK45', t_eval=np.linspace(0, 5, 100))

solve_ivpのオプション設定

solve_ivp関数には、さまざまなオプションを設定することができます。

主なオプションには以下のものがあります。

t_eval: 解を評価する時間の配列を指定します。
max_step: 最大ステップサイズを指定します。
rtol: 相対誤差の許容値を指定します。
atol: 絶対誤差の許容値を指定します。

これらのオプションを使用することで、解法の精度や計算の挙動を調整することができます。

実際のODE問題を解く例

以下は、実際のODE問題を解く例です。

ここでは、単振り子の運動方程式を解きます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# 単振り子の運動方程式
def pendulum(t, y):
    g = 9.81  # 重力加速度
    L = 1.0   # 振り子の長さ
    theta, omega = y
    dydt = [omega, -g/L * np.sin(theta)]
    return dydt
# 初期条件
y0 = [np.pi / 4, 0]  # 初期角度と初期角速度
t_span = (0, 10)     # 時間の範囲
# ODEを解く
solution = solve_ivp(pendulum, t_span, y0, method='RK45', t_eval=np.linspace(0, 10, 100))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Angle (theta)')
plt.plot(solution.t, solution.y[1], label='Angular Velocity (omega)')
plt.title('Pendulum Motion using solve_ivp')
plt.xlabel('Time (s)')
plt.ylabel('Values')
plt.grid()
plt.legend()
plt.show()

このコードでは、単振り子の運動方程式を定義し、初期条件を設定してsolve_ivpを使用して解を求めています。

結果は、時間に対する角度と角速度の変化を示すグラフとして表示されます。

応用例：物理シミュレーション

単振り子の運動方程式を解く

単振り子は、物理学における基本的なモデルであり、振り子の運動を理解するための良い例です。

単振り子の運動方程式は次のように表されます：

$\frac{d^{2} θ}{d t^{2}} + \frac{g}{L} \sin (θ) = 0$

ここで、 $θ$ は振り子の角度、 $g$ は重力加速度、 $L$ は振り子の長さです。

この方程式を解くために、solve_ivpを使用して数値的に解くことができます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# 単振り子の運動方程式
def pendulum(t, y):
    g = 9.81  # 重力加速度
    L = 1.0   # 振り子の長さ
    theta, omega = y
    dydt = [omega, -g/L * np.sin(theta)]
    return dydt
# 初期条件
y0 = [np.pi / 4, 0]  # 初期角度と初期角速度
t_span = (0, 10)     # 時間の範囲
# ODEを解く
solution = solve_ivp(pendulum, t_span, y0, method='RK45', t_eval=np.linspace(0, 10, 100))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Angle (theta)')
plt.title('Pendulum Motion')
plt.xlabel('Time (s)')
plt.ylabel('Angle (radians)')
plt.grid()
plt.legend()
plt.show()

このコードでは、単振り子の運動をシミュレーションし、時間に対する角度の変化をプロットしています。

惑星の軌道計算

惑星の運動は、万有引力の法則に基づいて記述されます。

2つの天体間の引力は次のように表されます：

$F = \frac{G m_{1} m_{2}}{r^{2}}$

ここで、 $G$ は万有引力定数、 $m_{1}$ と $m_{2}$ はそれぞれの天体の質量、 $r$ は天体間の距離です。

この力を用いて、惑星の運動方程式を数値的に解くことができます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# 惑星の運動方程式
def planet_motion(t, y):
    G = 6.67430e-11  # 万有引力定数
    m_sun = 1.989e30  # 太陽の質量
    x, y_pos, vx, vy_pos = y
    r = np.sqrt(x**2 + y_pos**2)
    ax = -G * m_sun * x / r**3
    ay = -G * m_sun * y_pos / r**3
    return [vx, vy_pos, ax, ay]
# 初期条件(地球の初期位置と速度)
y0 = [1.496e11, 0, 0, 29780]  # 地球の初期位置(m)と速度(m/s)
t_span = (0, 3.154e7)  # 1年(秒)
# ODEを解く
solution = solve_ivp(planet_motion, t_span, y0, method='RK45', t_eval=np.linspace(0, 3.154e7, 1000))
# 結果のプロット
plt.plot(solution.y[0], solution.y[1], label='Earth Orbit')
plt.title('Planetary Motion')
plt.xlabel('X Position (m)')
plt.ylabel('Y Position (m)')
plt.axis('equal')
plt.grid()
plt.legend()
plt.show()

このコードでは、地球の軌道をシミュレーションし、太陽の周りの軌道をプロットしています。

電気回路の微分方程式モデル

電気回路の動作は、オームの法則やキルヒホッフの法則に基づいて記述されます。

例えば、RLC回路(抵抗、インダクタ、キャパシタからなる回路)の微分方程式は次のように表されます：

$L \frac{d^{2} i}{d t^{2}} + R \frac{d i}{d t} + \frac{1}{C} i = V (t)$

ここで、 $i$ は電流、 $R$ は抵抗、 $L$ はインダクタンス、 $C$ はキャパシタンス、 $V (t)$ は外部電圧です。

この方程式を解くために、solve_ivpを使用して数値的に解くことができます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# RLC回路の微分方程式
def rlc_circuit(t, y):
    R = 1.0  # 抵抗(Ω)
    L = 1.0  # インダクタンス(H)
    C = 1.0  # キャパシタンス(F)
    V = 5.0  # 外部電圧(V)
    i, q = y
    di_dt = (V - R * i - q / C) / L
    dq_dt = i
    return [di_dt, dq_dt]
# 初期条件
y0 = [0, 0]  # 初期電流と初期電荷
t_span = (0, 10)  # 時間の範囲
# ODEを解く
solution = solve_ivp(rlc_circuit, t_span, y0, method='RK45', t_eval=np.linspace(0, 10, 100))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Current (i)')
plt.plot(solution.t, solution.y[1], label='Charge (q)')
plt.title('RLC Circuit Response')
plt.xlabel('Time (s)')
plt.ylabel('Values')
plt.grid()
plt.legend()
plt.show()

このコードでは、RLC回路の応答をシミュレーションし、時間に対する電流と電荷の変化をプロットしています。

これにより、電気回路の動作を視覚的に理解することができます。

応用例：生物学的モデル

ロトカ＝ヴォルテラ方程式(捕食者-被食者モデル)

ロトカ＝ヴォルテラ方程式は、捕食者と被食者の相互作用をモデル化するための微分方程式です。

このモデルは、被食者の個体数 $x$ と捕食者の個体数 $y$ の変化を次のように表現します：

$\frac{d x}{d t} = α x - β x y$

$\frac{d y}{d t} = δ x y - γ y$

ここで、 $α$ は被食者の成長率、 $β$ は捕食者による被食者の捕食率、 $δ$ は捕食者の成長率、 $γ$ は捕食者の死亡率です。

この方程式を数値的に解くことで、捕食者と被食者の動態をシミュレーションできます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# ロトカ＝ヴォルテラ方程式
def lotka_volterra(t, y):
    alpha = 0.1  # 被食者の成長率
    beta = 0.02  # 捕食者による被食者の捕食率
    delta = 0.01  # 捕食者の成長率
    gamma = 0.1  # 捕食者の死亡率
    x, y = y
    dxdt = alpha * x - beta * x * y
    dydt = delta * x * y - gamma * y
    return [dxdt, dydt]
# 初期条件
y0 = [40, 9]  # 被食者と捕食者の初期個体数
t_span = (0, 200)  # 時間の範囲
# ODEを解く
solution = solve_ivp(lotka_volterra, t_span, y0, method='RK45', t_eval=np.linspace(0, 200, 1000))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Prey (被食者)')
plt.plot(solution.t, solution.y[1], label='Predator (捕食者)')
plt.title('Lotka-Volterra Model')
plt.xlabel('Time')
plt.ylabel('Population')
plt.grid()
plt.legend()
plt.show()

このコードでは、捕食者と被食者の個体数の変化をシミュレーションし、時間に対する個体数の変化をプロットしています。

SIRモデルによる感染症の拡大シミュレーション

SIRモデルは、感染症の拡大をモデル化するための基本的なモデルです。

このモデルでは、個体を感受性者(S)、感染者(I)、回復者(R)の3つの状態に分類します。

モデルの方程式は次のように表されます：

$\frac{d S}{d t} = - β S I$

$\frac{d I}{d t} = β S I - γ I$

$\frac{d R}{d t} = γ I$

ここで、 $β$ は感染率、 $γ$ は回復率です。

この方程式を数値的に解くことで、感染症の拡大をシミュレーションできます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# SIRモデル
def sir_model(t, y):
    beta = 0.3  # 感染率
    gamma = 0.1  # 回復率
    S, I, R = y
    dSdt = -beta * S * I
    dIdt = beta * S * I - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]
# 初期条件
N = 1000  # 総人口
I0 = 1    # 初期感染者
R0 = 0    # 初期回復者
S0 = N - I0 - R0  # 初期感受性者
y0 = [S0, I0, R0]
t_span = (0, 160)  # 時間の範囲
# ODEを解く
solution = solve_ivp(sir_model, t_span, y0, method='RK45', t_eval=np.linspace(0, 160, 1000))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Susceptible (感受性者)')
plt.plot(solution.t, solution.y[1], label='Infected (感染者)')
plt.plot(solution.t, solution.y[2], label='Recovered (回復者)')
plt.title('SIR Model of Infectious Disease')
plt.xlabel('Time')
plt.ylabel('Population')
plt.grid()
plt.legend()
plt.show()

このコードでは、SIRモデルを用いて感染症の拡大をシミュレーションし、時間に対する各状態の個体数の変化をプロットしています。

化学反応速度論のモデル

化学反応速度論では、反応物と生成物の濃度の変化を微分方程式で表現します。

例えば、単純な一次反応 $A \to B$ の場合、反応速度は次のように表されます：

$\frac{d [A]}{d t} = - k [A]$

$\frac{d [B]}{d t} = k [A]$

ここで、 $k$ は反応速度定数です。

この方程式を数値的に解くことで、反応物と生成物の濃度の変化をシミュレーションできます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# 化学反応速度論のモデル
def reaction_kinetics(t, y):
    k = 0.1  # 反応速度定数
    A, B = y
    dA_dt = -k * A
    dB_dt = k * A
    return [dA_dt, dB_dt]
# 初期条件
A0 = 1.0  # 初期濃度
B0 = 0.0  # 初期濃度
y0 = [A0, B0]
t_span = (0, 50)  # 時間の範囲
# ODEを解く
solution = solve_ivp(reaction_kinetics, t_span, y0, method='RK45', t_eval=np.linspace(0, 50, 1000))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Concentration of A (反応物)')
plt.plot(solution.t, solution.y[1], label='Concentration of B (生成物)')
plt.title('Chemical Reaction Kinetics')
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.grid()
plt.legend()
plt.show()

このコードでは、化学反応の濃度変化をシミュレーションし、時間に対する反応物と生成物の濃度の変化をプロットしています。

これにより、化学反応の動態を視覚的に理解することができます。

応用例：経済学モデル

ソロー成長モデルのシミュレーション

ソロー成長モデルは、経済成長を説明するための古典的なモデルで、資本の蓄積、労働力の成長、技術進歩が経済成長に与える影響を考慮します。

このモデルでは、次のような微分方程式が用いられます：

$\frac{d k}{d t} = s f (k) - (n + δ) k$

ここで、 $k$ は資本の蓄積、 $s$ は貯蓄率、 $f (k)$ は生産関数、 $n$ は労働力の成長率、 $δ$ は資本の減耗率です。

この方程式を数値的に解くことで、経済成長のシミュレーションが可能です。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# ソロー成長モデル
def solow_growth(t, k):
    s = 0.3  # 貯蓄率
    n = 0.01  # 労働力の成長率
    delta = 0.05  # 資本の減耗率
    alpha = 0.3  # 生産関数の弾力性
    f_k = k**alpha  # 生産関数
    dk_dt = s * f_k - (n + delta) * k
    return [dk_dt]
# 初期条件
k0 = 1.0  # 初期資本
t_span = (0, 100)  # 時間の範囲
# ODEを解く
solution = solve_ivp(solow_growth, t_span, [k0], method='RK45', t_eval=np.linspace(0, 100, 1000))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Capital per Worker (k)')
plt.title('Solow Growth Model Simulation')
plt.xlabel('Time')
plt.ylabel('Capital per Worker')
plt.grid()
plt.legend()
plt.show()

このコードでは、ソロー成長モデルを用いて資本の蓄積をシミュレーションし、時間に対する資本の変化をプロットしています。

ロジスティック成長モデル

ロジスティック成長モデルは、個体群の成長をモデル化するための手法で、環境のキャパシティを考慮します。

このモデルは次のように表されます：

$\frac{d P}{d t} = r P (1 - \frac{P}{K})$

ここで、 $P$ は個体数、 $r$ は成長率、 $K$ は環境のキャパシティです。

この方程式を数値的に解くことで、個体群の成長をシミュレーションできます。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
# ロジスティック成長モデル
def logistic_growth(t, P):
    r = 0.1  # 成長率
    K = 100  # 環境のキャパシティ
    dP_dt = r * P * (1 - P / K)
    return [dP_dt]
# 初期条件
P0 = 10  # 初期個体数
t_span = (0, 100)  # 時間の範囲
# ODEを解く
solution = solve_ivp(logistic_growth, t_span, [P0], method='RK45', t_eval=np.linspace(0, 100, 1000))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='Population (P)')
plt.title('Logistic Growth Model Simulation')
plt.xlabel('Time')
plt.ylabel('Population')
plt.grid()
plt.legend()
plt.show()

このコードでは、ロジスティック成長モデルを用いて個体群の成長をシミュレーションし、時間に対する個体数の変化をプロットしています。

IS-LMモデルの動的シミュレーション

IS-LMモデルは、経済の総需要と総供給の関係を示すためのモデルで、利子率と国民所得の関係を考慮します。

IS曲線は投資と貯蓄の均衡を示し、LM曲線は貨幣市場の均衡を示します。

これらの曲線を用いて、経済の動的な変化をシミュレーションできます。

IS曲線は次のように表されます：

$Y = C (Y - T) + I (r) + G$

LM曲線は次のように表されます：

$M / P = L (Y, r)$

ここで、 $Y$ は国民所得、 $C$ は消費、 $T$ は税金、 $I$ は投資、 $G$ は政府支出、 $M$ は貨幣供給、 $P$ は物価、 $L$ は貨幣需要です。

このモデルを数値的に解くことで、経済の動的な変化をシミュレーションできます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# IS-LMモデル
def is_lm(t, y):
    Y, r = y
    C = 0.8 * (Y - 100)  # 消費関数
    I = 200 - 10 * r  # 投資関数
    G = 100  # 政府支出
    M = 1000  # 貨幣供給
    P = 1  # 物価
    L = 0.5 * Y - 20 * r  # 貨幣需要
    dY_dt = C + I + G - Y  # IS曲線
    dr_dt = (M / P - L) / 10  # LM曲線
    return [dY_dt, dr_dt]
# 初期条件
Y0 = 500  # 初期国民所得
r0 = 5    # 初期利子率
t_span = (0, 100)  # 時間の範囲
# ODEを解く
solution = solve_ivp(is_lm, t_span, [Y0, r0], method='RK45', t_eval=np.linspace(0, 100, 1000))
# 結果のプロット
plt.plot(solution.t, solution.y[0], label='National Income (Y)')
plt.plot(solution.t, solution.y[1], label='Interest Rate (r)')
plt.title('IS-LM Model Simulation')
plt.xlabel('Time')
plt.ylabel('Values')
plt.grid()
plt.legend()
plt.show()