SPICEの予測子について

---

## はじめに

アナログ回路シミュレータ “**SPICE**” は、回路方程式（非線形な常微分方程式や微分代数方程式）を数値的に「時間積分」して解析します。そのとき重要な役割を果たすのが**予測子 (Predictor)** と**修正子 (Corrector)** という仕組みです。  

数値解析の教科書では、**Adams–Bashforth + Adams–Moulton** のような予測子・修正子法がよく紹介されますが、実際のSPICE実装では**Adams法** は採用されていません。  
代わりに、  
- **Gear法 (BDF)**
- **台形公式 (Trapezoidal rule)**  
の2本柱で構成されているのが一般的です。

本記事では、SPICEの **「予測子」** に焦点を当て、以下のトピックスをまとめます:

1. 予測子・修正子法の概要  
2. Gear法(BDF)の可変次数と予測子の仕組み  
3. Adams法がなぜSPICEで採用されていないのか  
4. 台形公式の予測子例  
5. 高次Gear法(6次)の予測子例  
6. ウォームアップ (warm-up) とは何か  
7. BDF法は1～6次：7次以上は零安定性を失う  
8. **参考情報**: 零-安定性／A-安定性／L-安定性について

---

## 1. 予測子 (Predictor) とは何か

### 1.1 予測子・修正子法の流れ

SPICEの時間積分で非線形要素を含む方程式を解く場合、単なる陽的（明示的）ステップだけでは安定性や誤差面で不十分になることがあります。そこで、

1. **予測子 (Predictor)**  
   - 過去のステップの解・微分値などを用い、多項式外挿等で次時刻の解をざっくり推定。  
2. **修正子 (Corrector)**  
   - 予測子で得た近似解を初期値にして、陰的方程式（Newton法など）で本来の非線形方程式を解き、次ステップの解を確定。

という2段階手法（Predictor–Corrector法）が用いられます。  
これにより、**非線形方程式を解くNewton法**の初期値が改善され、収束速度が向上したり、ステップサイズを大きく取れる可能性が高まります。

---

## 2. Gear法 (BDF) と可変次数予測子

### 2.1 Gear法(BDF法)とは

- **Gear法 (BDF: Backward Difference Formula)** は、  
  過去ステップの解を“後退差分”で表す**陰的多段法**です。  
- 回路方程式のように **“硬い (stiff)”** 性質を持つ問題に対し安定性が高くBerkeley SPICE 以来、多くのSPICE実装で採用されています。

### 2.2 可変次数と予測子

- Berkeley SPICEやngspiceなどでは、**BDF1(1次)～BDF6(6次)** くらいまで可変次数をサポートし、シミュレーション中に「1次→2次→3次…」など**自動的に切り替え**られます。  
- 予測子(Predictor)も**同じ次数**の外挿を行い、修正子(Corrector)ステップでNewton法に渡す流れです。  
- 高次にすると**初期値精度**が上がる一方、ウォームアップが増える・数値振動のリスク・管理コストが増えるなどのデメリットもあるため、**SPICEは誤差評価をしながら“最適な次数”をその都度選択**する仕組みを備えています。

---

## 3. なぜAdams法はSPICEに実装されていないのか？

### 3.1 教科書的には「Adams–Bashforth + Adams–Moulton」

数値解析の一般的な教科書では、  
- **Adams–Bashforth (陽的)** を予測子  
- **Adams–Moulton (陰的)** を修正子  
という形で**予測子・修正子法**をよく解説しています。

### 3.2 SPICEでは採用されていない理由

1. **Gear法 (BDF) で十分に安定性を確保できる**  
   - 回路は硬い系が多く、BDFとの相性が良い。  
2. **台形公式との併用で大抵間に合う**  
   - Berkeley SPICEは “台形公式 or Gear(BDF)” だけで、回路解析の大半をカバーできる。  
3. **陰的Adams (Adams–Moulton) を導入する利点が小さい**  
   - BDF/台形公式 の2本立てで運用可能で、わざわざAdams–Moultonを実装する意義が薄い。

以上の様な理由から、**Adams法** は**SPICE系に実装されていない**（公式ドキュメントやソースなどにも見当たらない）と考えられます。

---

## 4. 台形公式 (Trapezoidal rule) の予測子

### 4.1 台形公式も陰的

- 台形公式は  
  $$
    x_{n+1} = x_n + \frac{\Delta t}{2}\bigl[f(t_n, x_n) + f(t_{n+1}, x_{n+1})\bigr]
  $$  
  という**陰的**方程式となります。  
- Newton法などの反復を使わないと解けないため、**初期値**（予測子）が必要になります。

> **補足**: 実は台形公式は**Adams–Moulton法の1つ (2次)** とも見なせますが、  
> SPICEでは「Trapezoidal rule」と独立した呼称で扱われています。

### 4.2 予測子としてシンプルな外挿

- BDFのような可変次数多段外挿は用いず、**1～2ステップの線形外挿**やオイラー的な方法で“初期値”を作るだけが一般的です。  
- たとえば、  
  $$
    x_{n+1}^{(\mathrm{pred})} = 2x_n - x_{n-1}
  $$  
  などの簡易外挿で初期値を作り、Newton法のスタートにします。

#### 4.2.1 サンプルコード

```cpp
#include <iostream>
#include <cmath>
#include <vector>

// 例: dx/dt = f(t, x) = -x
double f(double t, double x) {
    return -x;
}
double df_dx(double t, double x) {
    return -1.0;
}

int main() {
    double t0 = 0.0, x0 = 1.0;
    double dt = 0.1;
    double tEnd = 1.0;

    // 初期ステップ (例: オイラー法で一回進める)
    double x1 = x0 + dt * f(t0, x0);

    std::vector<double> xHistory(2);
    xHistory[0] = x1;   // x_n
    xHistory[1] = x0;   // x_{n-1}

    double t = t0 + dt;
    double x_n = x1;

    while(t < tEnd - 1e-12) {
        double tNext = t + dt;

        // (A) 予測子: 単なる線形外挿
        double xPred = 2.0 * xHistory[0] - xHistory[1];

        // (B) 修正子: トラペゾイダル(台形公式)をNewton法で解く
        double x_guess = xPred;
        const int MAX_ITER = 20;
        const double TOL = 1e-10;

        for(int iter = 0; iter < MAX_ITER; iter++) {
            double fx_guess = f(tNext, x_guess);
            double R = x_guess - x_n
                       - 0.5*dt*( f(t, x_n) + fx_guess );
            double dRdx = 1.0 - 0.5*dt*df_dx(tNext, x_guess);

            double dx = -R / dRdx;
            x_guess += dx;
            if(std::fabs(dx) < TOL) break;
        }

        double xNext = x_guess;

        std::cout << "t = " << tNext 
                  << ", x = " << xNext << std::endl;

        // 更新
        xHistory[1] = xHistory[0];
        xHistory[0] = xNext;

        t = tNext;
        x_n = xNext;
    }
    return 0;
}
```

---

## 5. 6次Gear法 (BDF6) の予測子例

次に、“高次ギア法”の予測子のイメージとして、サンプルコードを示します。

```cpp
#include <iostream>
#include <vector>
#include <cmath>

//--------------------------------------------------------------
// 6次ギア法(BDF6)の予測子の例
//   xHist[0] = x_n, xHist[1] = x_{n-1}, ..., xHist[5] = x_{n-5}
//--------------------------------------------------------------
double gearPredictor6(const std::vector<double>& xHist) {
    // ここでは binomial係数系の多項式外挿を例示
    double xPred = 6.0*xHist[0] 
                 -15.0*xHist[1]
                 +20.0*xHist[2]
                 -15.0*xHist[3]
                 + 6.0*xHist[4]
                 - 1.0*xHist[5];
    return xPred;
}

int main() {
    // 過去6ステップの値(仮)
    std::vector<double> xHist = {1.0, 0.9, 0.8, 0.6, 0.5, 0.4};

    double xPred = gearPredictor6(xHist);
    std::cout << "Predicted x_{n+1} = " << xPred << std::endl;
    // ここでは予測子のみ。実際にはこれをNewton法の初期値にして修正子ステップを行う。
    return 0;
}
```

- BDF6 では、過去6ステップの解を利用した**多項式外挿**を行い、修正子ステップでBDF6の本体方程式をNewton法で解きます。  
- 高次ほど**ウォームアップ(後述)** や数値振動リスクがあるため、 **SPICEは可変次数** で状況を見ながら最適化します。

---

## 6. ウォームアップ (warm-up) とは何か

### 6.1 多段法の「最初のステップ問題」

多段法 (BDF系、Adams系など) は「過去 n ステップ」の解を利用して次の解を求める仕組みです。たとえばBDF6なら、**過去6ステップ** の解がないと6次の多段外挿はできません。  
しかしシミュレーション開始直後は x(0) だけしかなく、x(1), x(2), ... は未知です。

### 6.2 ウォームアップの必要性

そこで、最初の数ステップは

- もっと低次の手法 (オイラー法やBDF1,2) で進める  
- 台形公式を最初だけ使う  

などして、**「過去 n ステップ分の解データを揃える」** 工程を踏みます。  
これが **「ウォームアップ (warm-up)」** です。

> 例:  
> - BDF6を本格適用したい → まず5ステップ程度はBDF2などで計算して x(1), x(2), ..., x(5) を得る → 6ステップ目からBDF6を開始

### 6.3 ウォームアップのデメリット

- **高次化**すると「n ステップ前」まで必要となり、**ウォームアップ期間**が長くなる。  
- シミュレーション開始時に余計な手順を踏む必要があるため、高次の予測子・修正子が常に得とは限らない要因の一つになります。

---

## 7. BDF法は1～6次まで：7次以上は零-安定性が失われる

### 7.1 零-安定性 (Zero-stability) の問題

- BDF法は**1～6次**までは零-安定性を持ち、**硬い(Stiff)問題にも安定**です。  
- しかし、7次以上にすると零安定性を失うため、数値解が不安定に崩れる可能性が高まることが理論的に示されています（Dahlquistらの結果）。  

### 7.2 SPICE実装への影響

- そのため、SPICE系の実装でも**最大6次**(BDF6)が上限として採用されています。  
- これが **「BDF法は1～6次まで」** という根拠で、Berkeley SPICEやngspiceでもBDF1～BDF6を可変次数でサポートしているわけです。

---

## 8. 参考情報: 零-安定性／A-安定性／L-安定性

**零-安定性 (Zero-stability)** と **A-安定性 (A-stability)**、**L-安定性 (L-stability)** は、いずれも「数値解法の安定性」を論じる概念ですが、扱う範囲や焦点が異なります。

1. **零-安定性** [(Zero-stability)](https://en.wikipedia.org/wiki/Zero_stability)
   - 多段法の係数設定が「ステップ数を増やし、$\(\Delta t\to 0\)$ にしたときに解が破綻しないか」を保証する概念。  
   - BDF7以上で失われるのは、この零-安定性。

2. **A-安定性** (A-stability)  
   - 剛性問題を表すテスト方程式  $\( y' = \lambda y \) (\(\mathrm{Re}(\lambda) < 0\))$ に対し、ステップサイズ無制限に安定解を保てるかを見る概念。  
   - A安定な手法は安定領域が左半平面全体を含む。

3. **L-安定性** [(L-stability)](https://en.wikipedia.org/wiki/L-stability)
   - A-安定性をさらに強化し、大きく負の固有値をもつ場合でも速やかに数値解を減衰させる性質。  
   - 後退オイラー (BDF1と同義) はL安定、台形公式はL安定ではない等の特徴がある。

### 8.1 手法ごとの対応一覧表

以下に、**BDF(1～7)、台形公式、後退オイラー法**に対して、「零-安定性／A-安定性／L-安定性」を満たすかどうかを一覧表にしました。  

| 手法  | 零-安定性 | A-安定性 | L-安定性 |
|------|:-----:|:-----:|:-----:|
| **台形公式 (Trapezoidal rule)** | ○ | ○ | × |
| **後退オイラー (BDF1と同義)**  | ○ | ○ | ○ |
| **BDF2**       | ○ | ○ | × |
| **BDF3〜6**    | ○ | × | × |
| **BDF7〜**     | × | - | - |

> - “○”は満たす、**“×”** は満たさない、**“−”** は零安定性を失うなどで論外。  
> - BDF1は後退オイラーと同義。  
> - 台形公式はA安定だがL安定ではなく、硬い問題で振動を起こしやすい場合もある。

---

## まとめ

1. **SPICEの予測子**  
   - 過去ステップの解を外挿し、陰的解法(修正子)の初期値を作る重要な役割。  
   - BDF(多段法)では可変次数外挿を使い、台形公式ではシンプルな線形外挿を行う。  
2. **Gear法(BDF)がメイン**  
   - 回路問題の硬い特性に適し、可変次数(1～6次)をサポート。  
   - 7次以上は零安定性を失うため実用にならない。  
3. **Adams–Moulton(陰的Adams)** はSPICEに実装されていない  
   - 教科書では「Adams–Bashforth (陽的) + Adams–Moulton (陰的)」が定番だが、  
     SPICEではBDF/台形公式で十分対応可能とされる。  
4. **台形公式**  
   - これも陰的法でNewton法が必要。  
   - 簡易外挿を予測子とすることが多く、SPICEで最も基本的な実装のひとつ。  
5. **ウォームアップ (warm-up)**  
   - 多段法の初期段階で「過去nステップの解」が足りないとき、低次手法等で数ステップ進めて解データを揃える工程。  
   - 高次ほどウォームアップが長くなり、計算コストとのトレードオフを考慮する必要がある。  

---

## 参考リンク

- [ngspice / Berkeley SPICE ドキュメント](http://ngspice.sourceforge.net/)  
- [Backward Difference Formula (Wikipedia, 英語)](https://en.wikipedia.org/wiki/Backward_differentiation_formula)  
- [Predictor–corrector method (Wikipedia, 英語)](https://en.wikipedia.org/wiki/Predictor%E2%80%93corrector_method)  
- Dahlquist, G. (1963), "A special stability problem for linear multistep methods", BIT 3: 27–43.
- [Zero stability (Wikipedia, 英語)](https://en.wikipedia.org/wiki/Zero_stability)
- [L-stability (Wikipedia, 英語)](https://en.wikipedia.org/wiki/L-stability)

以上、「SPICEの予測子」についての総まとめと、**零-安定性／A-安定性／L-安定性**の概念の違いや **BDF(1～6)、BDF7、台形公式、後退オイラー(=BDF1)** の安定性に関する一覧表でした。  

興味のある方は、ngspiceなどのソースを読み解くと、「ステップサイズ制御・誤差評価・Newton法反復」などが複雑に絡み合って、 **可変次数多段法×陰的手法** をうまく動かしている様子がわかります。