バネで結ばれた2つの質点の振動を記述するとき、各質点の変位は互いに影響し合います。一方の質点の位置が他方に力を及ぼし、他方の応答がさらに一方に返ってくる — このような「相互作用する量」のダイナミクスを記述するのが連立常微分方程式（system of ODEs）です。
単独の1変数ODEの世界では、特性方程式や変数分離法が威力を発揮しました。しかし複数の変数が結合した系では、これらの手法をそのまま適用することはできません。ここで救いの手を差し伸べるのが線形代数です。連立1階線形ODEは行列とベクトルの言葉で書くことができ、その解は行列指数関数 $e^{At}$ によって統一的に表現されます。
連立ODEの理論を理解すると、以下のような場面で活用できます。
- 振動工学: 多自由度振動系の固有振動数とモード解析
- 制御工学: 状態空間モデル $\dot{\bm{x}} = A\bm{x} + B\bm{u}$ の解の構造
- 生態学: ロトカ=ヴォルテラ方程式（捕食者-被食者モデル）
- 電気回路: RLC回路網の過渡応答解析
本記事の内容
- 連立1階ODEとベクトル表記
- 定数係数連立線形ODEの解法（固有値・固有ベクトル法）
- 行列指数関数の定義と性質
- 固有値による解の分類と位相図
- 非同次系と定数変化法
- Pythonでの実装と位相平面の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 固有値分解の理論 — 固有値・固有ベクトルの計算
- 2階線形常微分方程式の解法 — 単独ODEの解法
- 行列の微積分 — 行列関数の微分
高階ODEから連立1階ODEへの変換
なぜ連立1階系を考えるのか
2階ODE $y” + 3y’ + 2y = 0$ は、$x_1 = y$, $x_2 = y’$ と置くことで連立1階系に変換できます。
$$ \begin{cases} x_1′ = x_2 \\ x_2′ = -2x_1 – 3x_2 \end{cases} $$
これは一般的な現象です。任意の $n$ 階ODEは、$n$ 個の1階ODEの連立系に書き直せます。つまり、連立1階ODEの理論がすべてのODEの基盤になるのです。
ベクトル・行列表記
上の連立系をベクトルと行列で書くと、非常にコンパクトになります。
$$ \bm{x}’ = A\bm{x} $$
ここで
$$ \bm{x} = \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}, \quad A = \begin{pmatrix} 0 & 1 \\ -2 & -3 \end{pmatrix} $$
この形式は、変数がいくつあっても $\bm{x}’ = A\bm{x}$ という同じ形で表せます。2変数でも100変数でも、数学的な構造は同じです。
一般の $n$ 階ODEの変換
$n$ 階ODE
$$ y^{(n)} + a_{n-1}y^{(n-1)} + \cdots + a_1 y’ + a_0 y = 0 $$
に対して $x_k = y^{(k-1)}$（$k = 1, 2, \ldots, n$）と置くと
$$ \bm{x}’ = \begin{pmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ -a_0 & -a_1 & -a_2 & \cdots & -a_{n-1} \end{pmatrix} \bm{x} $$
この行列はコンパニオン行列（companion matrix）と呼ばれます。コンパニオン行列の特性多項式は元のODEの特性方程式と一致するため、固有値法との接続が自然に成り立ちます。
連立系の表記が手に入ったところで、次はこの系の解を具体的に求める方法を見ていきましょう。
固有値・固有ベクトル法による解法
スカラーODEとの類比
スカラーODE $x’ = ax$ の解は $x(t) = Ce^{at}$ です。この類推から、連立系 $\bm{x}’ = A\bm{x}$ の解として $\bm{x}(t) = \bm{v}e^{\lambda t}$（$\bm{v}$ は定ベクトル）を試してみましょう。
代入すると
$$ \lambda \bm{v} e^{\lambda t} = A \bm{v} e^{\lambda t} $$
$e^{\lambda t} \neq 0$ で割ると
$$ A\bm{v} = \lambda \bm{v} $$
これは固有値問題そのものです。つまり、$\lambda$ が $A$ の固有値、$\bm{v}$ が対応する固有ベクトルであれば、$\bm{x}(t) = \bm{v}e^{\lambda t}$ は解になります。
2次元系の一般解（実固有値の場合）
$A$ が2つの相異なる実固有値 $\lambda_1, \lambda_2$ を持ち、対応する固有ベクトルが $\bm{v}_1, \bm{v}_2$ であれば、一般解は
$$ \bm{x}(t) = c_1 \bm{v}_1 e^{\lambda_1 t} + c_2 \bm{v}_2 e^{\lambda_2 t} $$
定数 $c_1, c_2$ は初期条件 $\bm{x}(0) = \bm{x}_0$ から決定されます。
$$ \bm{x}_0 = c_1 \bm{v}_1 + c_2 \bm{v}_2 $$
これは $\bm{x}_0$ を固有ベクトル基底で展開する問題に帰着します。
具体例: 2つの実固有値
$$ A = \begin{pmatrix} -1 & 2 \\ 0 & -3 \end{pmatrix} $$
特性方程式 $\det(A – \lambda I) = 0$ を計算すると
$$ (-1 – \lambda)(-3 – \lambda) = 0 $$
よって $\lambda_1 = -1$, $\lambda_2 = -3$ です。
$\lambda_1 = -1$ に対する固有ベクトルは $(A + I)\bm{v} = \bm{0}$ を解いて
$$ \begin{pmatrix} 0 & 2 \\ 0 & -2 \end{pmatrix} \bm{v} = \bm{0} \implies \bm{v}_1 = \begin{pmatrix} 1 \\ 0 \end{pmatrix} $$
$\lambda_2 = -3$ に対しては $(A + 3I)\bm{v} = \bm{0}$ を解いて
$$ \begin{pmatrix} 2 & 2 \\ 0 & 0 \end{pmatrix} \bm{v} = \bm{0} \implies \bm{v}_2 = \begin{pmatrix} 1 \\ -1 \end{pmatrix} $$
一般解は
$$ \bm{x}(t) = c_1 \begin{pmatrix} 1 \\ 0 \end{pmatrix} e^{-t} + c_2 \begin{pmatrix} 1 \\ -1 \end{pmatrix} e^{-3t} $$
両方の固有値が負なので、解は $t \to \infty$ で原点に収束します。$e^{-t}$ の方が減衰が遅いため、十分な時間が経つと解は $\bm{v}_1$ の方向に沿って原点に向かいます。
複素固有値の場合
固有値が複素数 $\lambda = \alpha \pm \beta i$（$\beta \neq 0$）の場合、対応する固有ベクトルも複素数になります。$\bm{v} = \bm{p} + i\bm{q}$ として、オイラーの公式を使うと実数解を構成できます。
$\bm{v}e^{(\alpha + \beta i)t}$ の実部と虚部がそれぞれ独立な実数解になります。
$$ \bm{x}_1(t) = e^{\alpha t}(\bm{p}\cos\beta t – \bm{q}\sin\beta t) $$
$$ \bm{x}_2(t) = e^{\alpha t}(\bm{p}\sin\beta t + \bm{q}\cos\beta t) $$
一般解は $\bm{x}(t) = c_1 \bm{x}_1(t) + c_2 \bm{x}_2(t)$ です。$\alpha < 0$ なら減衰振動、$\alpha > 0$ なら発散振動、$\alpha = 0$ なら中心（周期軌道）になります。
重複固有値の場合
$A$ が重複固有値 $\lambda$ を持ち、独立な固有ベクトルが1つしかない場合（固有空間の次元が固有値の重複度より小さい場合）、一般化固有ベクトルが必要になります。
$\bm{v}_1$ を固有ベクトル、$\bm{v}_2$ を $(A – \lambda I)\bm{v}_2 = \bm{v}_1$ を満たす一般化固有ベクトルとすると
$$ \bm{x}(t) = c_1 \bm{v}_1 e^{\lambda t} + c_2 (\bm{v}_1 t + \bm{v}_2) e^{\lambda t} $$
スカラーの場合の $te^{\lambda t}$ に対応する項が現れていることに注目してください。
固有値法は直感的でわかりやすいですが、すべての場合を統一的に扱うにはより強力な道具が必要です。それが行列指数関数です。
行列指数関数
定義
スカラーの指数関数 $e^{at} = \sum_{k=0}^{\infty} \frac{(at)^k}{k!}$ をそのまま行列に拡張します。
$$ \begin{equation} e^{At} = \sum_{k=0}^{\infty} \frac{(At)^k}{k!} = I + At + \frac{A^2 t^2}{2!} + \frac{A^3 t^3}{3!} + \cdots \end{equation} $$
この級数は任意の正方行列 $A$ に対して絶対収束します（ノルムの意味で）。
基本性質
行列指数関数は以下の性質を持ちます。
$$ \frac{d}{dt}e^{At} = Ae^{At} = e^{At}A $$
この性質から、$\bm{x}’ = A\bm{x}$ の解が直ちに得られます。
$$ \begin{equation} \bm{x}(t) = e^{At}\bm{x}_0 \end{equation} $$
初期条件 $\bm{x}(0) = \bm{x}_0$ を確認すると $e^{A \cdot 0}\bm{x}_0 = I\bm{x}_0 = \bm{x}_0$ で確かに成り立ちます。
また、$e^{A(t+s)} = e^{At}e^{As}$ が成り立つため、行列指数関数は半群を形成します（ただし一般に $e^{(A+B)t} \neq e^{At}e^{Bt}$ であることに注意。等号が成り立つのは $AB = BA$ のときのみ）。
対角化による計算
$A$ が対角化可能、すなわち $A = PDP^{-1}$（$D$ は固有値の対角行列）と書けるとき
$$ A^k = PD^kP^{-1} $$
これを級数の定義に代入すると
$$ e^{At} = P\left(\sum_{k=0}^{\infty} \frac{D^k t^k}{k!}\right)P^{-1} = Pe^{Dt}P^{-1} $$
対角行列の指数関数は簡単に計算できます。
$$ e^{Dt} = \begin{pmatrix} e^{\lambda_1 t} & 0 & \cdots \\ 0 & e^{\lambda_2 t} & \cdots \\ \vdots & & \ddots \end{pmatrix} $$
つまり、$e^{At}$ の計算は固有値分解に帰着するのです。
ジョルダン標準形による計算
$A$ が対角化不可能な場合、ジョルダン標準形 $A = PJP^{-1}$ を使います。ジョルダンブロック
$$ J_k(\lambda) = \begin{pmatrix} \lambda & 1 & 0 & \cdots \\ 0 & \lambda & 1 & \cdots \\ \vdots & & \ddots & \ddots \\ 0 & \cdots & 0 & \lambda \end{pmatrix} $$
に対する行列指数関数は
$$ e^{J_k(\lambda)t} = e^{\lambda t} \begin{pmatrix} 1 & t & \frac{t^2}{2!} & \cdots \\ 0 & 1 & t & \cdots \\ \vdots & & \ddots & \ddots \\ 0 & \cdots & 0 & 1 \end{pmatrix} $$
$t$ の多項式が現れるのは、重複固有値の場合に $te^{\lambda t}$ の項が解に現れることの行列版です。
行列指数関数の理論を手にしたことで、解の構造を統一的に理解できるようになりました。次は、この解がどのような幾何学的振る舞いを示すかを分類しましょう。
2次元系の位相図による分類
固有値と解の挙動の対応
2次元定数係数系 $\bm{x}’ = A\bm{x}$ の平衡点（原点）の分類は、$A$ の固有値 $\lambda_1, \lambda_2$ のみで決まります。以下の表に主要な場合をまとめます。
| 固有値の種類 | 条件 | 平衡点の型 | 安定性 |
|---|---|---|---|
| 実数・異符号 | $\lambda_1 > 0 > \lambda_2$ | 鞍点（saddle） | 不安定 |
| 実数・共に負 | $\lambda_1 < \lambda_2 < 0$ | 安定結節点（stable node） | 漸近安定 |
| 実数・共に正 | $0 < \lambda_1 < \lambda_2$ | 不安定結節点（unstable node） | 不安定 |
| 複素数・負実部 | $\alpha \pm \beta i$, $\alpha < 0$ | 安定渦巻点（stable spiral） | 漸近安定 |
| 複素数・正実部 | $\alpha \pm \beta i$, $\alpha > 0$ | 不安定渦巻点（unstable spiral） | 不安定 |
| 純虚数 | $\pm \beta i$ | 中心（center） | 安定（漸近安定でない） |
| 重複・退化 | $\lambda_1 = \lambda_2$, 1固有ベクトル | 退化結節点（degenerate node） | $\lambda < 0$ で漸近安定 |
鞍点（saddle point）
固有値が異符号 $\lambda_1 > 0 > \lambda_2$ のとき、$\bm{v}_1$ 方向には指数的に発散し、$\bm{v}_2$ 方向には指数的に収束します。軌道は双曲線的な形状を描きます。
鞍点は不安定ですが、$\bm{v}_2$ 方向の安定多様体（stable manifold）と $\bm{v}_1$ 方向の不安定多様体（unstable manifold）は重要な概念です。これらは非線形系の安定性解析でも中心的な役割を果たします。
安定結節点（stable node）
固有値が共に負で実数の場合、すべての軌道は原点に収束します。$|\lambda_1| < |\lambda_2|$ のとき、$e^{\lambda_2 t}$ の方が速く減衰するため、十分な時間が経つと軌道は遅い固有方向 $\bm{v}_1$ に沿って原点に接近します。
渦巻点（spiral point）
複素固有値 $\alpha \pm \beta i$ の場合、解は振動しながら原点に近づく（$\alpha < 0$）か、離れていく（$\alpha > 0$）かします。振動の周期は $2\pi/\beta$ です。
トレースと行列式による分類
固有値を直接計算しなくても、2次元の場合は $A$ のトレース $\tau = \text{tr}(A) = \lambda_1 + \lambda_2$ と行列式 $\Delta = \det(A) = \lambda_1 \lambda_2$ から分類できます。
- $\Delta < 0$: 鞍点
- $\Delta > 0$ かつ $\tau^2 – 4\Delta > 0$: 結節点（$\tau < 0$ で安定、$\tau > 0$ で不安定）
- $\Delta > 0$ かつ $\tau^2 – 4\Delta < 0$: 渦巻点（$\tau < 0$ で安定、$\tau > 0$ で不安定）
- $\Delta > 0$ かつ $\tau = 0$: 中心
判別式 $\tau^2 – 4\Delta$ は特性方程式 $\lambda^2 – \tau\lambda + \Delta = 0$ の判別式そのものです。
位相図の理論を理解したところで、次はPythonでこれらを実際に可視化してみましょう。
非同次系と定数変化法
非同次系の一般解
非同次系
$$ \bm{x}’ = A\bm{x} + \bm{f}(t) $$
の解は、同次解（余関数）と特殊解の和で構成されます。
$$ \bm{x}(t) = e^{At}\bm{x}_0 + \int_0^t e^{A(t-s)}\bm{f}(s)\,ds $$
この公式は定数変化法（variation of parameters）から導かれます。
導出
$\bm{x}(t) = e^{At}\bm{c}(t)$ と置きます（定数 $\bm{c}$ を $t$ の関数に「変化」させる）。
$\bm{x}’ = A\bm{x} + \bm{f}(t)$ に代入すると
$$ Ae^{At}\bm{c}(t) + e^{At}\bm{c}'(t) = Ae^{At}\bm{c}(t) + \bm{f}(t) $$
第1項が相殺されて
$$ e^{At}\bm{c}'(t) = \bm{f}(t) $$
両辺に $e^{-At}$ を左から掛けて積分すると
$$ \bm{c}(t) = \bm{c}(0) + \int_0^t e^{-As}\bm{f}(s)\,ds $$
$\bm{x}(t) = e^{At}\bm{c}(t)$ に戻すと、先ほどの公式が得られます。
物理的解釈
$$ \bm{x}(t) = \underbrace{e^{At}\bm{x}_0}_{\text{自由応答}} + \underbrace{\int_0^t e^{A(t-s)}\bm{f}(s)\,ds}_{\text{強制応答}} $$
第1項は初期条件のみによる応答（外力がない場合の振る舞い）、第2項は外力 $\bm{f}(s)$ に対する応答です。$e^{A(t-s)}$ は時刻 $s$ に加えられたインパルスが時刻 $t$ でどれだけの効果を持つかを表すグリーン関数（インパルス応答）に対応します。
この積分は制御工学で畳み込み積分として頻繁に登場し、系の入出力関係を記述する基盤となります。
では、ここまでの理論をPythonで実装し、位相図を描いてみましょう。
Pythonによる実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm
# 4種類の平衡点を可視化
fig, axes = plt.subplots(2, 2, figsize=(12, 12))
# --- (a) 鞍点 (saddle): 固有値 2, -1 ---
A_saddle = np.array([[0.5, 1.5], [1.5, 0.5]])
# --- (b) 安定結節点 (stable node): 固有値 -1, -3 ---
A_node = np.array([[-1, 2], [0, -3]])
# --- (c) 安定渦巻点 (stable spiral): 固有値 -0.2 ± 2i ---
A_spiral = np.array([[-0.2, 2], [-2, -0.2]])
# --- (d) 中心 (center): 固有値 ±2i ---
A_center = np.array([[0, 2], [-2, 0]])
matrices = [A_saddle, A_node, A_spiral, A_center]
titles = [
f"(a) Saddle Point\n$\\lambda = 2, -1$",
f"(b) Stable Node\n$\\lambda = -1, -3$",
f"(c) Stable Spiral\n$\\lambda = -0.2 \\pm 2i$",
f"(d) Center\n$\\lambda = \\pm 2i$"
]
for idx, (A, title) in enumerate(zip(matrices, titles)):
ax = axes[idx // 2][idx % 2]
# 初期条件を円周上に配置
n_traj = 16
thetas = np.linspace(0, 2 * np.pi, n_traj, endpoint=False)
r0 = 2.0
t_span = np.linspace(0, 5, 500)
for theta in thetas:
x0 = r0 * np.array([np.cos(theta), np.sin(theta)])
# 行列指数関数で解を計算
traj = np.array([expm(A * t) @ x0 for t in t_span])
# 表示範囲内のみプロット
mask = (np.abs(traj[:, 0]) < 4) & (np.abs(traj[:, 1]) < 4)
ax.plot(traj[mask, 0], traj[mask, 1], "b-", linewidth=0.7, alpha=0.6)
# 固有値・固有ベクトルの表示（実固有値の場合）
eigvals, eigvecs = np.linalg.eig(A)
if np.all(np.isreal(eigvals)):
for i in range(2):
v = eigvecs[:, i].real
ax.annotate("", xy=v * 3, xytext=-v * 3,
arrowprops=dict(arrowstyle="-", color="red",
linewidth=2, linestyle="--"))
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title(title, fontsize=13)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.plot(0, 0, "ko", markersize=5)
plt.tight_layout()
plt.savefig("ode_system_phase_portraits.png", dpi=150, bbox_inches="tight")
plt.show()
この4つの位相図から、固有値と解の幾何学的性質の関係が明確に読み取れます。
-
左上（鞍点）: 2つの固有ベクトル方向（赤い破線）に沿って、一方では発散し他方では収束しています。軌道は双曲線的な形状を描いており、安定多様体（$\lambda = -1$ の方向）に沿ってのみ原点に到達できることがわかります
-
右上（安定結節点）: すべての軌道が原点に収束しています。$\lambda_2 = -3$ 方向の成分が速く減衰するため、軌道は最終的に遅い方向（$\lambda_1 = -1$ の固有ベクトル）に沿って原点に接近しています
-
左下（安定渦巻点）: 軌道が反時計回りに渦を巻きながら原点に収束しています。虚部 $\beta = 2$ が振動の角周波数を、実部 $\alpha = -0.2$ が減衰率を決めています
-
右下（中心）: 軌道が閉じた楕円を描いており、原点には収束も発散もしません。これは $\alpha = 0$（純虚数固有値）に対応し、エネルギーが保存される保存系の典型的な振る舞いです
次に、行列指数関数の計算と非同次系の解を数値的に確認してみましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm
from scipy.integrate import solve_ivp
# 非同次系: x' = Ax + f(t) の解を比較
# A = [[-1, 2], [0, -3]], f(t) = [sin(t), 0]
A = np.array([[-1, 2], [0, -3]])
x0 = np.array([1.0, 0.5])
def system(t, x):
f = np.array([np.sin(t), 0])
return A @ x + f
# 数値解（Runge-Kutta法）
t_span = (0, 10)
t_eval = np.linspace(0, 10, 500)
sol = solve_ivp(system, t_span, x0, t_eval=t_eval, rtol=1e-10)
# 行列指数関数による解（数値積分で検証）
def exact_solution(t, A, x0, dt=0.001):
"""行列指数関数を使った厳密解の数値近似"""
# 同次部分
x_homo = expm(A * t) @ x0
# 非同次部分（数値積分）
n_steps = int(t / dt) + 1
s_vals = np.linspace(0, t, n_steps)
integral = np.zeros(2)
for i in range(1, len(s_vals)):
s = s_vals[i]
ds = s_vals[i] - s_vals[i-1]
f_s = np.array([np.sin(s), 0])
integral += expm(A * (t - s)) @ f_s * ds
return x_homo + integral
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (a) 時系列
ax = axes[0]
ax.plot(sol.t, sol.y[0], "b-", linewidth=2, label="$x_1(t)$ (RK45)")
ax.plot(sol.t, sol.y[1], "r-", linewidth=2, label="$x_2(t)$ (RK45)")
# 行列指数関数で数点を検証
t_check = np.array([1, 2, 3, 5, 7, 10])
x_exact = np.array([exact_solution(t, A, x0) for t in t_check])
ax.plot(t_check, x_exact[:, 0], "bo", markersize=8, label="$x_1$ (matrix exp)")
ax.plot(t_check, x_exact[:, 1], "rs", markersize=8, label="$x_2$ (matrix exp)")
ax.set_xlabel("$t$", fontsize=12)
ax.set_ylabel("$x(t)$", fontsize=12)
ax.set_title("Non-homogeneous System: Time Series", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
# (b) 位相平面
ax = axes[1]
ax.plot(sol.y[0], sol.y[1], "b-", linewidth=2, label="Trajectory")
ax.plot(x0[0], x0[1], "go", markersize=10, label="Initial condition")
# 同次系の軌道も比較
def system_homo(t, x):
return A @ x
sol_homo = solve_ivp(system_homo, t_span, x0, t_eval=t_eval, rtol=1e-10)
ax.plot(sol_homo.y[0], sol_homo.y[1], "r--", linewidth=1.5, alpha=0.7,
label="Homogeneous (no forcing)")
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title("Phase Plane", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
plt.tight_layout()
plt.savefig("ode_system_nonhomogeneous.png", dpi=150, bbox_inches="tight")
plt.show()
この結果から、非同次系の振る舞いが読み取れます。
-
左図（時系列）: 数値解（RK45, 実線）と行列指数関数による解（マーカー）が完全に一致しており、理論の正しさが数値的に検証されています。$x_2(t)$ は同次系の減衰のみ受けるため急速に0に収束し、$x_1(t)$ は外力 $\sin(t)$ の影響で定常振動に漸近しています
-
右図（位相平面）: 同次系（赤の破線）は原点に直接収束しますが、非同次系（青の実線）は原点の周りの楕円的な軌道に漸近しています。これは外力による定常応答が同次解の減衰後に支配的になるためです
トレース-行列式平面による分類
import numpy as np
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1, figsize=(9, 7))
# トレース-行列式平面
tau = np.linspace(-4, 4, 500)
# 判別式 = 0 の境界: delta = tau^2 / 4
delta_boundary = tau**2 / 4
# 領域の塗り分け
delta_range = np.linspace(-3, 6, 500)
T, D = np.meshgrid(tau, delta_range)
# 色の設定
colors = np.zeros((*T.shape, 4))
# 鞍点: Delta < 0
mask_saddle = D < 0
colors[mask_saddle] = [1, 0.8, 0.8, 0.5]
# 安定結節点: Delta > 0, tau < 0, tau^2 - 4*Delta > 0
mask_stable_node = (D > 0) & (T < 0) & (T**2 - 4*D > 0)
colors[mask_stable_node] = [0.8, 0.8, 1, 0.5]
# 不安定結節点: Delta > 0, tau > 0, tau^2 - 4*Delta > 0
mask_unstable_node = (D > 0) & (T > 0) & (T**2 - 4*D > 0)
colors[mask_unstable_node] = [1, 0.9, 0.7, 0.5]
# 安定渦巻点: Delta > 0, tau < 0, tau^2 - 4*Delta < 0
mask_stable_spiral = (D > 0) & (T < 0) & (T**2 - 4*D < 0)
colors[mask_stable_spiral] = [0.7, 1, 0.7, 0.5]
# 不安定渦巻点: Delta > 0, tau > 0, tau^2 - 4*Delta < 0
mask_unstable_spiral = (D > 0) & (T > 0) & (T**2 - 4*D < 0)
colors[mask_unstable_spiral] = [1, 0.7, 1, 0.5]
# 中心: Delta > 0, tau = 0 (帯で表示)
mask_center = (D > 0) & (np.abs(T) < 0.08)
colors[mask_center] = [0.5, 1, 1, 0.7]
ax.imshow(colors, extent=[-4, 4, -3, 6], origin="lower", aspect="auto")
# 境界線
ax.plot(tau, delta_boundary, "k-", linewidth=2, label="$\\tau^2 - 4\\Delta = 0$")
ax.axhline(0, color="k", linewidth=1.5, linestyle="--")
ax.axvline(0, color="k", linewidth=1, linestyle=":")
# ラベル
ax.text(-3, -1.5, "Saddle", fontsize=14, fontweight="bold", color="darkred")
ax.text(-3, 2, "Stable\nNode", fontsize=12, fontweight="bold", color="darkblue")
ax.text(1.5, 2, "Unstable\nNode", fontsize=12, fontweight="bold", color="darkorange")
ax.text(-1.8, 4.5, "Stable\nSpiral", fontsize=12, fontweight="bold", color="darkgreen")
ax.text(0.5, 4.5, "Unstable\nSpiral", fontsize=12, fontweight="bold", color="purple")
ax.text(-0.5, 5.5, "Center", fontsize=10, fontweight="bold", color="teal", rotation=90)
ax.set_xlabel("Trace $\\tau = \\mathrm{tr}(A)$", fontsize=13)
ax.set_ylabel("Determinant $\\Delta = \\det(A)$", fontsize=13)
ax.set_title("Classification of 2D Linear Systems", fontsize=14)
ax.set_xlim(-4, 4)
ax.set_ylim(-3, 6)
ax.legend(fontsize=11, loc="upper left")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ode_system_trace_det.png", dpi=150, bbox_inches="tight")
plt.show()
このトレース-行列式平面の図は、2次元線形系の完全な分類を一目で把握できるようにしています。
-
水平軸（トレース） が安定性を決めます。$\tau < 0$（左半分）が安定、$\tau > 0$（右半分）が不安定です。これは直感的にも理解できます。トレースは固有値の和であり、これが負なら少なくとも「全体として」減衰する傾向があります
-
$\Delta = 0$ の水平線（破線）が鞍点の境界です。$\Delta < 0$ は固有値が異符号であることを意味し、必ず鞍点になります
-
放物線 $\tau^2 = 4\Delta$ が結節点と渦巻点の境界です。放物線の外側（判別式 > 0）では実固有値、内側（判別式 < 0）では複素固有値になります
-
$\tau = 0$ の縦線上で $\Delta > 0$ の場合が中心（純虚数固有値）です。これは構造的に不安定で、わずかな摂動で渦巻点に変わります
まとめ
本記事では、連立常微分方程式の理論を行列指数関数を中心に解説しました。
- 任意の $n$ 階ODEは連立1階系 $\bm{x}’ = A\bm{x}$ に変換でき、解の理論が統一される
- 固有値・固有ベクトル法では、$A$ の固有値 $\lambda$ と固有ベクトル $\bm{v}$ から解 $\bm{v}e^{\lambda t}$ を構成する
- 行列指数関数 $e^{At}$ は解 $\bm{x}(t) = e^{At}\bm{x}_0$ を与え、対角化またはジョルダン標準形で計算できる
- 2次元系の平衡点は固有値の種類によって鞍点・結節点・渦巻点・中心に分類される
- 非同次系の解は定数変化法により $\bm{x}(t) = e^{At}\bm{x}_0 + \int_0^t e^{A(t-s)}\bm{f}(s)\,ds$ で与えられる
次のステップとして、以下の記事も参考にしてください。
- 常微分方程式の安定性解析 — 非線形系のリアプノフ安定性
- 状態空間モデルの理論 — 制御工学での応用
- 固有値分解の理論 — 行列の対角化の詳細