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

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

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

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

ストーク

シママ、ODEって聞いた瞬間に身構える人が多かたい。
「結局どう計算すればよか？」を一言で言うなら、どげん？

シママ

「時間を小さく区切って、次の値を順番に作る」だけ。
その“次の値の作り方”として、まずRK4を覚える。これで十分スタートできる。

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

ストーク

Δtを小さくしたら勝ち、って思いがちたい。
でも現場は「速さ」も要るけん、加減が難しか。

シママ

そう。だから「検算の型」を最初に持つのがコスパ良い。
同じ条件で、ΔtとΔt/2で回して“差”を見る。差が小さければ、まず合格。

ODEは、基本この形にします。

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

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

ストーク

なるほど。
「式」より先に、まず yを何にするか を決めるのが肝たいね。

シママ

そう。y（状態）が決まると、計算は“作業”になる。
あとは f(t, y) を関数として書ければ、ルンゲクッタに流し込める。

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

ストーク

これ、式だけ見ると厳つかたい。
感覚的には「途中で何回か様子見して、平均して前に進む」…で合っとる？

シママ

合ってる。
オイラー法は“最初の傾き”だけで突っ切るけど、RK4は途中も覗いて、曲がり方を拾う感じ。

最小コード（Python例）

def rk4_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の持ち方（ステップ/時系列）」を先に揃えると迷いが減ります。

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

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

② 状態（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 の結果差が十分小さいか
- パラメータを振ったとき、直感と逆を向かないか

ストーク

「検算＝ΔtとΔt/2」…これが一番効きそうたい。
見た目に騙されんための“保険”になるけん。

シママ

うん。まずそこだけは固定でいい。
精度の議論をする前に、同じモデルが“刻み方”で崩れないか確認する。それだけで事故が減る。

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

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

ストーク

「条件を揃える」…耳が痛か。
式を綺麗にして満足しがちたい。

シママ

式は大事。でも“回って初めて”価値になる。
だから、テンプレに沿って状態と条件を先に固定する。そこから。

ストーク

まとめるばい。ODEは「刻む」、型はRK4、検算はΔt/2。
これで“式が回る”状態を作る、たいね。

シママ

そう。次は「複数状態（ベクトル）」と「入力が時系列のときの扱い」をやれば、実務の幅が一気に広がる。

※ 次の資料案：2槽タンク／反応速度式（Arrhenius）を含む熱収支／安定性（発散しやすい系でのΔt設計）。