二階常微分方程式ソルバー

二階常微分方程式ソルバー

二階常微分方程式ソルバーは、定数実係数を持つ a·y″ + b·y′ + c·y = g(x) の形式の線形常微分方程式を入力として受け取り、自動的に特性方程式を導出し、減衰条件（過減衰、臨界減衰、不足減衰、無減衰、または不安定）を分類し、記号的な閉形式解と高精度の数値解の両方を作成します。インタラクティブな出力には、y(x) と y′(x) の二重曲線時間プロットと (y, y′) の位相面軌道が組み合わされており、不足減衰の巻き込み、過減衰のノードへの接近、無減衰の閉ループ、不安定の巻き出しなど、挙動を一目で把握できます。

定数係数の二階線形常微分方程式とは何ですか？

実定数係数を持つ二階線形常微分方程式とは、次のような形式の方程式です。

ここで a ≠ 0 であり、b, c は実定数、g(x) は強制項です。2つの初期条件 y(x₀) = y₀ と y′(x₀) = y′₀ を与えることで、これは x₀ の近傍で一意な解を持つ初期値問題となります。これは、等価な一階系に適用されるピカール・リンデレフの定理から導かれます。

g(x) = 0 の場合、方程式は斉次（homogeneous）です。それ以外の場合は非斉次（non-homogeneous）であり、完全な解は次のように分解されます。

ここで y_h は関連する斉次方程式の一般解（2つの任意定数を含む）であり、y_p は元の完全な方程式の任意の特殊解です。2つの初期条件を適用することで、2つの任意定数が固定されます。

特性方程式

斉次方程式において y = e^(r·x) と推測すると、特性方程式（または補助方程式）が得られます。

これは 2 次方程式であり、その判別式 Δ = b² − 4ac が全体の定性的な挙動を制御します。

3つの根のケースと減衰条件

未定係数法（非斉次ケース）

g(x) が以下の単純な形式のいずれかをとる場合、未定係数法により、同じ形式の試行関数を仮定して未知の係数を解くことで特殊解が得られます。

• 定数 g(x) = k。試行関数：y_p = K。c = 0 の場合は x を掛け、さらに b = 0 であれば再び x を掛けます。
• n次の多項式。試行関数：一般的な n 次の多項式。定数項または一次項が共振する場合は、x または x² を掛けます。
• 指数関数 g(x) = A·e^(k·x)。試行関数：y_p = K·e^(k·x)。k が特性根と一致する場合、x（単根）または x²（重根）を掛けます。これを共振と呼びます。
• 正弦波/余弦波 g(x) = A·cos(ω·x) + B·sin(ω·x)。試行関数：y_p = K₁·cos(ω·x) + K₂·sin(ω·x)。iω が特性根である場合は x を掛けます（純粋周波数共振）。
• 積と和は、線形性と積の微分の法則に従います。

位相面の読み方

等価な一階系は u = y, v = y′ であり、u′ = v および v′ = (g(x) − b·v − c·u)/a です。x をパラメータとして u に対して v をプロットすると、位相面軌道が得られます。斉次の自励系（g に x が含まれない）の場合、軌道は開始点 (y₀, y′₀) によって一意に決まり、挙動が一目でわかります。

• 不足減衰：軌道は原点に向かって螺旋状に巻き込まれます。
• 過減衰：軌道は不変直線（遅い固有ベクトル）に沿って原点に接近します。
• 臨界減衰：退化したノードとなり、軌道は単一の固有ベクトルに接します。
• 無減衰：原点を囲む閉じた楕円 — 永久的な振動。
• 不安定：軌道は螺旋を描くか、無限遠に向かって外側に逃げます。

計算例：強制減衰調和振動子

方程式 y″ + 2·y′ + 5·y = 10 で y(0) = 0, y′(0) = 0（強制不足減衰系）を考えます。

1. 特性方程式: r² + 2r + 5 = 0 → Δ = 4 − 20 = −16 → r = −1 ± 2i。
2. 斉次解: y_h = e^(−x)·(C₁·cos 2x + C₂·sin 2x)。
3. 特殊解: 定数強制項 g = 10 に対して y_p = K と置くと、5K = 10 より y_p = 2。
4. 初期条件の適用: y(0) = 0 → C₁ + 2 = 0 → C₁ = −2。y′(0) = 0 → −C₁ + 2C₂ = 0 → C₂ = −1。
5. 最終回答: y(x) = 2 − e^(−x)·(2·cos 2x + sin 2x) — 減衰しながら振動し、極限 y → 2 に収束します。

この電卓の使い方

1. 係数 a, b, c を入力する。最上段に入力します。a は 0 以外である必要があります（0 の場合、方程式は一階となります）。
2. 強制項 g(x) を入力する。斉次問題の場合は 0 のままにします。定数、2次までの多項式、および共振ケースを含む単一の指数関数 A·e^(k·x) に対して閉形式の特殊解が導出されます。
3. 初期条件を入力する (x₀, y₀, y′₀)。二階方程式であるため、x₀ における y と y′ の両方を指定する必要があります。
4. プロット範囲を選択する。ソルバーは RK4 を使用して x₀ から両方の x 方向に積分します。
5. 「解く ＆ 可視化」をクリックする。複素平面上の根を含む特性方程式、減衰条件の分類、斉次および特殊解の閉形式、y と y′ の二重曲線時間プロット、および位相面軌道が表示されます。

主な用途

• 機械的なばね-質量-ダンパ系: m·x″ + c·x′ + k·x = F(t)。過減衰、臨界減衰、不足減衰は、減衰比 ζ = c/(2·√(m·k)) に対応します。
• RLC電気回路: 直列 RLC 回路は L·Q″ + R·Q′ + Q/C = V(t) に従います — 記号は異なりますが構造は同じです。
• 振り子（微小角）: θ″ + (g/L)·θ = 0 は単振動を与え、空気抵抗を加えると減衰振動になります。
• 地震に対する建物の応答: 地盤加速度を強制項とする一自由度構造物。
• PID制御サーボシステム: 閉ループ誤差動特性は二階常微分方程式に集約され、その減衰比がオーバーシュートを左右します。
• 慣性を伴う人口モデル: 資本蓄積の遅れを伴う経済成長、あるいは応答の遅れがある生態学的モデル。

数値解法 — 2次元システムに対する古典的ルンゲ＝クッタ法（RK4）

ツールは a·y″ + b·y′ + c·y = g(x) を次の一階系に縮小します。

ここで u(x₀) = y₀, v(x₀) = y′₀ です。その後、ベクトル状態 (u, v) に対して 4 段のルンゲ＝クッタ法が適用されます。RK4 の局所打切誤差は O(h⁵)、大域誤差は O(h⁴) です。各方向のデフォルト 400 サブステップにより、非スティフな問題に対して概ね 6 桁の精度が得られます。

よくある質問

定数係数の二階線形常微分方程式とは何ですか？

定数係数の二階線形常微分方程式は a·y″ + b·y′ + c·y = g(x) の形式を持ちます。ここで a, b, c は実定数、g(x) は強制（非斉次）項です。2つの初期条件 y(x₀) = y₀ と y′(x₀) = y′₀ により、解は一意に定まります。斉次の場合（g(x) = 0）は常に特性方程式 a·r² + b·r + c = 0 を通じて閉形式の解が得られ、非斉次の場合は y(x) = y_h(x) + y_p(x) として解かれます。

特性方程式とは何ですか？

a·y″ + b·y′ + c·y = 0 に対して y = e^(r·x) を代入すると、a·r² + b·r + c = 0 という特性方程式（補助方程式）が得られます。その根によって斉次解の形式が決まります：異なる2つの実数根は y_h = C₁·e^(r₁·x) + C₂·e^(r₂·x)、重根 r は y_h = (C₁ + C₂·x)·e^(r·x)、共役複素数根 α ± β·i は y_h = e^(α·x)·(C₁·cos(β·x) + C₂·sin(β·x)) となります。

不足減衰、臨界減衰、過減衰とはどういう意味ですか？

これらの用語は、ばね-質量-ダンパモデル m·x″ + c·x′ + k·x = 0 に由来します。過減衰（判別式 > 0、2つの実数根）は、振動せずにゆっくりと平衡状態に戻ることを意味します。臨界減衰（判別式 = 0、重根）は、オーバーシュートせずに最も早く戻る状態です。不足減衰（判別式 < 0、複素数根）は、減衰しながら振動します。無減衰（b = 0, c/a > 0）は、永遠に純粋な正弦波振動を続けます。

未定係数法とは何ですか？

単純な強制項 g(x)（定数、多項式、指数関数、正弦波、余弦波、およびそれらの積）に対して、特殊解 y_p が g と同じ形式を持つと仮定し、方程式に代入して係数を決定する方法です。g(x) が特性根と共振する場合、試行関数に x（または二重根の場合は x²）を掛ける必要があります。

位相面とは何ですか？

二階方程式を2次元システム (y, y′) に縮小し、x の進行に伴う y′ 対 y をプロットしたものです。位相面の解曲線は、不足減衰の減衰螺旋、過減衰の内部ノード、無減衰調和運動の閉じた楕円、不安定振動の外部螺旋など、系の挙動を一目で明らかにします。これは特性方程式の根の図の幾何学的な対応物です。

このツールはどのような数値解法を使用していますか？

等価な一階系 u = y, v = y′（u′ = v, v′ = (g(x) − b·v − c·u)/a）に対して、古典的な四次のルンゲ＝クッタ法（RK4）が適用されます。RK4の局所打切誤差は O(h⁵) であり、デフォルトの方向あたり400サブステップにより、選択したウィンドウ内の非スティフな方程式に対して、概ね6桁の精度が得られます。

参考文献

• 線形微分方程式 — Wikipedia
• 特性方程式 — Wikipedia
• 未定係数法 — Wikipedia (英語)
• 調和振動子 — Wikipedia
• 位相平面 — Wikipedia
• ルンゲ＝クッタ法 — Wikipedia