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

常微分方程式は、解析解がきれいに求まらない場面が多くあります。だから実務では、
時間を小さく区切って、次の状態を近似的に更新していくやり方が基本になります。

その中でもルンゲクッタ法（4次）は、精度と扱いやすさの釣り合いがよい標準的な方法です。
一次遅れ、タンク収支、反応速度式、熱収支のようなモデルに広く使えます。

最初に意識するべきなのは、①何を状態として置くか ②刻み幅 $\Delta t$ をどう決めるか ③半分刻みで検算するかの3点です。
この3つが揃うと、式が「読めるもの」から「回せるもの」に変わります。

• 式は書けたのに、何を状態として更新するのかが曖昧で、計算手順に落ちない
• $\Delta t$ を小さくすれば安心だと思い込み、必要以上に計算を重くする
• 波形がそれらしく見えたことで満足して、半分刻みの検算をしないまま採用する

数値計算では、まず対象を
「時間に対して何がどう変化するか」
という形に置きます。基本形は次のとおりです。

ここで $y$ は状態です。温度、濃度、液面、圧力のように、
今の時刻における系の持ち物を表します。
$f(t,y)$ は、その状態が今この瞬間にどちら向きへ、どれくらいの速さで変わろうとしているかを与える関数です。

この形に落ちていれば、数値計算は「傾きを読んで、少し先の状態を作る」という作業に整理できます。
逆に言えば、式が複雑でも、状態と傾きの関係さえ定まれば、計算手順自体は共通化できます。

たとえば一次遅れのモデルなら、

となります。入力 $u$ に対して、状態 $y$ が時定数 $\tau$ で遅れて追従する形です。
差 $(u-y)$ が大きいほど変化は速く、追いついてくると変化は小さくなります。

タンクの物質収支なら、濃度 $C$ を状態として

と書けます。流量 $F$、体積 $V$、入口濃度 $C_{\mathrm{in}}$ が与えられれば、
今の濃度 $C$ がどの向きにどれだけ変わるかが決まります。
ここでも本質は同じで、状態と、その変化率の対応を作っているだけです。

最も単純な前進オイラー法では、今の時刻の傾きだけを使って1歩進みます。
ただしそれだけだと、傾きの変化が大きい場面では誤差が出やすくなります。

ルンゲクッタ法（4次）は、その1歩の中で途中の傾きも確認し、
1点の情報ではなく、区間全体の傾きの様子を反映して次の値を作る方法です。
スカラーでもベクトルでも、考え方は同じです。

$k_1$ は今いる場所での傾き、$k_2$ と $k_3$ は半歩先を仮に見に行ったときの傾き、
$k_4$ は1歩先まで行ってみたときの傾きです。
それらを重み付きで平均することで、単純に1方向へ突っ走るよりも、区間内の変化を丁寧に拾えます。

直感的には、「この1歩のあいだに傾きがどれくらい変わりそうかを途中で観察してから進む」
という方法です。だから、前進オイラー法よりも少ない刻みで十分な精度が得られることが多く、
実務での最初の選択肢になりやすいわけです。

最小コード（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

※ 実際に回すときは、単位、初期値、入力の持ち方（定数・ステップ・時系列）を先に揃えると、後からの混乱が減ります。

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

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

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

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

④ 出力の確認
- y(t) は物理的にあり得る範囲か
- Δt と Δt/2 の差は十分小さいか
- パラメータを変えたとき、挙動は直感と大きく逆転しないか

このテンプレの役割は、計算の前に条件を固定することです。
数値計算で迷いやすいのは、解法そのものよりも、
何を入力し、何を初期値とし、何を比較してよしとするかが曖昧なときです。

先にこの枠を埋めておくと、途中で「何を見たかったんだっけ」と戻りにくくなります。
ルンゲクッタ法はあくまで更新の仕組みなので、結果の良し悪しは、条件設定と検算のほうで決まる部分が大きいです。

$\Delta t$ を決める前に、最適解を探し始めてしまう。
刻み幅は、最初から完璧である必要はありません。
まず粗く回し、半分刻みとの差を見て、必要なら細かくする。
その順番で十分です。いきなり最適な $\Delta t$ を当てにいくと、計算そのものが始まりません。

単位、初期値、入力の扱いが揃っていないまま回してしまう。
数値計算は、式の見た目より条件の整合性に敏感です。
たとえば時定数だけ秒なのに、時間軸を分で刻んでいたり、入口濃度の定義が途中で変わっていたりすると、
計算は回っても意味のある結果にはなりません。

波形が滑らかだから正しいと判断してしまう。
数値解は、見た目がきれいでも誤差を含みます。
だからこそ、$\Delta t$ と $\Delta t/2$ の比較や、物理的な範囲の確認が必要です。
「滑らかに見える」は採用条件ではなく、せいぜい出発点にすぎません。

常微分方程式を数値計算で扱うときに必要なのは、特別な魔法ではありません。
状態を決め、変化率を書き、刻み幅を置き、検算する。
その流れを丁寧に踏めば、式はちゃんと手元で動きます。

ルンゲクッタ法（4次）は、その流れを無理なく回しやすい標準形です。
まずは一次遅れやタンク収支のような身近なモデルで試し、
「状態を更新する感覚」を掴むところから始めるのがいちばん確実です。