各手法の精度オーダー
ホイン: O(h²) 2次
RK4: O(h⁴) 4次
刻み幅を½にすると
RK4の誤差は1/16倍
オイラー・ホイン・RK2・RK4・RK45法を同一問題で比較。刻み幅による精度・安定性の違いを可視化。ロジスティック成長・ファン・デル・ポール振動子・ローレンツ方程式のプリセット付き。
このシミュレーターで解く一般的な常微分方程式の初期値問題は、次のように表されます。
$$ \frac{dy}{dt}= f(t, y), \quad y(t_0) = y_0 $$ここで、$y$は求めたい未知関数(例:人口、振動の変位)、$t$は時間、$f(t, y)$は変化率を表す既知の関数です。$y_0$は初期状態です。
代表的な解法である4次ルンゲ・クッタ法(RK4)の更新式は次の通りです。1ステップ($t_n$から$t_{n+1}=t_n+h$)進めるために、4種類の傾き$k_1, k_2, k_3, k_4$を計算します。
$$ \begin{aligned}k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1) \\ k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_2) \\ k_4 &= f(t_n + h, y_n + h k_3) \\ y_{n+1}&= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned}$$$h$は刻み幅です。中点や終点での傾きを予測し、それらを重み付き平均することで、オイラー法($y_{n+1}=y_n+h f(t_n, y_n)$)よりもはるかに高精度な解が得られます。この「4次」というのは、局所誤差が$O(h^5)$、大域誤差が$O(h^4)$であることを意味します。
機械・自動車工学:エンジンのピストン運動やサスペンションの振動解析に利用されます。リアルタイム制御が必要なECU(エンジン制御ユニット)内の計算では計算コストの低いオイラー法が、設計段階での詳細な性能評価にはRK4やRK45がよく用いられます。
電気・電子回路シミュレーション:トランジスタやオペアンプを含む非線形回路の過渡応答を計算するために必須です。スイッチングの瞬間など電圧・電流が急激に変化する部分では、RK45のような刻み幅制御機能を持つ手法が威力を発揮します。
化学反応・生物増殖モデリング:複数の化学物質が関与する反応速度論や、人口増加のロジスティック方程式の予測に使われます。反応速度が速い物質と遅い物質が混在する「硬い方程式」では、解法の選択が計算の成否を分けます。
気象・流体シミュレーション(CFD):ナビエ-ストークス方程式を時間方向に離散化する際の時間積分スキームとして応用されます。乱流の計算などでは膨大な計算コストがかかるため、精度と効率を両立させる高度な解法が研究されています。
まず、「高次=常に正義」ではないという点を押さえよう。確かにRK4はオイラー法より圧倒的に精度が高い。しかし、例えば制御システムのリアルタイムシミュレーションのように、計算速度が最優先で、多少の誤差がフィードバックループで吸収される場合は、オイラー法が採用されることもある。このツールで「ロジスティック成長」を選び、刻み幅h=1.0で各手法を比較してみてほしい。RK4とオイラー法の結果は大きく違うが、h=0.1まで細かくすると、オイラー法でも実用に足る解に近づく。つまり、「刻み幅」と「次数」はトレードオフの関係にあるんだ。
次に、「安定性」を見落とさないでほしい。例えば「ファン・デル・ポール振動子」で、刻み幅をh=2.0など非常に大きく設定してオイラー法を実行すると、解が発散してグラフが吹き飛ぶことがある。これは数値的不安定性と呼ばれる現象で、手法によって許容できる刻み幅の上限が違う。実務では、計算が破綻しない安定した刻み幅を事前に見積もることが必須だ。
最後に、RK45の「自動刻み幅制御」を過信しないこと。これは万能ツールではなく、許容誤差(ツール内の「許容誤差」パラメータ)という目標精度をユーザーが設定する必要がある。許容誤差を甘く(例えば1e-3)設定すると計算は速いが粗い解に、厳しく(1e-9)設定すると高精度だが計算コストが跳ね上がる。常に「どの程度の精度が必要か」という要求仕様から逆算してパラメータを決めよう。
このシミュレーターで扱うODE解法は、あらゆる「時間変化する現象」のシミュレーションの根幹を成す。例えば、電気回路の過渡応答解析では、コイルやコンデンサの電圧・電流の時間変化が連立ODEで記述される。スイッチを入れた瞬間の振る舞いをRK4などで解くことで、適切なサージ保護素子を設計できる。
また、化学工学における反応器の設計でも必須だ。連続攪拌槽型反応器(CSTR)内での化学種の濃度変化は非線形ODEでモデル化される。反応速度が速い部分ではRK45で刻み幅を細かく、遅い部分では粗く計算することで、効率的に安定操作条件を探れる。
さらに高度な応用として、構造力学の動的解析(過渡応答解析)がある。自動車が段差を越える時の車体の振動や、風力発電のブレードに働く変動荷重への応答は、巨大な連立ODEシステムとして離散化される。ここでは計算コストが膨大になるため、陰解法など別の解法ファミリーも組み合わせつつ、RK系の手法が部分的なサブシステムのシミュレーションに活用されるんだ。
まず次の一歩は、「連立ODE」と「高階ODE」への挑戦だ。このツールの「ローレンツ方程式」は実は3本の連立ODEだ。多くの実問題は複数の変数が相互に影響し合う連立系になる。また、振動問題でおなじみの $m \frac{d^2x}{dt^2} + c \frac{dx}{dt} + kx = F$ のような2階の方程式は、$v = \frac{dx}{dt}$ と置くことで、$x$ と $v$ に関する2本の連立1階ODEに変換できる。この変換技術をマスターすれば、解ける問題の幅が一気に広がる。
数学的背景を深めたいなら、各解法の「収束性」と「安定性領域」を調べよう。収束性は「刻み幅を無限に小さくすれば真の解に近づくか」を理論保証するもの。安定性領域は、テスト方程式 $\frac{dy}{dt} = \lambda y$ ($\lambda$ は複素数)に対して、数値解が発散しない刻み幅 $h$ の条件を複素平面上に描いた領域だ。オイラー法の安定性領域は小さく、RK4のそれはより広い。この概念は、振動や減衰を含む問題で数値解が暴走する理由を理解する鍵になる。
最終的には、「陰解法」の世界に進むことを勧める。オイラー法やRK4は「陽解法」と呼ばれ、現在の値から明示的に未来の値を計算する。一方、陰解法(例えば後退オイラー法やクランク・ニコルソン法)は方程式を解く(反復計算が必要)ことで未来の値を求める。剛性(スティッフ)な問題、つまり現象に非常に速い時定数と遅い時定数が混在する場合、陽解法では極端に刻み幅を細かくしないと不安定になるが、陰解法では大幅な計算効率向上が可能だ。このツールで学んだ「刻み幅と精度・安定性の関係」のその先にある、実務計算の核心の一つだ。