過渡現象を通して学ぶSimulinkソルバーの特性（Python実装付き）

シミュレーション技術の進展により、制御工学や物理モデリングの分野では、モデルを数値的に解析するための**ソルバー（数値解法アルゴリズム）の選定が重要な役割を果たします。中でもSimulink®に代表される動的システムシミュレータでは、モデルの特性に応じて剛性、非剛性、離散、微分代数方程式（DAE）**などのタイプに適したソルバーを選ぶ必要があります。

本記事では、**一次遅れ系（first-order lag system）**を題材とし、Simulinkに搭載されている主要なソルバー（例：ode45, ode15s, ode23t, FixedStepDiscreteなど）の挙動を、Pythonによる実装を通して模倣・可視化します。

# プログラム名: first_order_all_solvers_and_variants.py
# 概要: 一次遅れ系のステップ応答を複数の数値ソルバーとパラメータで比較
# Program to simulate and compare step responses of a first-order lag system
# using various numerical solvers and parameter sets (gain and time constant)

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# --- 一次遅れ系の微分方程式 / First-order lag differential equation ---
def first_order_ode(t, y, K, tau):
    # dy/dt = (-y + K) / tau
    return (-y + K) / tau

# --- パラメータセット定義 / Define multiple parameter sets ---
parameter_sets = [
    {"K": 1.0, "tau": 2.0, "label": "Base"},      # 基本条件 / Base condition
    {"K": 2.0, "tau": 1.0, "label": "Variant"}    # バリアント条件 / Variant condition
]

# --- 時間設定 / Common time settings ---
y0 = [0.0]                               # 初期値 / Initial state
t_span = (0, 10)                         # 時間範囲 / Simulation time range
t_eval = np.linspace(*t_span, 500)      # 評価時刻 / Evaluation time points

# --- 使用するソルバー一覧 / List of solver methods ---
solver_methods = {
    "ode45 (RK45, Non-stiff)": 'RK45',
    "ode15s (BDF, Stiff)": 'BDF',
    "ode23t (Radau, DAE-like)": 'Radau',
    "odeX (LSODA, auto-switch)": 'LSODA',
    "odeX (DOP853, High-order)": 'DOP853',
    "odeX (RK23, Low-order)": 'RK23'
}

# --- 結果保存用ディクショナリ / Dictionary to store all results ---
all_solutions = {}

# --- 各条件ごとにソルバーを適用 / Run solvers for each parameter set ---
for params in parameter_sets:
    K = params["K"]
    tau = params["tau"]
    label_suffix = f"[K={K}, τ={tau}]"

    # 可変ステップソルバーでのシミュレーション / Variable-step solvers
    for label, method in solver_methods.items():
        sol = solve_ivp(first_order_ode, t_span, y0, method=method, t_eval=t_eval, args=(K, tau))
        all_solutions[f"{label} {label_suffix}"] = sol.y[0]

    # 固定ステップEuler法（ode3の模倣） / Fixed-step Euler method (ode3-like)
    dt = t_eval[1] - t_eval[0]
    y_euler = [y0[0]]
    for i in range(1, len(t_eval)):
        y_prev = y_euler[-1]
        dy = first_order_ode(t_eval[i-1], y_prev, K, tau) * dt
        y_euler.append(y_prev + dy)
    all_solutions[f"ode3 (Euler Fixed-Step) {label_suffix}"] = np.array(y_euler)

# --- グラフ描画 / Plotting all solver responses ---
plt.figure(figsize=(14, 8))

for label, y in all_solutions.items():
    linestyle = '--' if "Variant" in label else '-'
    plt.plot(t_eval, y, linestyle=linestyle, label=label)

# --- 補助線 / Reference lines ---
plt.axhline(1.0, color='gray', linestyle='--', label='Steady State K=1.0')
plt.axhline(2.0, color='gray', linestyle=':', label='Steady State K=2.0')
plt.axvline(2.0, color='red', linestyle='--', label='Time Constant τ=2.0')
plt.axvline(1.0, color='blue', linestyle=':', label='Time Constant τ=1.0')

# --- 図の装飾 / Final figure settings ---
plt.title('First-Order System Step Response using All Solvers and Parameter Variants')
plt.xlabel('Time [s]')
plt.ylabel('Output')
plt.grid(True)
plt.legend(loc='best', fontsize=9)
plt.tight_layout()
plt.show()

結果