2階線形常微分方程式の解法

バネにつないだ重りを引っ張って放すと振動が始まります。空気抵抗がなければ永遠に振動し続け（単振動）、空気抵抗があれば徐々に減衰します。この振動現象を記述するのが2階線形常微分方程式です。

なぜ「2階」なのでしょうか。ニュートンの運動方程式 $F = ma$ において、加速度 $a = d^2x/dt^2$ は位置の2階微分です。力が位置 $x$ や速度 $dx/dt$ に依存する場合、運動方程式は自然に2階の微分方程式になります。したがって、2階ODEは力学系のダイナミクスを記述する最も基本的な数学的枠組みなのです。

2階線形ODEの一般的な形は

$$ \frac{d^2y}{dx^2} + p(x)\frac{dy}{dx} + q(x)y = g(x) $$

です。左辺は $y$ とその導関数の線形結合であり、$g(x) = 0$ の場合を斉次（homogeneous）、$g(x) \neq 0$ の場合を非斉次（inhomogeneous）と呼びます。

2階ODEは物理学と工学の最も基本的な方程式であり、驚くほど多くの現象を記述します。

• 力学: バネ-質量-ダンパ系の振動（単振動、減衰振動、強制振動）
• 電気工学: RLC回路の過渡応答
• 量子力学: シュレーディンガー方程式（時間非依存）
• 構造力学: 梁のたわみ
• 制御工学: 2次系の過渡応答

特筆すべきは、力学のバネ-質量系と電気回路のRLC回路が全く同じ数学的構造を持つことです。質量はインダクタンスに、摩擦係数は抵抗に、バネ定数はキャパシタンスの逆数に対応します。2階ODEの解法を一度身につければ、分野を超えて応用できるのが大きな魅力です。

本記事の内容

• 2階線形ODEの一般論（重ね合わせの原理）
• 定数係数の場合の特性方程式による解法
• 3つのケース（過減衰・臨界減衰・振動減衰）
• 非斉次方程式の特殊解（未定係数法・定数変化法）
• Pythonによる振動系の解析と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 1階常微分方程式の解法まとめ — 1階ODEの基本
• 固有値と固有ベクトル — 特性方程式との関連

1階ODEの変数分離法や積分因子法を土台として、2階ODEでは「2つの独立な解の重ね合わせ」という新しい概念が登場します。また、特性方程式の解が固有値と密接に関連していることも見ていきます。

2階線形ODEの一般論

重ね合わせの原理

1階ODEの一般解は1つの任意定数を含んでいました。2階ODEでは、初期条件として位置 $y(x_0)$ と速度 $y'(x_0)$ の2つを指定する必要があるため、一般解は2つの任意定数を含みます。

これは解空間の次元と関係しています。2階線形斉次方程式 $y” + p(x)y’ + q(x)y = 0$ の解空間は2次元のベクトル空間です。線形代数の言葉で言えば、この空間の「基底」に相当する2つの線形独立な解を見つければ、すべての解はそれらの線形結合で表されます。

2つの線形独立な解 $y_1(x), y_2(x)$ が見つかれば、一般解は

$$ y = c_1 y_1(x) + c_2 y_2(x) $$

です。2つの定数 $c_1, c_2$ は初期条件 $y(x_0) = y_0, y'(x_0) = y_0’$ で決まります。初期条件を代入すると、$c_1, c_2$ に関する連立1次方程式

$$ \begin{cases} c_1 y_1(x_0) + c_2 y_2(x_0) = y_0 \\ c_1 y_1′(x_0) + c_2 y_2′(x_0) = y_0′ \end{cases} $$

が得られ、$y_1, y_2$ が線形独立であれば一意に解けます。

ロンスキアン

では、$y_1, y_2$ が「線形独立」であることをどう判定すればよいでしょうか。そのための道具がロンスキアン（Wronskian）です。

$$ W(y_1, y_2) = \begin{vmatrix} y_1 & y_2 \\ y_1′ & y_2′ \end{vmatrix} = y_1 y_2′ – y_2 y_1′ \neq 0 $$

ロンスキアンが恒等的にゼロでなければ、$y_1$ と $y_2$ は線形独立です。直感的には、ロンスキアンは「2つの解がどれだけ異なる振る舞いをしているか」を測る量です。例えば、$y_1 = e^x$ と $y_2 = e^{2x}$ のロンスキアンは $W = e^x \cdot 2e^{2x} – e^{2x} \cdot e^x = e^{3x} \neq 0$ なので線形独立です。一方、$y_1 = e^x$ と $y_2 = 3e^x$ なら $W = e^x \cdot 3e^x – 3e^x \cdot e^x = 0$ で線形従属です（片方が他方の定数倍）。

アーベルの公式

ロンスキアンの時間発展に関する有用な結果として、アーベルの公式があります。$y” + p(x)y’ + q(x)y = 0$ の任意の2つの解 $y_1, y_2$ に対して

$$ W(x) = W(x_0) \exp\left(-\int_{x_0}^x p(t) \, dt\right) $$

が成り立ちます。この公式から、ロンスキアンがある点でゼロでなければ全ての点でゼロでないこと（また、ある点でゼロなら全ての点でゼロであること）がわかります。

非斉次方程式の構造

非斉次方程式 $y” + py’ + qy = g(x)$ の一般解は

$$ y = y_h + y_p = (c_1 y_1 + c_2 y_2) + y_p $$

つまり、斉次の一般解 + 非斉次の特殊解という構造を持ちます。これは「全体的な解 = 自由な振る舞い + 外力に駆動された振る舞い」という物理的な解釈に対応しています。斉次解 $y_h$ は系の「固有の」振る舞い（自由振動）を表し、特殊解 $y_p$ は外力 $g(x)$ に対する「強制的な」応答を表します。

一般論の枠組みができたところで、具体的な解法に進みましょう。最も重要なのは、係数 $p, q$ が定数の場合です。この場合、特性方程式の代数的な解法に帰着できます。

定数係数の場合

特性方程式

定数係数の2階線形ODE

$$ \begin{equation} ay” + by’ + cy = 0 \end{equation} $$

を考えます。1階ODEの $y’ + \alpha y = 0$ の解が $y = e^{-\alpha x}$ であったことを思い出しましょう。2階の場合も指数関数 $y = e^{rx}$ を試みるのは自然な発想です。

$y = e^{rx}$ を代入すると $y’ = re^{rx}$, $y” = r^2 e^{rx}$ なので

$$ ar^2 e^{rx} + bre^{rx} + ce^{rx} = (ar^2 + br + c)e^{rx} = 0 $$

$e^{rx} \neq 0$ であるため

$$ ar^2 + br + c = 0 $$

これが特性方程式（characteristic equation）です。2階のODEが2次の代数方程式に帰着されました。2次方程式の解の公式から

$$ r = \frac{-b \pm \sqrt{b^2 – 4ac}}{2a} $$

判別式 $D = b^2 – 4ac$ の符号によって3つのケースに分かれ、それぞれが質的に異なる物理的振る舞いに対応します。

ケース1: $D > 0$（2つの異なる実数解）

特性方程式が2つの異なる実数解 $r_1 \neq r_2$ を持つ場合

$$ r_1, r_2 = \frac{-b \pm \sqrt{D}}{2a}, \quad y = c_1 e^{r_1 x} + c_2 e^{r_2 x} $$

$r_1, r_2$ がともに負であれば、解は単調にゼロに減衰します。物理的には過減衰（overdamping）に対応し、系は振動することなくゆっくりと平衡状態に戻ります。ドアクローザー（自動的にドアを閉めるダンパ機構）が過減衰に設計されている好例です — ドアがバタンと振動するのではなく、静かに閉まるように減衰を強くしています。

ケース2: $D = 0$（重根）

判別式がちょうどゼロの場合、特性方程式は重根 $r = r_1 = r_2$ を持ちます。

$$ r = -\frac{b}{2a} $$

しかし、$y_1 = e^{rx}$ は1つの解しか与えません。2階ODEの一般解には2つの独立な解が必要です。そこで、$y_2 = xe^{rx}$ が2番目の独立な解になることを確認しましょう。$y_2 = xe^{rx}$ を方程式に代入すると

$$ y_2′ = e^{rx} + rxe^{rx} = (1 + rx)e^{rx} $$

$$ y_2” = re^{rx} + r(1 + rx)e^{rx} = (2r + r^2 x)e^{rx} $$

これを $ay” + by’ + cy = 0$ に代入すると

$$ a(2r + r^2 x)e^{rx} + b(1 + rx)e^{rx} + cxe^{rx} = 0 $$

$e^{rx}$ でくくると

$$ (ar^2 + br + c)x + (2ar + b) = 0 $$

$r$ は特性方程式の根なので $ar^2 + br + c = 0$ です。また、重根の条件 $r = -b/(2a)$ から $2ar + b = 0$ です。したがって両方の項がゼロになり、$y_2 = xe^{rx}$ は確かに解です。

一般解は

$$ y = (c_1 + c_2 x)e^{rx} $$

です。物理的には臨界減衰（critical damping）に対応します。臨界減衰は「振動せずに最も速く平衡に戻る」状態であり、地震計や自動車のサスペンションなど、素早い応答が求められるシステムの設計目標となります。

ケース3: $D < 0$（複素数解）

判別式が負の場合、特性方程式の解は複素数になります。

$$ r = \alpha \pm i\beta, \quad \alpha = -\frac{b}{2a}, \quad \beta = \frac{\sqrt{|D|}}{2a} $$

$e^{(\alpha + i\beta)x} = e^{\alpha x}(\cos\beta x + i\sin\beta x)$ なので、実部と虚部を取ると2つの実数解が得られます。

$$ y = e^{\alpha x}(c_1 \cos\beta x + c_2 \sin\beta x) $$

この解は $e^{\alpha x}$ による振幅の変調（エンベロープ）と、$\cos\beta x$, $\sin\beta x$ による振動を組み合わせた形です。

• $\alpha < 0$ なら減衰振動（underdamping）: 振動しながら振幅が指数的に減衰します。ギターの弦をはじいた後の振動が典型例です
• $\alpha = 0$ なら単振動（harmonic oscillation）: 振幅が一定の永続的な振動です。理想的なバネ-質量系に対応します
• $\alpha > 0$ なら振幅が増大する不安定な振動です。実際の系では非線形効果が発生して振幅が制限されます

角周波数は $\beta$（振動の速さ）、減衰率は $|\alpha|$（振幅が減る速さ）です。振動の周期は $T = 2\pi/\beta$ であり、減衰のない場合の周期 $T_0 = 2\pi/\omega_0$ より長くなります。

ここまでで3つのケースの解を求めました。次に、これらの結果を力学系の具体例に適用してみましょう。

バネ-質量-ダンパ系

2階ODEの最も重要な物理的応用が、バネ-質量-ダンパ系です。質量 $m$ の物体が、バネ定数 $k$ のバネと減衰係数 $c$ のダンパに接続されている場合、ニュートンの運動方程式は

$$ m\ddot{x} + c\dot{x} + kx = 0 $$

です。ここで $m\ddot{x}$ は慣性力、$c\dot{x}$ は粘性抵抗力（速度に比例）、$kx$ は復元力（変位に比例）を表します。

この方程式を無次元化して見通しをよくしましょう。$\omega_0 = \sqrt{k/m}$（固有角周波数 — バネと質量だけで決まる振動の速さ）と $\zeta = c/(2\sqrt{mk})$（減衰比 — 減衰の強さを表す無次元量）を定義すると

$$ \ddot{x} + 2\zeta\omega_0\dot{x} + \omega_0^2 x = 0 $$

と書けます。系の振る舞いは $\zeta$ のたった1つのパラメータで完全に特徴づけられるのです。特性方程式は

$$ r^2 + 2\zeta\omega_0 r + \omega_0^2 = 0 $$

であり、解は $r = \omega_0(-\zeta \pm \sqrt{\zeta^2 – 1})$ です。判別式の符号は $\zeta^2 – 1$ で決まります。

減衰比 $\zeta$ は工学設計において極めて重要なパラメータです。例えば自動車のサスペンションでは $\zeta \approx 0.2$〜$0.4$ に設定されます。小さすぎると路面の凹凸で車体がいつまでも揺れ続け、大きすぎると路面の衝撃が直接乗員に伝わります。

RLC回路との対応

力学系のバネ-質量-ダンパ系と電気回路のRLC直列回路は、数学的に同じ構造を持ちます。RLC回路の回路方程式は

$$ L\ddot{q} + R\dot{q} + \frac{q}{C} = 0 $$

ここで $q$ は電荷、$L$ はインダクタンス、$R$ は抵抗、$C$ はキャパシタンスです。対応関係をまとめると

この対応により、力学系で得た解をそのまま電気回路に転用できます。電気回路の固有周波数は $\omega_0 = 1/\sqrt{LC}$、減衰比は $\zeta = R/(2\sqrt{L/C})$ です。

斉次方程式の自由振動を理解したところで、次に外力が加わる場合（非斉次方程式）の解法を見ていきましょう。強制振動や共振現象の理解に直結します。

非斉次方程式の特殊解

非斉次方程式 $ay” + by’ + cy = g(x)$ は、物理的には外力 $g(x)$ が加わった強制振動に対応します。一般解は $y = y_h + y_p$（斉次の一般解 + 特殊解）なので、特殊解 $y_p$ を見つけることが課題です。

未定係数法

$g(x)$ が多項式、指数関数、三角関数またはそれらの積の場合、特殊解の形を「推測」して係数を決定します。この方法は計算が簡単で、多くの実用的な問題に適用できます。

推測の指針は以下の通りです。

推測した形が斉次解と重複する場合は $x$ を掛けます。これは「共鳴」（resonance）に対応しており、外力の周波数が系の固有周波数に一致する場合に起こります。重根の場合は $x^2$ を掛ける必要があることもあります。

未定係数法の例

$y” + 4y = \cos 2x$ を解きましょう。斉次解は $y_h = c_1\cos 2x + c_2\sin 2x$ です。$g(x) = \cos 2x$ に対して $y_p = A\cos 2x + B\sin 2x$ を推測したいところですが、これは斉次解と重複します。そこで $x$ を掛けて $y_p = x(A\cos 2x + B\sin 2x)$ と推測します。

$y_p$ を方程式に代入して係数を比較すると（計算過程は省略）、$A = 0$, $B = 1/4$ が得られるので

$$ y_p = \frac{x}{4}\sin 2x $$

$x$ が掛かっているため、振幅が時間とともに線形に増大します。これは外力の周波数が固有周波数と一致する共振（resonance）現象です。

定数変化法

未定係数法が使えない一般の $g(x)$ に対しては、定数変化法（variation of parameters）を使います。このアプローチは、斉次解の定数 $c_1, c_2$ を関数 $u_1(x), u_2(x)$ に置き換えるものです。

特殊解を $y_p = u_1(x)y_1 + u_2(x)y_2$ とおきます。2つの未知関数 $u_1, u_2$ に対して条件が1つ（元の方程式）しかないので、追加条件として $u_1’y_1 + u_2’y_2 = 0$ を課します。この条件により $y_p$ の1階微分が簡単になります。

2つの条件から $u_1′, u_2’$ を求めると、クラメルの公式より

$$ u_1′ = \frac{-y_2 g}{aW}, \quad u_2′ = \frac{y_1 g}{aW} $$

ここで $W = y_1 y_2′ – y_2 y_1’$ はロンスキアンです。$u_1′, u_2’$ を積分すれば $u_1, u_2$ が得られ、特殊解が決まります。

定数変化法の利点は、$g(x)$ がどんな関数でも（積分が実行できる限り）適用できることです。一方、計算量は未定係数法より多くなりがちです。

強制振動と共振

非斉次方程式の重要な応用が強制振動です。$\ddot{x} + 2\zeta\omega_0\dot{x} + \omega_0^2 x = F_0\cos\omega t$ の定常解は

$$ x_p(t) = A(\omega)\cos(\omega t – \phi) $$

ここで振幅と位相は

$$ A(\omega) = \frac{F_0/\omega_0^2}{\sqrt{(1 – (\omega/\omega_0)^2)^2 + (2\zeta\omega/\omega_0)^2}} $$

$$ \phi = \arctan\frac{2\zeta\omega/\omega_0}{1 – (\omega/\omega_0)^2} $$

です。振幅 $A(\omega)$ は $\omega \approx \omega_0$ 付近で極大を取ります。これが共振（resonance）です。減衰比 $\zeta$ が小さいほど共振ピークは鋭く高くなり、$\zeta \to 0$ では $A \to \infty$（無限大の振幅）となります。

共振は工学において極めて重要な現象です。1940年のタコマナローズ橋の崩壊は、風による強制振動が橋の固有振動数と共振して振幅が増大した結果です。一方、ラジオの同調回路やMRI装置の共鳴吸収は、共振を積極的に利用した例です。

ここまでの理論をPythonで実装し、各ケースの振る舞いを可視化して理解を深めましょう。

Pythonでの実装

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# (a) 3つの減衰ケース
ax = axes[0, 0]
omega0 = 2*np.pi  # 固有周波数
t = np.linspace(0, 5, 1000)

for zeta, label, color in [(0.0, "Undamped (ζ=0)", "blue"),
                            (0.2, "Underdamped (ζ=0.2)", "green"),
                            (1.0, "Critical (ζ=1)", "red"),
                            (2.0, "Overdamped (ζ=2)", "purple")]:
    sol = solve_ivp(lambda t, y: [y[1], -2*zeta*omega0*y[1] - omega0**2*y[0]],
                    [0, 5], [1.0, 0.0], t_eval=t, max_step=0.01)
    ax.plot(sol.t, sol.y[0], color=color, linewidth=2, label=label)

ax.set_xlabel("Time", fontsize=12)
ax.set_ylabel("Displacement", fontsize=12)
ax.set_title("Damped Oscillator: $\\ddot{x} + 2\\zeta\\omega_0\\dot{x} + \\omega_0^2 x = 0$", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) 強制振動と共振
ax = axes[0, 1]
zeta = 0.1
omega_d_ratios = np.linspace(0.01, 3.0, 500)
amplitude = 1.0 / np.sqrt((1 - omega_d_ratios**2)**2 + (2*zeta*omega_d_ratios)**2)

for z, color in [(0.05, "red"), (0.1, "blue"), (0.3, "green"), (0.7, "purple")]:
    amp = 1.0 / np.sqrt((1 - omega_d_ratios**2)**2 + (2*z*omega_d_ratios)**2)
    ax.plot(omega_d_ratios, amp, color=color, linewidth=2, label=f"ζ = {z}")

ax.set_xlabel("$\\omega/\\omega_0$", fontsize=12)
ax.set_ylabel("Amplitude ratio", fontsize=12)
ax.set_title("Frequency Response (Resonance)", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(0, 12)

# (c) RLC回路
ax = axes[1, 0]
R_vals = [10, 100, 200, 632]  # 632Ω ≈ 臨界減衰
L = 0.1  # H
C_val = 1e-6  # F
E = 5.0  # V

for R in R_vals:
    sol = solve_ivp(lambda t, y: [y[1], (E - R*y[1] - y[0]/C_val)/L],
                    [0, 0.005], [0, 0], t_eval=np.linspace(0, 0.005, 1000))
    ax.plot(sol.t*1000, sol.y[0]*1e6, linewidth=2, label=f"R = {R} Ω")

ax.set_xlabel("Time (ms)", fontsize=12)
ax.set_ylabel("Charge (μC)", fontsize=12)
ax.set_title("RLC Circuit Step Response", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (d) 位相平面
ax = axes[1, 1]
zeta_val = 0.3
for x0 in [1, 2, 3]:
    for v0 in [-5, 0, 5]:
        sol = solve_ivp(lambda t, y: [y[1], -2*zeta_val*omega0*y[1] - omega0**2*y[0]],
                        [0, 5], [x0, v0], t_eval=np.linspace(0, 5, 500))
        ax.plot(sol.y[0], sol.y[1], linewidth=0.8, alpha=0.7)

ax.plot(0, 0, "ro", markersize=8, label="Equilibrium")
ax.set_xlabel("$x$ (displacement)", fontsize=12)
ax.set_ylabel("$\\dot{x}$ (velocity)", fontsize=12)
ax.set_title(f"Phase Portrait (ζ = {zeta_val})", fontsize=12)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")

plt.tight_layout()
plt.savefig("ode_second_order.png", dpi=150, bbox_inches="tight")
plt.show()

このグラフから、2階ODEの豊かな振る舞いが読み取れます。

1. 減衰振動の4つのケース（左上）: 減衰比 $\zeta$ による4つの振る舞いが明確に区別されています。非減衰（$\zeta = 0$、青）は一定振幅で振動し続け、減衰振動（$\zeta = 0.2$、緑）は振動しながら振幅が指数的に減衰しています。臨界減衰（$\zeta = 1$、赤）は振動せずに最も速く平衡に到達しており、過減衰（$\zeta = 2$、紫）は臨界減衰よりもさらにゆっくりとゼロに近づいています。「最速で平衡に戻る」という臨界減衰の特性が、赤い曲線が紫い曲線よりも先にゼロ付近に到達していることから視覚的に確認できます

2. 周波数応答と共振（右上）: 周波数応答曲線で、$\omega/\omega_0 = 1$ 付近で共振ピークが発生しています。$\zeta = 0.05$（赤）では振幅が約10倍に増幅されており、共振の危険性が見て取れます。減衰が大きくなるにつれて共振ピークは低く広くなり、$\zeta = 0.7$ ではピークがほぼ消失しています。$\omega/\omega_0 \gg 1$ ではすべての曲線が急速にゼロに近づいており、高周波の外力に対して系が応答できないことを示しています。これはローパスフィルタとしての2次系の特性そのものです

3. RLC回路の過渡応答（左下）: 抵抗値 $R$ を変えることで、同一のRLC回路が振動的（低抵抗）から過減衰（高抵抗）まで変化する様子が見えます。$R = 632$ $\Omega$ 付近が臨界減衰に対応しており、振動せずに目標値に最速で到達しています。低抵抗（$R = 10$ $\Omega$）の場合は激しく振動しながら定常値に近づいており、これは力学系の減衰振動と完全に対応しています

4. 位相平面図（右下）: 位相平面（$x$-$\dot{x}$ 平面）では、異なる初期条件からの軌道がすべて原点（平衡点）に向かってスパイラル状に収束しており、安定な減衰振動を示しています。スパイラルが反時計回りに巻いているのは、変位が正のとき速度が負に転じて変位が減少するという振動のサイクルを反映しています。初期条件に依存せず全ての軌道が同じ平衡点に収束するのは、系が大域的に安定であることを意味しています

具体的な数値計算例

理論の理解を確実にするために、具体的な数値を用いた計算例を見ていきましょう。

例題1: 特性方程式の解法

$y” + 5y’ + 6y = 0$ を解きます。

特性方程式は $r^2 + 5r + 6 = 0$ です。因数分解すると $(r + 2)(r + 3) = 0$ なので $r_1 = -2$, $r_2 = -3$ です。判別式 $D = 25 – 24 = 1 > 0$ なので、2つの異なる実数解（ケース1）です。

一般解は

$$ y = c_1 e^{-2x} + c_2 e^{-3x} $$

初期条件 $y(0) = 1$, $y'(0) = 0$ を適用しましょう。$y(0) = c_1 + c_2 = 1$ と $y'(0) = -2c_1 – 3c_2 = 0$ から

$$ c_1 + c_2 = 1, \quad -2c_1 – 3c_2 = 0 $$

2番目の式から $c_1 = -3c_2/2$ を1番目に代入すると $-3c_2/2 + c_2 = 1$ より $c_2 = -2$、$c_1 = 3$ です。

$$ y = 3e^{-2x} – 2e^{-3x} $$

$x \to \infty$ で $y \to 0$ であり、振動することなく（両方の根が負の実数なので）平衡に戻る過減衰の解です。

例題2: 減衰振動の解

$y” + 2y’ + 5y = 0$, $y(0) = 2$, $y'(0) = 0$ を解きます。

特性方程式 $r^2 + 2r + 5 = 0$ の解は

$$ r = \frac{-2 \pm \sqrt{4 – 20}}{2} = \frac{-2 \pm \sqrt{-16}}{2} = -1 \pm 2i $$

$D = -16 < 0$ なのでケース3（複素数解）です。$\alpha = -1$, $\beta = 2$ として

$$ y = e^{-x}(c_1\cos 2x + c_2\sin 2x) $$

初期条件を適用します。$y(0) = c_1 = 2$ です。$y’ = -e^{-x}(c_1\cos 2x + c_2\sin 2x) + e^{-x}(-2c_1\sin 2x + 2c_2\cos 2x)$ なので

$$ y'(0) = -c_1 + 2c_2 = -2 + 2c_2 = 0 \implies c_2 = 1 $$

$$ y = e^{-x}(2\cos 2x + \sin 2x) $$

この解は角周波数 $\omega = 2$ で振動しながら、$e^{-x}$ のエンベロープで減衰します。振動の周期は $T = 2\pi/2 = \pi \approx 3.14$ です。

例題3: バネ-質量系の設計

質量 $m = 2$ kg の物体をバネ（$k = 200$ N/m）とダンパに取り付けます。

固有周波数: $\omega_0 = \sqrt{k/m} = \sqrt{200/2} = 10$ rad/s、周波数 $f_0 = \omega_0/(2\pi) \approx 1.59$ Hz

臨界減衰のダンパ係数: $\zeta = 1$ とすると

$$ c = 2\zeta\sqrt{mk} = 2 \times 1 \times \sqrt{2 \times 200} = 2\sqrt{400} = 40 \text{ N$\cdot$s/m} $$

$c < 40$ N$\cdot$s/m なら振動的、$c > 40$ N$\cdot$s/m なら過減衰です。

例えば $c = 8$ N$\cdot$s/m（$\zeta = 0.2$）の場合、減衰振動の角周波数は

$$ \omega_d = \omega_0\sqrt{1 – \zeta^2} = 10\sqrt{1 – 0.04} = 10 \times 0.98 = 9.80 \text{ rad/s} $$

減衰によって振動周波数がわずかに低下しますが、$\zeta = 0.2$ 程度ではほとんど変化しません。

例題4: RLC回路の臨界減衰条件

$L = 10$ mH、$C = 1$ $\mu$F のRLC直列回路で臨界減衰を実現する抵抗値を求めます。

$$ \omega_0 = \frac{1}{\sqrt{LC}} = \frac{1}{\sqrt{10 \times 10^{-3} \times 10^{-6}}} = \frac{1}{\sqrt{10^{-8}}} = \frac{1}{10^{-4}} = 10^4 \text{ rad/s} $$

周波数にすると $f_0 = 10^4/(2\pi) \approx 1592$ Hz です。

臨界減衰の条件 $\zeta = 1$ から

$$ R = 2\sqrt{\frac{L}{C}} = 2\sqrt{\frac{10 \times 10^{-3}}{10^{-6}}} = 2\sqrt{10^4} = 2 \times 100 = 200 \text{ $\Omega$} $$

$R < 200$ $\Omega$ では電荷が振動的に変化し、$R > 200$ $\Omega$ では過減衰になります。

例題5: 未定係数法による特殊解

$y” + 3y’ + 2y = e^{-x}$ の特殊解を求めます。

まず斉次解を求めます。特性方程式 $r^2 + 3r + 2 = 0$ を解くと $(r+1)(r+2) = 0$ より $r = -1, -2$ です。$y_h = c_1 e^{-x} + c_2 e^{-2x}$ です。

$g(x) = e^{-x}$ に対して $y_p = Ae^{-x}$ を推測しますが、$e^{-x}$ は斉次解 $y_1 = e^{-x}$ と重複しています。そこで $y_p = Axe^{-x}$ を推測します。

$y_p = Axe^{-x}$ を微分すると

$$ y_p’ = A(1 – x)e^{-x} $$

$$ y_p” = A(-1 – (1-x))e^{-x} = A(x – 2)e^{-x} $$

元の方程式に代入すると

$$ A(x-2)e^{-x} + 3A(1-x)e^{-x} + 2Axe^{-x} = e^{-x} $$

$e^{-x}$ でくくると

$$ A[(x-2) + 3(1-x) + 2x] = 1 $$

$$ A[x – 2 + 3 – 3x + 2x] = A \cdot 1 = 1 $$

よって $A = 1$ であり、特殊解は $y_p = xe^{-x}$ です。

一般解は

$$ y = c_1 e^{-x} + c_2 e^{-2x} + xe^{-x} $$

特殊解の項 $xe^{-x}$ は斉次解 $e^{-x}$ に $x$ が掛かった形をしており、共振に対応する振る舞い（振幅の線形増大）を示しています。

まとめ

本記事では、2階線形常微分方程式の解法について、一般論から具体的な応用まで解説しました。

• 解空間: 2階線形斉次ODEの解空間は2次元のベクトル空間であり、2つの線形独立な解の線形結合が一般解となります。ロンスキアンで独立性を判定できます
• 特性方程式: 定数係数の場合、$ar^2 + br + c = 0$ の判別式 $D$ で3つのケース（過減衰 $D>0$・臨界減衰 $D=0$・振動減衰 $D<0$）に分類されます
• 減衰比 $\zeta$ が系の振動特性を完全に決定し、$\zeta = 1$ が臨界減衰です。力学系と電気回路の間には数学的な完全な対応関係があります
• 非斉次方程式の特殊解は、$g(x)$ が簡単な関数なら未定係数法、一般の場合は定数変化法で求めます
• 共振現象: 強制振動の周波数が固有周波数に近いとき振幅が増大し、工学的に重要な注意を要する現象です

2階ODEは偏微分方程式を変数分離法で解く際にも頻出し、フーリエ級数やラプラス変換とも密接に関連しています。また、より高階のODEや連立ODEは行列の固有値問題に帰着されるため、線形代数との接点も多い分野です。

次のステップとして、以下の記事も参考にしてください。

• 級数解法とべき級数展開による微分方程式の解法 — 変数係数の場合の解法を学びます
• 熱伝導方程式の導出と解法 — 変数分離法で2階ODEが自然に現れます
• 常微分方程式の数値解法まとめ — 解析解が得られない場合の数値的アプローチです