振り子を真下に垂らした状態は安定ですが、真上に逆さまにした状態は不安定です。わずかな擾乱で倒れてしまうからです。しかし「安定」「不安定」を直感だけで判断するのには限界があります。変数が3つ以上ある非線形系では、何が安定で何が不安定かは、系統的な数学的手法なしには判断できません。
安定性解析は、微分方程式の解が平衡点の近傍でどのように振る舞うかを調べる理論です。19世紀末にロシアの数学者リアプノフが確立した方法論は、今日でも制御工学、生態学、経済学など幅広い分野で使われています。
安定性解析を理解すると、以下のような場面で活用できます。
- 制御工学: フィードバック制御系の安定判別（リアプノフ安定性）
- 生態学: 捕食者-被食者系の平衡点の安定性（種の共存条件）
- 電気工学: 電力系統の過渡安定度解析
- 機械工学: 非線形振動系のリミットサイクルの存在判定
本記事の内容
- 安定性の数学的定義（リアプノフの意味の安定性）
- 線形化による安定性解析（ハートマン=グロブマンの定理）
- リアプノフの直接法（エネルギー関数法）
- ロトカ=ヴォルテラ方程式の安定性解析
- リミットサイクルとポアンカレ=ベンディクソンの定理
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 連立常微分方程式の理論 — 線形系の位相図
- 固有値分解の理論 — 固有値の計算
- テイラー展開の基礎 — 線形化の基礎
安定性の定義
非線形自律系
一般的な非線形自律系は次のように書けます。
$$ \dot{\bm{x}} = \bm{f}(\bm{x}), \quad \bm{x} \in \mathbb{R}^n $$
ここで $\bm{f}: \mathbb{R}^n \to \mathbb{R}^n$ は十分滑らかな関数です。平衡点（equilibrium point）とは $\bm{f}(\bm{x}^*) = \bm{0}$ を満たす点です。
平衡点に対して、以下の3段階の安定性が定義されます。
リアプノフの意味の安定性（Lyapunov stability）
平衡点 $\bm{x}^*$ が安定（stable in the sense of Lyapunov）であるとは、任意の $\varepsilon > 0$ に対して $\delta > 0$ が存在し
$$ \|\bm{x}(0) – \bm{x}^*\| < \delta \implies \|\bm{x}(t) - \bm{x}^*\| < \varepsilon, \quad \forall t \geq 0 $$
が成り立つことです。直感的には「近くから出発すれば、いつまでも近くにいる」ことを意味します。
漸近安定性（asymptotic stability）
安定であり、かつ
$$ \lim_{t \to \infty} \bm{x}(t) = \bm{x}^* $$
が成り立つとき、漸近安定と言います。「近くから出発すれば、やがて平衡点に収束する」ということです。
不安定性
安定でないとき不安定です。すなわち、どんなに近くから出発しても、ある軌道が $\varepsilon$ 以上離れてしまいます。
安定と漸近安定の違い
線形系の中心（純虚数固有値の場合）は安定ですが漸近安定ではありません。軌道は閉じた楕円を描き、原点から離れることはありませんが、近づくこともありません。一方、安定渦巻点は漸近安定です。
この区別は非線形系で特に重要になります。線形化で中心と判定された場合、非線形項の効果で安定にも不安定にもなりうるからです。
では、安定性を判定するための具体的な手法を見ていきましょう。最初は最も基本的な線形化法です。
線形化による安定性解析
ヤコビ行列と線形化
平衡点 $\bm{x}^*$ の近傍で $\bm{f}(\bm{x})$ をテイラー展開します。$\bm{y} = \bm{x} – \bm{x}^*$ と置くと
$$ \dot{\bm{y}} = \bm{f}(\bm{x}^* + \bm{y}) = \underbrace{\bm{f}(\bm{x}^*)}_{= \bm{0}} + J\bm{y} + O(\|\bm{y}\|^2) $$
ここで $J$ はヤコビ行列（Jacobian matrix）です。
$$ J = \frac{\partial \bm{f}}{\partial \bm{x}}\bigg|_{\bm{x} = \bm{x}^*}, \quad J_{ij} = \frac{\partial f_i}{\partial x_j}\bigg|_{\bm{x} = \bm{x}^*} $$
高次項を無視すると、線形化系 $\dot{\bm{y}} = J\bm{y}$ が得られます。
ハートマン=グロブマンの定理
定理（Hartman-Grobman）: ヤコビ行列 $J$ のすべての固有値の実部が非零（双曲型平衡点）であるとき、非線形系の位相図は線形化系の位相図と定性的に同じ（位相同型）である。
この定理は強力な結果です。双曲型平衡点に限れば、固有値を計算するだけで安定性が完全にわかります。
- $J$ の全固有値の実部が負 → 漸近安定
- $J$ の固有値のうち1つでも実部が正 → 不安定
ただし、実部が0の固有値がある場合（中心の場合）、この定理は適用できません。線形化では安定性を決定できず、高次項の解析が必要になります。これがリアプノフの直接法の出番です。
具体例: 単振り子
減衰のない単振り子の運動方程式は
$$ \ddot{\theta} + \frac{g}{l}\sin\theta = 0 $$
$x_1 = \theta$, $x_2 = \dot{\theta}$ と置くと
$$ \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = -\frac{g}{l}\sin x_1 \end{cases} $$
平衡点1: $(0, 0)$（真下）。ヤコビ行列は
$$ J = \begin{pmatrix} 0 & 1 \\ -\frac{g}{l} & 0 \end{pmatrix} $$
固有値は $\lambda = \pm i\sqrt{g/l}$（純虚数）→ 線形化では中心。ハートマン=グロブマンの定理は適用できませんが、エネルギー保存則からリアプノフの意味で安定であることが示せます（後述）。
平衡点2: $(\pi, 0)$（真上）。$\sin(\pi + y) = -\sin y \approx -y$ なので
$$ J = \begin{pmatrix} 0 & 1 \\ \frac{g}{l} & 0 \end{pmatrix} $$
固有値は $\lambda = \pm\sqrt{g/l}$（実数・異符号）→ 鞍点 → 不安定。
線形化で安定性が決まらない場合に備えて、次はリアプノフの直接法を学びましょう。
リアプノフの直接法
基本的な考え方
物理系のエネルギーが単調に減少するなら、系は安定です。リアプノフの直接法は、この「エネルギー的な考え方」を一般化したものです。
リアプノフ関数 $V(\bm{x})$ とは、以下の性質を持つスカラー関数です。
- $V(\bm{x}^*) = 0$
- $V(\bm{x}) > 0$ for $\bm{x} \neq \bm{x}^*$（正定値）
- $\dot{V}(\bm{x}) = \nabla V \cdot \bm{f}(\bm{x}) \leq 0$（軌道に沿って非増加）
リアプノフの安定性定理
定理: 平衡点 $\bm{x}^*$ の近傍でリアプノフ関数 $V$ が存在するとき
- $\dot{V} \leq 0$ → $\bm{x}^*$ はリアプノフの意味で安定
- $\dot{V} < 0$（$\bm{x} \neq \bm{x}^*$ で狭義負定値） → $\bm{x}^*$ は漸近安定
$\dot{V}$ は軌道に沿った $V$ の時間微分で、連鎖律から
$$ \dot{V} = \sum_{i=1}^n \frac{\partial V}{\partial x_i} f_i(\bm{x}) = \nabla V \cdot \bm{f}(\bm{x}) $$
と計算できます。微分方程式を解かずに安定性を判定できるのが、この方法の大きな利点です。
証明の直感
$V$ は正定値なので、$V = c$（定数）の等高線は平衡点を囲む閉曲線を描きます。$\dot{V} \leq 0$ は、軌道が常に $V$ の等高線の内側に向かう（または等高線上にとどまる）ことを意味します。したがって、$V < \varepsilon^2$ の領域に入った軌道はそこから出ることができず、安定性が保証されます。
単振り子への適用
単振り子のエネルギー
$$ V(\theta, \dot{\theta}) = \frac{1}{2}\dot{\theta}^2 + \frac{g}{l}(1 – \cos\theta) $$
を考えます。$V(0, 0) = 0$, $V > 0$（$(0,0)$ 以外で）を確認できます。
$$ \dot{V} = \dot{\theta}\ddot{\theta} + \frac{g}{l}\sin\theta \cdot \dot{\theta} = \dot{\theta}\left(-\frac{g}{l}\sin\theta\right) + \frac{g}{l}\sin\theta \cdot \dot{\theta} = 0 $$
$\dot{V} = 0$ なので、リアプノフの意味で安定（エネルギーが保存される）。ただし漸近安定ではありません。
減衰振り子の場合
減衰項 $-b\dot{\theta}$ を加えると $\ddot{\theta} + b\dot{\theta} + \frac{g}{l}\sin\theta = 0$ となり
$$ \dot{V} = \dot{\theta}(-b\dot{\theta} – \frac{g}{l}\sin\theta) + \frac{g}{l}\sin\theta \cdot \dot{\theta} = -b\dot{\theta}^2 \leq 0 $$
$\dot{V} \leq 0$ ですが、$\dot{\theta} = 0$ かつ $\theta \neq 0$ の点で $\dot{V} = 0$ となるため、狭義負定値ではありません。ここでLaSalleの不変性原理を使うと、$\dot{V} = 0$ の集合上に留まり続ける軌道は $\dot{\theta} \equiv 0$ すなわち $\ddot{\theta} \equiv 0$ を満たす必要があり、運動方程式から $\sin\theta = 0$ → $\theta = 0$。つまり原点以外に $\dot{V} = 0$ に留まり続ける軌道は存在せず、漸近安定が示されます。
リアプノフ関数の構成
リアプノフの直接法の難しさは、適切なリアプノフ関数を見つけることです。一般的なガイドラインとして
- 物理的エネルギー: 力学系では運動エネルギー + ポテンシャルエネルギーが候補
- 二次形式: $V(\bm{x}) = \bm{x}^T P \bm{x}$（$P$ は正定値対称行列）は線形系の解析に有効
- Krasovskiiの方法: $V = \bm{f}^T \bm{f}$ を試す
線形系 $\dot{\bm{x}} = A\bm{x}$ に対しては、$V = \bm{x}^T P\bm{x}$ として $\dot{V} = \bm{x}^T(A^T P + PA)\bm{x}$ を計算します。$A^TP + PA = -Q$（$Q$ は正定値）を満たす $P$ が存在すれば漸近安定です。この方程式はリアプノフ方程式と呼ばれ、$A$ が安定（全固有値の実部が負）であれば必ず正定値解 $P$ が存在します。
リアプノフの直接法は一般的で強力ですが、実際の応用ではまず線形化を試み、それで決まらない場合にリアプノフ関数を探すのが一般的です。次は、非線形系の具体例としてロトカ=ヴォルテラ方程式の安定性を解析してみましょう。
ロトカ=ヴォルテラ方程式
モデル
ロトカ=ヴォルテラ方程式（捕食者-被食者モデル）は
$$ \begin{cases} \dot{x} = \alpha x – \beta xy \\ \dot{y} = \delta xy – \gamma y \end{cases} $$
ここで $x$ は被食者（ウサギ）、$y$ は捕食者（キツネ）の個体数、$\alpha, \beta, \gamma, \delta > 0$ はパラメータです。
平衡点
$\dot{x} = 0, \dot{y} = 0$ を解くと、2つの平衡点が得られます。
- $(x^*, y^*) = (0, 0)$ — 両種絶滅
- $(x^*, y^*) = (\gamma/\delta, \alpha/\beta)$ — 共存
各平衡点の安定性解析
ヤコビ行列:
$$ J = \begin{pmatrix} \alpha – \beta y & -\beta x \\ \delta y & \delta x – \gamma \end{pmatrix} $$
平衡点 $(0, 0)$:
$$ J(0,0) = \begin{pmatrix} \alpha & 0 \\ 0 & -\gamma \end{pmatrix} $$
固有値は $\lambda_1 = \alpha > 0$, $\lambda_2 = -\gamma < 0$ → 鞍点（不安定）。被食者は増殖し、捕食者は餓死する方向に動く。
平衡点 $(\gamma/\delta, \alpha/\beta)$:
$$ J = \begin{pmatrix} 0 & -\beta\gamma/\delta \\ \delta\alpha/\beta & 0 \end{pmatrix} $$
固有値は $\lambda = \pm i\sqrt{\alpha\gamma}$（純虚数）→ 線形化では中心。ハートマン=グロブマンの定理は適用できません。
リアプノフ関数による安定性
ロトカ=ヴォルテラ方程式には保存量
$$ V(x, y) = \delta x – \gamma \ln x + \beta y – \alpha \ln y $$
が存在します。$\dot{V} = 0$ を直接計算で確認できます。
$\dot{V}$ を計算すると
$$ \dot{V} = \delta \dot{x} – \frac{\gamma}{x}\dot{x} + \beta \dot{y} – \frac{\alpha}{y}\dot{y} $$
各項に $\dot{x} = \alpha x – \beta xy$, $\dot{y} = \delta xy – \gamma y$ を代入すると
$$ \dot{V} = \delta(\alpha x – \beta xy) – \gamma(\alpha – \beta y) + \beta(\delta xy – \gamma y) – \alpha(\delta x – \gamma) $$
展開して整理すると
$$ \dot{V} = \delta\alpha x – \delta\beta xy – \gamma\alpha + \gamma\beta y + \beta\delta xy – \beta\gamma y – \alpha\delta x + \alpha\gamma = 0 $$
すべての項が完全に相殺されます。$V$ は保存量であり、共存平衡点の近傍でリアプノフ関数の条件を満たします（平衡点で最小値をとる）。$\dot{V} = 0$ なので、共存平衡点はリアプノフの意味で安定ですが、漸近安定ではありません。
軌道は $V = \text{const.}$ の等高線上を周期的に巡ります。これが「ウサギが増えるとキツネも増え、キツネが増えるとウサギが減り…」という個体数の周期的変動を数学的に説明しています。
次は、もう一つの重要な現象であるリミットサイクルについて見てみましょう。
リミットサイクルとファンデルポル振動子
リミットサイクルとは
リミットサイクル（limit cycle）とは、位相平面上の孤立した閉軌道のことです。近傍の軌道がリミットサイクルに漸近する場合、安定リミットサイクルと呼びます。
ロトカ=ヴォルテラ方程式の閉軌道は保存系による周期軌道であり、リミットサイクルではありません（孤立していないため）。リミットサイクルは散逸系に特有の現象で、エネルギーの供給と散逸のバランスから自発的に周期運動が生まれます。
ファンデルポル振動子
ファンデルポル振動子は
$$ \ddot{x} – \mu(1 – x^2)\dot{x} + x = 0 \quad (\mu > 0) $$
$|x| < 1$ のとき減衰が負（エネルギーが供給される）、$|x| > 1$ のとき減衰が正（エネルギーが散逸する）。このバランスにより、振幅に依存しない自励振動が発生します。
1階系に変換すると
$$ \begin{cases} \dot{x}_1 = x_2 \\ \dot{x}_2 = \mu(1 – x_1^2)x_2 – x_1 \end{cases} $$
原点のヤコビ行列は
$$ J(0,0) = \begin{pmatrix} 0 & 1 \\ -1 & \mu \end{pmatrix} $$
固有値は $\lambda = \frac{\mu}{2} \pm \frac{1}{2}\sqrt{\mu^2 – 4}$ です。$\mu > 0$ のとき固有値の実部は正なので、原点は不安定です。
ポアンカレ=ベンディクソンの定理
定理（Poincare-Bendixson）: $\mathbb{R}^2$ の有界な正不変集合が平衡点を含まないならば、その中にリミットサイクルが存在する。
ファンデルポル振動子では、原点が不安定（軌道が外に向かう）かつ十分遠くでは軌道が内に向かうことを示せるため、この定理からリミットサイクルの存在が保証されます。
この定理は2次元系にのみ成り立つことに注意してください。3次元以上ではカオスが可能になるため、この定理の拡張は成立しません。
Pythonによる実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
fig, axes = plt.subplots(2, 2, figsize=(14, 12))
# --- (a) 単振り子の位相図 ---
ax = axes[0, 0]
g_over_l = 1.0
def pendulum(t, state):
theta, omega = state
return [omega, -g_over_l * np.sin(theta)]
for E_level in np.linspace(0.1, 3.5, 12):
# 初期条件: theta=0, omega = sqrt(2E) (エネルギーレベル)
omega0 = np.sqrt(2 * E_level)
if omega0 > 4:
continue
for sign in [1, -1]:
sol = solve_ivp(pendulum, [0, 20], [0, sign * omega0],
t_eval=np.linspace(0, 20, 2000), rtol=1e-10)
ax.plot(sol.y[0], sol.y[1], "b-", linewidth=0.6, alpha=0.6)
# セパラトリクス
omega_sep = np.sqrt(2 * g_over_l * 2)
for sign in [1, -1]:
sol = solve_ivp(pendulum, [0, 50], [0.001, sign * omega_sep * 0.9999],
t_eval=np.linspace(0, 50, 5000), rtol=1e-12)
mask = np.abs(sol.y[0]) < np.pi * 1.2
ax.plot(sol.y[0][mask], sol.y[1][mask], "r-", linewidth=1.5, alpha=0.8)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.plot(0, 0, "go", markersize=8)
ax.plot(np.pi, 0, "rx", markersize=10, markeredgewidth=2)
ax.plot(-np.pi, 0, "rx", markersize=10, markeredgewidth=2)
ax.set_xlabel(r"$\theta$", fontsize=12)
ax.set_ylabel(r"$\dot{\theta}$", fontsize=12)
ax.set_title("(a) Simple Pendulum Phase Portrait", fontsize=13)
ax.grid(True, alpha=0.3)
# --- (b) ロトカ=ヴォルテラ ---
ax = axes[0, 1]
alpha_lv, beta_lv, gamma_lv, delta_lv = 1.0, 0.5, 0.8, 0.3
def lotka_volterra(t, state):
x, y = state
return [alpha_lv * x - beta_lv * x * y,
delta_lv * x * y - gamma_lv * y]
# 共存平衡点
x_eq = gamma_lv / delta_lv
y_eq = alpha_lv / beta_lv
for r in [0.3, 0.6, 1.0, 1.5, 2.0]:
for angle in np.linspace(0, 2 * np.pi, 4, endpoint=False):
x0 = x_eq + r * np.cos(angle)
y0 = y_eq + r * np.sin(angle)
if x0 > 0.1 and y0 > 0.1:
sol = solve_ivp(lotka_volterra, [0, 30], [x0, y0],
t_eval=np.linspace(0, 30, 3000), rtol=1e-10)
ax.plot(sol.y[0], sol.y[1], "b-", linewidth=0.6, alpha=0.6)
ax.plot(x_eq, y_eq, "go", markersize=8, label=f"Equilibrium ({x_eq:.1f}, {y_eq:.1f})")
ax.set_xlabel("Prey $x$", fontsize=12)
ax.set_ylabel("Predator $y$", fontsize=12)
ax.set_title("(b) Lotka-Volterra Phase Portrait", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 8)
ax.set_ylim(0, 6)
# --- (c) ファンデルポル振動子（mu=1） ---
ax = axes[1, 0]
mu_vdp = 1.0
def van_der_pol(t, state):
x, v = state
return [v, mu_vdp * (1 - x**2) * v - x]
# 様々な初期条件
for r0 in [0.3, 0.5, 1.0, 3.0, 4.0]:
for angle in np.linspace(0, 2 * np.pi, 4, endpoint=False):
x0 = r0 * np.cos(angle)
v0 = r0 * np.sin(angle)
sol = solve_ivp(van_der_pol, [0, 30], [x0, v0],
t_eval=np.linspace(0, 30, 3000), rtol=1e-10)
ax.plot(sol.y[0], sol.y[1], "b-", linewidth=0.6, alpha=0.5)
# リミットサイクル（十分な時間後）
sol_lc = solve_ivp(van_der_pol, [0, 100], [0.1, 0],
t_eval=np.linspace(80, 100, 2000), rtol=1e-10)
ax.plot(sol_lc.y[0], sol_lc.y[1], "r-", linewidth=2.5, label="Limit cycle")
ax.plot(0, 0, "kx", markersize=10, markeredgewidth=2)
ax.set_xlabel("$x$", fontsize=12)
ax.set_ylabel("$\\dot{x}$", fontsize=12)
ax.set_title(f"(c) Van der Pol Oscillator ($\\mu={mu_vdp}$)", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 5)
# --- (d) ファンデルポルの時間波形 ---
ax = axes[1, 1]
for mu_val, color, ls in [(0.5, "blue", "-"), (1.0, "green", "--"), (3.0, "red", ":")]:
def vdp_mu(t, state, mu=mu_val):
x, v = state
return [v, mu * (1 - x**2) * v - x]
sol = solve_ivp(vdp_mu, [0, 40], [0.1, 0],
t_eval=np.linspace(0, 40, 2000), rtol=1e-10)
ax.plot(sol.t, sol.y[0], color=color, linestyle=ls, linewidth=1.5,
label=f"$\\mu = {mu_val}$")
ax.set_xlabel("$t$", fontsize=12)
ax.set_ylabel("$x(t)$", fontsize=12)
ax.set_title("(d) Van der Pol Time Series", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ode_stability_analysis.png", dpi=150, bbox_inches="tight")
plt.show()
この4つの図から、非線形系の安定性に関する重要な知見が読み取れます。
-
左上（単振り子）: 原点（緑丸、真下の平衡点）の周りの軌道は閉曲線を描いており、リアプノフの意味で安定（中心）です。一方、$(\pm\pi, 0)$（赤×、真上の平衡点）は鞍点で、セパラトリクス（赤い曲線）が安定多様体と不安定多様体を形成しています
-
右上（ロトカ=ヴォルテラ）: 共存平衡点の周りの軌道は閉曲線で、保存量 $V$ の等高線上を周期的に巡っています。すべての軌道が閉じていることは、このモデルが保存系であることの直接的な反映です
-
左下（ファンデルポル振動子）: 原点（黒×）は不安定で、内側の軌道は外に向かいます。一方、外側の軌道は内に向かいます。両者がリミットサイクル（赤い太線）に漸近しており、安定リミットサイクルの存在が確認できます
-
右下（ファンデルポル時間波形）: $\mu$ が大きくなると波形が正弦波から離れ、緩和振動（急激な変化とゆっくりした変化の繰り返し）になっていきます。$\mu = 0.5$ ではほぼ正弦波、$\mu = 3.0$ では明確な緩和振動が見られます
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import solve_lyapunov
# リアプノフ方程式の解と安定性判定
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# --- (a) リアプノフ関数の等高線（2次元線形系） ---
ax = axes[0]
A = np.array([[-1, 2], [-3, -2]])
# リアプノフ方程式 A^T P + P A = -Q を解く
Q = np.eye(2)
P = solve_lyapunov(A.T, -Q)
print(f"A の固有値: {np.linalg.eigvals(A)}")
print(f"P = \n{P}")
print(f"P の固有値: {np.linalg.eigvals(P)} (全て正なら P は正定値)")
x_range = np.linspace(-3, 3, 200)
y_range = np.linspace(-3, 3, 200)
X, Y = np.meshgrid(x_range, y_range)
# V(x) = x^T P x
V = P[0, 0] * X**2 + (P[0, 1] + P[1, 0]) * X * Y + P[1, 1] * Y**2
# V_dot = -x^T Q x = -(x1^2 + x2^2)
V_dot = -(X**2 + Y**2)
ax.contour(X, Y, V, levels=np.linspace(0.1, 8, 15), cmap="Blues", linewidths=1)
ax.contour(X, Y, V_dot, levels=np.linspace(-8, -0.1, 10),
cmap="Reds_r", linewidths=0.8, linestyles="--")
# 軌道
from scipy.integrate import solve_ivp
for angle in np.linspace(0, 2 * np.pi, 8, endpoint=False):
x0 = [2.5 * np.cos(angle), 2.5 * np.sin(angle)]
sol = solve_ivp(lambda t, x: A @ x, [0, 5], x0,
t_eval=np.linspace(0, 5, 500), rtol=1e-10)
ax.plot(sol.y[0], sol.y[1], "k-", linewidth=0.8, alpha=0.5)
ax.plot(0, 0, "go", markersize=8)
ax.set_xlabel("$x_1$", fontsize=12)
ax.set_ylabel("$x_2$", fontsize=12)
ax.set_title("(a) Lyapunov Function $V = \\mathbf{x}^T P \\mathbf{x}$\n(blue) and $\\dot{V}$ (red dashed)", fontsize=12)
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
# --- (b) リアプノフ関数としてのエネルギー（振り子） ---
ax = axes[1]
theta_range = np.linspace(-np.pi, np.pi, 200)
omega_range = np.linspace(-3, 3, 200)
THETA, OMEGA = np.meshgrid(theta_range, omega_range)
# エネルギー V = 0.5*omega^2 + (1 - cos(theta))
V_pendulum = 0.5 * OMEGA**2 + (1 - np.cos(THETA))
ax.contourf(THETA, OMEGA, V_pendulum, levels=20, cmap="YlOrRd", alpha=0.7)
cs = ax.contour(THETA, OMEGA, V_pendulum, levels=[0.5, 1.0, 1.5, 2.0],
colors="black", linewidths=1)
ax.clabel(cs, fmt="%.1f", fontsize=9)
# 軌道
for E in [0.3, 0.7, 1.2, 1.8]:
omega0 = np.sqrt(2 * E)
sol = solve_ivp(lambda t, s: [s[1], -np.sin(s[0])], [0, 20], [0, omega0],
t_eval=np.linspace(0, 20, 2000), rtol=1e-10)
ax.plot(sol.y[0], sol.y[1], "w-", linewidth=1.5, alpha=0.8)
ax.plot(0, 0, "go", markersize=8)
ax.set_xlabel(r"$\theta$", fontsize=12)
ax.set_ylabel(r"$\dot{\theta}$", fontsize=12)
ax.set_title("(b) Pendulum Energy as Lyapunov Function", fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("ode_lyapunov_functions.png", dpi=150, bbox_inches="tight")
plt.show()
この図から、リアプノフ関数の性質が視覚的に理解できます。
-
左図（二次形式リアプノフ関数）: 青い等高線は $V = \bm{x}^T P\bm{x}$ の値を示しており、楕円形です。赤い破線は $\dot{V} = -\bm{x}^T Q\bm{x}$ の等高線で、常に負（原点以外で）です。軌道（黒線）が $V$ の等高線を内側に横切りながら原点に収束していることが確認できます
-
右図（振り子のエネルギー）: 色の濃さがエネルギー $V$ の大きさを表しています。白い軌道は $V = \text{const.}$ の等高線上を動いており、エネルギーが保存されていることが直接確認できます。原点（緑丸）はエネルギーの最小値をとる点で、リアプノフ関数の条件を満たしています
まとめ
本記事では、常微分方程式の安定性を解析する3つの主要な手法を解説しました。
- 安定性の3段階は、リアプノフの意味の安定性（近くにとどまる）、漸近安定性（収束する）、不安定性で定義される
- 線形化法はヤコビ行列の固有値で安定性を判定する。ハートマン=グロブマンの定理により、双曲型平衡点では非線形系と線形化系の位相図が定性的に一致する
- リアプノフの直接法は微分方程式を解かずに安定性を判定でき、物理的にはエネルギー関数の単調減少に対応する
- ロトカ=ヴォルテラ方程式の共存平衡点は保存量の存在によりリアプノフ安定だが漸近安定ではない
- リミットサイクルは散逸系に特有の自励振動で、ポアンカレ=ベンディクソンの定理により2次元系での存在が保証される
次のステップとして、以下の記事も参考にしてください。
- 連立常微分方程式の理論 — 線形系の解法と位相図
- リアプノフ安定性（制御工学） — 制御系設計への応用
- 常微分方程式の数値解法 — ルンゲ=クッタ法