› Runge-Kutta法 安定性シミュレーター 戻る
数値解析

Runge-Kutta法 安定性シミュレーター

常微分方程式を解く数値積分法（オイラー法・RK2・RK4）の「絶対安定性」を可視化するツールです。刻み幅 h と固有値 λ を変えると、増幅率 R(z) と複素平面の安定領域、そして数値解が発散するかどうかがリアルタイムで分かります。

パラメータ設定
数値積分法
安定関数 R(z) の次数を決める
テスト方程式の λ
y'=λy の固有値。λ<0 で真の解は減衰
刻み幅 h
1ステップで進める時間幅
計算時間 t_end
s
積分を行う総時間
計算結果
—
z = hλ
—
|R(z)| 増幅率
—
安定判定
—
数値解の最終値
—
厳密解の最終値
—
最大誤差
—
複素平面の絶対安定領域 |R(z)| ≤ 1

横軸が z の実部、縦軸が虚部。塗りつぶした領域が選択した手法の安定領域です。現在の動作点 z=hλ（実軸上）は緑＝安定／赤＝不安定で表示されます。

数値解 vs 厳密解 y(t)
増幅率 |R(z)| と刻み幅 h の関係
理論・主要公式

$$y'=\lambda y,\qquad y_{n+1}=R(z)\,y_n,\quad z=h\lambda$$

テスト方程式 y'=λy（厳密解 y(t)=y₀e^(λt)）。1ステップの数値解は安定関数 R(z) を掛けるだけで進み、n ステップ後は y_n = R(z)ⁿ·y₀ となる。

$$R_{\text{Euler}}=1+z,\qquad R_{\text{RK4}}=1+z+\tfrac{z^{2}}{2}+\tfrac{z^{3}}{6}+\tfrac{z^{4}}{24}$$

各手法の安定関数。RK2 は $R=1+z+z^{2}/2$。次数が上がるほど指数関数 eᶻ に近づき、安定領域も広がる。

$$|R(z)|\le 1 \;\Longleftrightarrow\; \text{絶対安定}$$

絶対安定の条件。これを満たす z の集合が「絶対安定領域」。減衰系（λ<0）では z=hλ がこの領域に入るよう h を選ぶ。

Runge-Kutta法の安定性とは

🙋
微分方程式をコンピュータで解くとき、刻み幅 h を小さくすればするほど正確になるって習いました。じゃあ h を適当に小さくしておけば安心ですよね？
🎓
精度の話としてはその通り。でも「精度」と「安定性」は別物なんだ。安定性というのは、h が大きすぎると数値解がどんどん発散して、真の答えとはまったく違う方向にぶっ飛んでいく現象のこと。例えば真の解が静かに減衰してゼロに向かう問題なのに、数値解だけが振動しながら無限大に発散することがある。h を小さくすれば確かに直るけど、「どこまで大きくしても大丈夫か」を知らないと無駄に計算が重くなるんだ。
🙋
発散するかどうかって、何で決まるんですか？左のスライダーで λ をすごくマイナスにしたら、急に「不安定」って赤くなりました。
🎓
いいところに気づいたね。鍵は z = hλ という1つの無次元量だ。テスト方程式 y'=λy を1ステップ解くと、数値解は必ず y_{n+1} = R(z)·y_n という形になる。この R(z) を「安定関数」と呼ぶ。n ステップ進めれば y_n = R(z)ⁿ。だから |R(z)| が 1 より大きいと、ステップを重ねるたびに値が膨らんで発散する。|R(z)| ≤ 1 なら減衰する。λ をマイナスに振ると z=hλ が安定領域の外に出て、それで赤くなったんだよ。
🙋
なるほど！じゃあオイラー法とRK4で安定領域が違うって、上の複素平面の図で見えるやつですか？
🎓
そう、それだ。前進オイラー法は R(z)=1+z だから、安定領域は複素平面で中心 -1・半径1の小さな単位円。実軸でいうと z∈[-2,0] しか安定じゃない。RK4 は R(z)=1+z+z²/2+z³/6+z⁴/24 で、もっと大きな葉っぱ型の領域になって、実軸では z≈-2.78 まで持つ。だから同じ λ の問題でも、RK4 のほうが刻み幅 h を大きく取れる。手法を上手に選ぶこと自体が、計算コストの削減になるんだ。
🙋
それなら、いつでもRK4で h を最大まで攻めればお得じゃないですか？
🎓
気持ちは分かるけど、そこが落とし穴。安定領域の端ギリギリだと、確かに発散はしないけど精度はガタ落ちになる。|R(z)| が 1 に近いと、本当は速く減衰すべき解がいつまでも残ってしまうからね。安定性は『発散しない最低条件』、精度は『どれだけ正解に近いか』。実務では安定限界の半分くらいの h にして、精度に余裕を持たせることが多いよ。特に剛性方程式——速い成分と遅い成分が混じった系——では、陽解法だとこの制約がきつくて、後退オイラー法のような陰解法に切り替えるのが定石なんだ。
🙋
剛性方程式って、具体的にはどんな場面で出てくるんですか？
🎓
化学反応シミュレーションが典型例だね。速い反応はミリ秒で平衡に達するのに、遅い反応は何時間もかかる。両方を1つの系で解こうとすると、速い成分の安定条件に引っ張られて h を極端に小さくせざるを得ない。回路解析、燃焼、生体内の薬物動態なんかも同じ。このツールで λ を -30 まで振ってみると、陽解法だと h を 0.07 くらいまで下げないと安定しないのが分かる。これが剛性の厳しさだよ。

よくある質問

絶対安定性とは、テスト方程式 y'=λy（λ<0 で真の解は減衰）を数値的に解いたとき、数値解が発散せずに減衰し続ける性質のことです。1ステップあたりの増幅率を安定関数 R(z)（z=hλ）とすると、数値解は n ステップ後に y_n = R(z)^n となります。|R(z)| ≤ 1 のときだけステップを重ねても解が大きくならず、絶対安定です。|R(z)| > 1 になると、真の解が減衰しているのに数値解だけが指数的に発散してしまいます。
実軸上で見ると、前進オイラー法 R(z)=1+z は z∈[-2,0] でのみ安定（複素平面では中心 -1・半径1の単位円）です。RK4 R(z)=1+z+z²/2+z³/6+z⁴/24 は実軸上でおよそ z∈[-2.78,0] まで安定で、安定領域はより大きな葉状になります。つまり同じ減衰率 λ の問題でも、RK4 のほうが大きな刻み幅 h を取れます。h を 2/|λ| 程度まで上げてオイラー法では発散する条件でも、RK4 なら 2.78/|λ| までは持ちこたえます。
安定性の観点では、z=hλ が選んだ手法の絶対安定領域の中に入るように h を決めます。減衰する系（λ<0）なら、オイラー法は h < 2/|λ|、RK4 は h < 約2.78/|λ| が安定の上限です。ただし安定なだけでは精度は保証されないため、実務では安定限界よりもさらに小さい h を、要求精度（局所打ち切り誤差）から決めます。安定性は『発散しない最低条件』、精度は『どれだけ正しいか』で、両者は別物です。
剛性方程式とは、非常に速く減衰する成分（|λ| が大きい固有値）と、ゆっくり変化する成分が混在する系です。解の見た目はゆっくりでも、速い成分を安定に解くために h を 2/|λ_max| 以下に抑える必要があり、計算が非効率になります。オイラー法やRK4のような陽解法は安定領域が有限なので剛性問題に弱く、後退オイラー法など A-安定な陰解法（安定領域が左半平面全体）を使うのが定石です。本ツールで λ を -30 など大きくすると、陽解法の限界が体感できます。

実世界での応用

構造振動・過渡応答解析：建物や機械の動的応答を時間積分で求めるとき、固有振動数の高いモードほど |λ| が大きく、陽解法の安定条件 h<2/|λ| を厳しく要求します。メッシュを細かくすると最高固有振動数が上がり、必要な刻み幅も小さくなります。陽解法（中心差分など）を使う衝突・爆発解析では、この安定限界が「クーラン条件」として現れ、要素サイズが時間刻みを決めてしまいます。

化学反応・燃焼シミュレーション：反応速度定数が桁違いに異なる多数の素反応を同時に解く問題は、典型的な剛性方程式です。速い反応の時定数に合わせて h を小さくすると、遅い反応を追うのに膨大なステップ数が必要になります。実務ではRK4のような陽解法ではなく、後退微分公式（BDF）などの陰解法ソルバを使い、安定領域の制約から解放するのが標準です。

電気回路・制御系の数値解析：RC回路やオペアンプ回路の過渡解析では、時定数の小さい寄生容量が大きな |λ| を生みます。SPICE系シミュレータが台形法や後退オイラー法を採用しているのは、まさにこの剛性に対する安定性を確保するためです。制御系の状態方程式を離散化する際も、サンプリング周期 h と固有値の積 z=hλ が安定領域に入るかが設計の前提になります。

CAEソルバの設定確認：陽解法の動的解析で結果が突然発散したとき、多くの場合は時間刻みが安定限界を超えています。本ツールのような安定領域の理解があれば、「メッシュを細かくしたら発散した」「材料を硬くしたら発散した」という現象が、固有値 λ が大きくなって z=hλ が領域外に出たためだと診断できます。逆に陰解法に切り替えれば、安定性のためだけに h を小さくする必要はなくなります。

よくある誤解と注意点

まず最大の誤解が、「高次の手法ほど常に安定」という思い込みです。RK4 はオイラー法より精度が高く、安定領域も実軸方向には広い（-2 対 -2.78）ですが、「高次だから無条件に安定」ではありません。陽的Runge-Kutta法はどれだけ次数を上げても安定領域は有限のままで、左半平面全体をカバーすることは決してありません。剛性方程式に対しては、次数を上げても本質的な解決にならず、A-安定な陰解法に切り替える必要があります。「RK4にしたのに発散する」という相談の多くは、この点の誤解から来ています。

次に、「安定なら正しい答え」だと考えてしまうこと。|R(z)| ≤ 1 はあくまで「数値解が発散しない」という条件であって、精度とは無関係です。例えば z が安定領域の端付近にあると、|R(z)| が 1 に近く、本来速やかに減衰すべき成分がいつまでも残ったり、振動的に減衰したりします。安定でも厳密解とは似ても似つかない解になることは珍しくありません。安定性チェックに合格しても、必ず刻み幅を半分にして解が変わらないか（収束したか）を別途確認してください。

最後に、「テスト方程式 y'=λy は単純すぎて実問題には関係ない」という誤解です。一見おもちゃのような線形スカラー方程式ですが、実際の連立非線形ODEを平衡点まわりで線形化すると、ヤコビ行列の固有値それぞれが λ の役割を果たします。つまり安定性解析は「ヤコビ行列の全固有値 λ_i について z=hλ_i が安定領域に入るか」という形で、そのまま実問題に適用されます。最も絶対値の大きい固有値が刻み幅を支配する——これが剛性の本質であり、テスト方程式の安定性理論が実務で役立つ理由です。

使い方ガイド

  1. 固有値λ（lamNum/lamRange）を-5から+5の範囲で設定し、刻み幅h（hNum/hRange）を0.01から0.5の範囲で入力します
  2. シミュレーターが複素平面z=hλを計算し、各Runge-Kutta法（オイラー法・RK2・RK4）の増幅率|R(z)|を表示します
  3. 安定領域内（|R(z)|≤1）か外かを判定し、常微分方程式dy/dt=λyの数値解と厳密解exp(λt)を比較して最大誤差（teNum/teRange）を確認します

具体的な計算例

λ=-2.0、h=0.1、積分時間t=1.0の場合、z=hλ=-0.2です。RK4法ではz=-0.2が安定領域内（複素平面の左半平面に位置）なので|R(z)|≈0.98となり、厳密解e^(-2.0)≈0.135に対して数値解は0.134で誤差は0.001以下に収束します。一方オイラー法では|R(z)|=|1-0.2|=0.8で同じ結果になりますが、λ=-3.0、h=0.15でz=-0.45の場合、オイラー法は不安定（|R(z)|=1.45>1）となり解が発散しますが、RK4法は|R(z)|≈0.997で安定を保ちます。

実務での注意点

  1. 刚性問題（stiff problem）で負の大きな固有値を含む場合、オイラー法の安定領域は実軸に沿って-2までに限定されるため、陰的Runge-Kutta法の導入が必要です
  2. 複素固有値λ=α±βiを持つ連立系では、z=h(α±βi)が複素平面の安定領域内に完全に収まることを確認し、刻み幅hを制限する必要があります
  3. RK4法は安定領域が広いため実務では推奨されますが、計算コスト（1ステップ4評価）を考慮し、低精度要件ではRK2（2評価）の検討も行います