数値計算で常微分方程式を解く：4次ルンゲクッタ法の基本と実装

常微分方程式は、解析解が取れないことが多い。だから現場では
「時間を刻んで近似する」が基本になる。

ルンゲクッタ法（特に4次）は、精度と手軽さのバランスが良い「標準の型」。
まずはここから始めれば、一次遅れ・タンク収支・反応・熱収支まで一気に扱える。

ポイントは3つだけ：①状態の置き方 ②刻み幅 $\Delta t$ ③検算（半分刻み）。

• 式は書けたのに、計算が回らない（状態が定義できていない）
• $\Delta t$ を小さくするほど良いと思って、計算が重くなる
• 結果がそれっぽく見えて、検算せずに採用してしまう

常微分方程式は、基本この形にします。

たとえば一次遅れ（入力 $u$ に追従するモデル）は：

タンクの物質収支（濃度 $C$）は：

ルンゲクッタ法（4次）は、「次の値」をこのルールで作ります（スカラーでもベクトルでも同じ）。

最小コード（Python例）

def runge_kutta_4_step(f, t, y, dt):
    k1 = f(t, y)
    k2 = f(t + dt/2, y + dt*k1/2)
    k3 = f(t + dt/2, y + dt*k2/2)
    k4 = f(t + dt,   y + dt*k3)
    return y + dt*(k1 + 2*k2 + 2*k3 + k4)/6

# 例：一次遅れ dy/dt = (u - y)/tau
def f(t, y, u=1.0, tau=5.0):
    return (u - y)/tau

※ 現場用途では「単位」「初期値」「入力 $u$ の持ち方（定数／ステップ／時系列）」を先に揃えると迷いが減ります。

数値計算メモ（コピペ用）

① 目的（何を見る？）
- 対象：一次遅れ／タンク収支／反応／熱収支（など）
- 見たいもの：立ち上がり、遅れ、偏り、オーバーシュート、定常値

② 状態（y）と式（f）
- 状態 y：＿＿＿＿（例：温度、濃度、圧力、レベル）
- 常微分方程式：dy/dt = f(t, y)
- パラメータ：＿＿＿＿（例：τ, F, V, k など）
- 入力/外乱の持ち方：定数／ステップ／時系列

③ 計算条件
- 初期値：y(t0)=＿＿＿＿
- 時間範囲：t0=＿＿ 〜 t_end=＿＿
- 刻み幅：Δt=＿＿（まずは粗めでOK）
- 検算：Δt と Δt/2 を両方回して差を見る

④ 出力（チェック項目）
- y(t) が物理的にあり得る範囲か（負にならない等）
- Δt と Δt/2 の結果差が十分小さいか
- パラメータを振ったとき、直感と逆を向かないか

$\Delta t$ を決めずに悩み続ける。
最初は「粗く回す → 検算する → 必要なら細かくする」で十分です。
いきなり最適な刻み幅を探しに行くと、手が止まります。

単位・初期値・入力の扱いが曖昧なまま回す。
数値計算は正直で、曖昧さがそのまま波形に出ます。式より先に、条件を揃えるのが早道です。