はじめに

数値解析や回路シミュレーションの教科書を見ると、Gear法（BDF法）の安定領域を示す図がしばしば載っていますが、実際にプログラムでどのように描くのかが解説された例はあまりありません。
そこで、Python で 1次～6次Gear法（BDF法）の安定領域を可視化してみるサンプルコードを作成しました。

安定性領域とは

線形多段階法と安定性

常微分方程式の数値解法（線形多段階法）の安定性を議論するために、テスト方程式と呼ばれる以下の簡単な微分方程式を考えます。

$ \displaystyle \frac{dy}{dt} = \lambda y, \quad \lambda \in \mathbb{C} $

ここで、$ \lambda $ は複素数です。線形多段階法をこのテスト方程式に適用したとき、数値解が安定となる条件を調べることで、その解法の安定性を評価できます。

安定性多項式と根の条件

線形多段階法をテスト方程式に適用すると、差分方程式が得られます。この差分方程式の特性多項式が安定性多項式です。安定性多項式の全ての根 $ \zeta_i $ が単位円内（$|\zeta_i|\le1$）にあるとき、その解法は安定になるといえます（根条件 (root condition)）。

安定性領域の定義

安定性領域とは、複素平面上の領域であって、そこに含まれる $ \lambda h $（$ \lambda $ はテスト方程式の固有値、$ h $ はステップ幅）に対して数値解法が安定となる領域のことです。
ふつうは、横軸を Re($ \lambda h $)、縦軸を Im($ \lambda h $) と置き、実部と虚部で平面を描きます。

Gear法と安定性多項式

Gear法（後退差分公式）の概要

Gear法（別名：後退差分公式、BDF: Backward Differentiation Formula）は、過去の複数ステップの値を利用して現在の値を計算する陰的多段法です。
Stiffな常微分方程式に強く、1次Gear法(後退オイラー)や2次Gear法(BDF2)はA-安定として知られています。3次以上(BDF3,4,5,6)になるとA-安定を失いますが、依然として剛性問題で有効に働く場合が多いです。

Gear法の係数とPythonコードでの対応

1次から6次までのGear法において、それぞれの係数 $ \beta_0 $ と $ \alpha_i $ は以下のように定義されます。

次数 ($p$)	$ \beta_0 $	$ \alpha_1 $	$ \alpha_2 $	$ \alpha_3 $	$ \alpha_4 $	$ \alpha_5 $	$ \alpha_6 $
1	1	1
2	2/3	4/3	-1/3
3	6/11	18/11	-9/11	2/11
4	12/25	48/25	-36/25	16/25	-3/25
5	60/137	300/137	-300/137	200/137	-75/137	12/137
6	60/147	360/147	-450/147	400/147	-225/147	72/147	-10/147

(出典: C. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, 1971, Table 8.1 を参照)

Python辞書との対応

後述のPythonコードでは、これらの係数を gear_coefficients という辞書で保持しています。

gear_coefficients = {
    1: {'beta0': 1, 'alpha': [1]},
    2: {'beta0': 2/3, 'alpha': [4/3, -1/3]},
    ...
}

beta0 が $ \beta_0 $ に相当し、alpha が $ \alpha_1, \alpha_2, ... $ に対応します。

安定性多項式

安定性多項式の導出

$p$ 次 Gear 法 (BDF) は、以下の形式で表されます (Gear, 1971)。

$$ y_n = \sum_{i=1}^p \alpha_i y_{n-i} + h \beta_0 y'_n $$

ここで $y_n$ はステップ $n$ での近似解、$y_{n-i}$ は過去のステップでの近似解、$y'_n = f(t_n, y_n)$ です。

この式をテスト方程式 $ y' = \lambda y $ (すなわち $ y'_n = \lambda y_n $) に適用します。

$$ y_n = \sum_{i=1}^p \alpha_i y_{n-i} + h \beta_0 (\lambda y_n) $$

$ z = \lambda h $ とおくと、
$$ y_n = \sum_{i=1}^p \alpha_i y_{n-i} + z \beta_0 y_n $$
この線形定数係数差分方程式の特性解を求めるため、$ y_n = \zeta^n $ と仮定して代入します。
$$ \zeta^n = \sum_{i=1}^p \alpha_i \zeta^{n-i} + z \beta_0 \zeta^n $$
両辺を $ \zeta^{n-p} $ で割り、$ \zeta^p $ の係数を持つ項を左辺にまとめると、
$$ \zeta^p = \sum_{i=1}^p \alpha_i \zeta^{p-i} + z \beta_0 \zeta^p $$
$$ \zeta^p - z \beta_0 \zeta^p - \sum_{i=1}^p \alpha_i \zeta^{p-i} = 0 $$
$$ (1 - \beta_0 z) \zeta^p - (\alpha_1 \zeta^{p-1} + \alpha_2 \zeta^{p-2} + \dots + \alpha_p \zeta^0) = 0 $$
これが $p$ 段 Gear 法の安定性多項式 $ P(\zeta, z) $ となります。

安定性多項式の定義

$ p $段Gear法の安定性多項式 $ P(\zeta, z) $ は、上で導出したように、次のような形で書かれます。（$ z = \lambda h $）
$$
P(\zeta, z) = (1 - \beta_0 z),\zeta^p ;-; \alpha_1 \zeta^{p-1} ;-; \dots ;-; \alpha_p ;=; 0
$$
ここで、$ \zeta $ は特性根（多段法を適用した差分方程式の解）を示し、この根が$|\zeta|\le 1$ であれば安定とみなします（根条件）。
Pythonコードでは、この多項式の係数配列を作り np.roots で根を求め、根がすべて単位円内にあるかどうかで安定か不安定かを判定します。

Pythonコードによる安定性領域の描画

以下に、1次〜6次Gear法（BDF法） それぞれの安定性領域を複素平面上に描画するPythonコードを示します。
ここでは描画範囲を -10〜10（実軸）× -10〜10（虚軸）、メッシュ数を200に設定しています。お好みで変更可能です。

注意: 実行には numpy と matplotlib が必要です。未インストールの場合、
pip install numpy matplotlib

import numpy as np
import matplotlib.pyplot as plt

gear_coefficients = {
    1: {'beta0': 1, 'alpha': [1]},
    2: {'beta0': 2/3, 'alpha': [4/3, -1/3]},
    3: {'beta0': 6/11, 'alpha': [18/11, -9/11, 2/11]},
    4: {'beta0': 12/25, 'alpha': [48/25, -36/25, 16/25, -3/25]},
    5: {'beta0': 60/137, 'alpha': [300/137, -300/137, 200/137, -75/137, 12/137]},
    6: {'beta0': 60/147, 'alpha': [360/147, -450/147, 400/147, -225/147, 72/147, -10/147]}
}

def is_stable(z, p):
    """
    次数 p のGear法において、複素数 z (= λh) が安定領域に含まれるかを判定。
    """
    coeffs = gear_coefficients[p]
    beta0 = coeffs['beta0']
    alpha = coeffs['alpha']

    # 安定性多項式 (1 - beta0*z)*zeta^p - alpha_1 zeta^(p-1) - ... = 0
    # poly_coeffs[0] が zeta^p の係数、[1] が zeta^(p-1)、…、[p] が zeta^0
    poly_coeffs = [(1 - beta0 * z)] + [-a for a in alpha]
    roots = np.roots(poly_coeffs)

    # 全ての根が単位円内かどうか (誤差1e-10加味)
    return np.all(np.abs(roots) <= 1.0 + 1e-10)

def stability_region_gear(p):
    """
    p段Gear法の安定性領域を複素平面に描画。
    """
    # 描画範囲・解像度
    real_axis = np.linspace(-10, 10, 200)
    imag_axis = np.linspace(-10, 10, 200)
    X, Y = np.meshgrid(real_axis, imag_axis)
    Z = X + 1j * Y

    # 安定or不安定をベクトル化した関数で判定
    stable_check = np.vectorize(is_stable, otypes=[bool])
    stable_map = stable_check(Z, p)

    plt.figure(figsize=(6, 6))
    # contourfで True/False を2段階に塗り分け (True=安定=白, False=不安定=灰色)
    plt.contourf(X, Y, stable_map, levels=[0.5, 1], colors=['lightgray', 'white'])
    # 境界線を引く
    plt.contour(X, Y, stable_map, levels=[0.5], colors='black', linewidths=0.8)

    # 凡例を追加 (プロットが見えるようにダミーで作成)
    stable_patch = plt.Rectangle((0, 0), 1, 1, fc='white', ec='black', label='Stable Region')
    unstable_patch = plt.Rectangle((0, 0), 1, 1, fc='lightgray', label='Unstable Region')
    plt.legend(handles=[stable_patch, unstable_patch], loc='upper right', fontsize=9)

    plt.gca().set_aspect('equal', 'box')
    plt.xlim([-10, 10])
    plt.ylim([-10, 10])
    plt.axhline(0, color='k', linewidth=0.8)
    plt.axvline(0, color='k', linewidth=0.8)
    plt.title(f"Stability Region of {p}th-order Gear Method (BDF{p})")
    plt.xlabel("Re(z) = Re(λh)")
    plt.ylabel("Im(z) = Im(λh)")
    plt.grid(True, linestyle='--', alpha=0.6)
    plt.show()

if __name__ == "__main__":
    # 1次~6次まで順番に描画
    for p in range(1,7):
        stability_region_gear(p)

コードの解説

gear_coefficients
- 1〜6次Gear法に対応する $ \beta_0 $ と $ \alpha_i $ をディクショナリに記述。
is_stable(z, p)
- 複素数 $ z = \lambda h $ と次数 $ p $ を引数に取り、安定性多項式を組み立てて np.roots で根を求め、全ての根が $|\zeta|\le1$ かどうかで安定/不安定を判定（根条件のチェック）。
- エッジケースとして、$|\zeta|=1$ となる根を誤差の影響で微妙に超えるかもしれないので、1e-10 程度の許容を加味。
stability_region_gear(p)
- 複素平面のメッシュ（-10〜10, -10〜10に200分割）を用意して is_stable を呼び出し、contourf で塗りつぶし表示。
- レベル levels=[0.5, 1] で二色に分割し、安定域（True）を白、不安定域（False）を灰色で表示。境界線も描画。
- 凡例を追加し、どちらの色が安定領域かを示します。
- グラフのタイトル、軸ラベル、グリッドなどを設定。

実行方法

このコードを gear_stability.py などに保存します。
以下を実行:
```
python gear_stability.py
```
1次〜6次までのGear法について、安定性領域が順次表示されます。

実行時の注意

numpy と matplotlib が必要です。
メッシュ数(200) や描画範囲(-10〜10, -10〜10)は適当に決めています。もっと広く見たい、あるいはもっと細かく見たい場合は、パラメータを調整してください。
3次以上のBDF法ではA-安定を失うため、左半平面の一部に不安定領域が出るのが観察できます。

実行結果

BDF1(後退オイラー) とBDF2 はA-安定なので、Re($z$)<0 の左半平面を完全にカバーする大きな安定領域を持ちます。
BDF3,4,5,6 はA-安定ではないため、左半平面の一部にも不安定領域が生じ、全域をカバーできません。ただし、実務上はこの領域でも大きくステップを取りにくくなるだけで、剛性問題に有効な面も多々あります。

補足説明1: 左半平面と右半平面の安定性

左半平面 (Re(z)<0):
物理的には減衰解に対応する領域。A-安定な数値解法は、この左半平面を丸ごと安定領域に含める。BDF1(後退オイラー)やBDF2がA-安定なのは、この左半平面を完全にカバーするからです。一方、BDF3以上は左半平面の一部が不安定になり、全域をカバーできないためA-安定を失います。
右半平面 (Re(z)>0):
真の解が指数発散する領域。ここが安定か不安定かを数値解法で議論しても、真の連続解自体が安定でないのであまり意味がありません。多くの教科書で「数値解がA-安定か否か」を論じる際は、左半平面をすべて含むかどうかが重要視されます。右半平面については、発散するのが当たり前なので、そこを不安定と呼んでも目新しい話ではないのです。

補足説明2: 教科書や論文によっては図が左右反転しているのはなぜ？

教科書の中には、テスト方程式を $ y' = -\lambda y $ の形で定義したり、$ z = -\lambda h $ などの置換をしている場合があります。その結果、実際には「Re($z$)>0」が減衰解に対応することになり、左右が反転します。
左右反転すると Re($z$)>0 が安定領域になり、その方が直感的なイメージと合うので、あえてそうしているのではないかと思います。

ある文献では
- $ \frac{dy}{dt} = - \lambda y, \quad \text{Re}(\lambda) > 0$
- と定義すると、真の解が減衰するためには $ \text{Re}(\lambda) > 0 $ であり、$z = \lambda h$ とおくと数値解法は Re($z$)>0 側を安定領域に含む必要がある……という形式になります。（この場合、$z = \lambda h$ の定義が文献によって異なる場合もあるので注意が必要です）
そういった符号の取り方やパラメータの定義で、図が左右反転するパターンがあります。
つまり「減衰解に対応する領域」と「発散解に対応する領域」をどちらの符号で定義したかによって、左半平面 or 右半平面が安定領域になるかが変わるのです。

具体例:

電子回路シミュレ-ション (現代非線形科学シリーズ 7), 牛田明夫・田中衞, 2002
- この本の中ではテスト方程式の符号設定やパラメータの導入が異なり、Gear法の安定領域の図が「右半平面」を対象に描かれています。
- 実質的には「減衰解に対応する領域」が右側に来るように符号を取っているだけで、数値的本質は同じです。

結論として、「減衰に対応する領域」が左なのか右なのかは、符号設定次第で逆転しうるため、教科書や論文によって図が左右反転している場合があります。
なお、Gear法の本家本元である、C. Gear氏の著書、“Numerical Initial Value Problems in Ordinary Differential Equations”、は本記事と同じく左半平面を対象としています（左右反転していません）。

まとめ

Pythonで1次〜6次Gear法（BDF法）の安定性領域を可視化する例を示しました。
教科書にある「BDF法はA-安定だ/そうでない」という記述を、実際に複素平面で可視化することで、回路シミュレーションなどでの安定性を直感的に把握できます。
BDF法は陰的多段法であり、剛性問題に強いという利点がありますが、3次以上ではA-安定を失う点を理解しておくことが重要です。
他の多段法（Adams法など）やRunge-Kutta法の安定領域と比較すると、手法ごとの特性の違いがより明確になります。
StiffなODEを実際に数値解し、どの程度大きいステップ幅で安定に解けるか、安定性領域と照らし合わせると面白いです。
教科書や論文によっては符号の取り方で図が左右反転している場合があります。

参考文献

C. Gear, “Numerical Initial Value Problems in Ordinary Differential Equations”, Prentice-Hall, 1971. (特に Chapter 8)
E. Hairer, S. P. Nørsett, G. W. Wanner, “Solving Ordinary Differential Equations I: Nonstiff Problems”, 2nd rev. ed., Springer, 1993. (特に Chapter III.2 Root condition)
E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”, 2nd rev. ed., Springer, 1996. (特に Chapter V.1 BDF methods)
電子回路シミュレ-ション (現代非線形科学シリーズ 7), 牛田明夫・田中衞, 2002.

以上