理論上の数値積分誤差とSPICE実装の違い

---

## はじめに
数値解析の教科書などでは「台形法」や「Gear法(BDF)」といった微分方程式を解くための数値積分法が数多く紹介されています。
特に **Gear法（*Backward Differentiation Formula, BDF*）** は、高次のバックワード（後退）差分を使って精度と安定性を両立しようとする手法で、SPICE系シミュレータに採用されている手法です。

- **後退オイラー法 (Backward Euler)** はBDF1に相当し、そこからさらに高次のBDF (BDF2, BDF3, …) が派生していきます。
- 一方、**台形法（Trapezoidal Rule）** はAdams-Moulton法の1つ（2次）として知られ、BDF（後退差分法）とは別系統の多段法です。

ここで注意したいのは、「理論上の誤差評価」と「実際のSPICE実装の誤差推定」が必ずしも同じ仕組みではない点です。SPICEは予測子-修正子法（Predictor-Corrector）によってローカル誤差を推定・ステップ幅を自動調整しており、汎用的な非線形回路解析に適したアプローチを取っています。本記事では、Gear法などの数値積分法の理論と、SPICEでの実装・誤差推定の違いを概観します。

---

## 1. Gear法(BDF)の理論的背景

### 1.1 微分方程式を数値的に解く考え方

```math
\frac{d\mathbf{x}(t)}{dt} = \mathbf{f}\bigl(\mathbf{x}(t), t\bigr)
```

Gear法は、このような初期値問題を解くための**Backward Differentiation Formula (BDF)** 系列の一種です。
「後退（Backward）差分」を用いて未知の $x_{n+1}$ を右辺に含む（**陰的**）ステップで計算するのが特徴です。

- 例: **BDF1（後退オイラー法）**
  ```math
  x_{n+1} - x_n = h \cdot f\bigl(x_{n+1}, t_{n+1}\bigr)
  ```
  ニュートン法などで $x_{n+1}$ を陰的に解く必要があります。

- **BDF2, BDF3, …**: 高次化によりローカル誤差を理論上は小さくできる一方、安定性や計算コストなどのトレードオフがあります。

### 1.2 理論上の誤差評価（ローカル離散誤差とグローバル誤差）

BDF法の**ローカル離散誤差 (LTE)** は以下のように評価されることが多いです。

```math
\mathrm{LTE} = \|\mathbf{x}(t_{n+1}) - \mathbf{x}_{n+1}^{(true)}\| = O(h^{p+1})
```
（ここで $\mathbf{x}_{n+1}^{(true)}$ は $t_n$ から出発した真の解の $t_{n+1}$ での値）

これが積み重なった**グローバル誤差 (GTE)** は、しばしば最終時刻 $T$ における誤差 $\|\mathbf{x}(T) - \mathbf{x}_N\|$ (ここで $T=Nh$) で評価され、以下のように評価されることが多いです。
```math
\mathrm{GTE} = O(h^p)
```

ここで $p$ はBDFの**次数**（k段BDF法の場合 $p=k$）。
ただし、「高次導関数が有界」「ステップ幅 $h$ が十分小さい」などの仮定が必要で、現実の回路シミュレーションでは破れることも多いです（例: スイッチング回路など）。

---

### 1.3 なぜ陰的解法が安定しやすいのか（数値解析の視点）

**Dahlquistの試験方程式**
```math
\frac{dx}{dt} = \lambda x
```
で $\Re(\lambda)<0$ の剛性がある場合、**陽的オイラー法**では $x_{n+1} = (1+h\lambda)x_n$ となり、ステップ幅が大きいと数値振動・発散を起こしやすい。

**後退オイラー法(BDF1)** は
```math
x_{n+1} = \frac{x_n}{1 - h\lambda}
```
となり、広い範囲の $\lambda$ に対して減衰方向へ進みやすい（A-安定）。
こうして**陰的法**は剛性問題に向いているわけです。

> **数値的ダンピング**
> 陰的法は数値振動が起きにくい反面、物理的高周波成分をやや抑える“ダンピング効果”があります。Gear法が台形法よりRingingしにくい一方、高周波成分を少し消してしまうのはこのためです。

---

### 1.4 BDF法のローカル離散誤差
BDF法の**ローカル離散誤差 (LTE)** は $O(h^{p+1})$ と評価されます（ここで $p$ は次数）。

**k段BDF法** では:

```math
\sum_{j=0}^k \alpha_j\, x_{n+j}
  =
  h\, \beta_k\, f\bigl(x_{n+k},\, t_{n+k}\bigr)
```

という差分方程式を用い、真の解 \(x(t)\) を代入した際の**ズレ**（ローカル離散誤差, LTE）を

```math
\tau_{n+k}
  =
  \frac{1}{h\,\beta_k}
    \sum_{j=0}^k \alpha_j \, x\bigl(t_{n+j}\bigr)
  -
  f\bigl(x(t_{n+k}),\, t_{n+k}\bigr)
```
と定義すると、**$\tau_{n+k} = O(h^{p+1})$** （$p=k$）になるよう係数が設計されています。
こうして、1ステップあたりの理論的な誤差生成が $h^{p+1}$ レベルで、積み重なった**グローバル誤差**が $O(h^p)$ となるわけです。

:::note info
**O記号(Big-O記法)**は、ある関数の**成長率や減少率をざっくり示す**ための表現です。
たとえば「$f(h) = O(h^p)$」と書く場合は、
“$h$を小さくしていったとき、$f(h)$が$h^p$と同程度かそれ以下のオーダーで増減する”
という意味になります。
微分方程式の誤差評価で「$\mathrm{LTE}=O(h^{p+1})$」というときも、**ステップ幅 $h$ が十分小さいときに誤差が $h^{p+1}$ に準じて減るオーダーである** ことを示しています。
:::

---

## 2. SPICEが解く問題の特徴

### 2.1 多変数・非線形の大規模回路

SPICEは抵抗・コンデンサ・トランジスタなど多彩なデバイスモデルが混在する**大規模連立微分方程式**を解きます。
非線形が強く、急峻なスイッチングなどがあると「高次導関数が有界」という仮定が破れ、理論どおりの誤差公式は適用困難です。

### 2.2 シミュレータ内の陰的方程式解法

Gear法(BDF)は**陰的ステップ**をとるので、SPICEは各ステップでニュートン法を回し

```math
\mathbf{x}_{n+1}^{(\mathrm{corr})}
```

を得ます。回路規模が大きいと1ステップの計算が重く、**誤差の上限**を事前に定めるのは現実的ではありません。

> **ポイント**: 次の章で説明する通り、SPICEは各ステップごとに**誤差を “実測”** し、ステップ幅を動的に最適化するという実用的な手段をとっています。

---

## 3. 予測子-修正子法による誤差推定

### 3.1 “実際の近似解同士”の比較で誤差を知る

SPICEはGear法や台形法を下地としつつ、**予測子-修正子法 (Predictor-Corrector)** でローカル誤差を推定します:

1.  **予測 (Predict)**: 多項式外挿などで $\mathbf{x}_{n+1}^{(\mathrm{pred})}$ を算出
2.  **修正 (Correct)**: 予測値を初期とし、ニュートン法を回して $\mathbf{x}_{n+1}^{(\mathrm{corr})}$ を得る
3.  **差分と誤差推定**:
    ```math
    \mathbf{e}_{n+1}
      =
      \mathbf{x}_{n+1}^{(\mathrm{pred})}
      -
      \mathbf{x}_{n+1}^{(\mathrm{corr})}
    ```
    この差分 $\mathbf{e}_{n+1}$ は、次数 $p$ に依存する係数 $C_p$ を用いてローカル離散誤差 (LTE) を推定するためにも使われます（LTE $\approx C_p \mathbf{e}_{n+1}$）。SPICEはこの推定誤差が許容範囲内に収まるようにステップ幅を調整します。誤差が大きければステップ幅を縮め、小さければ拡大して次ステップへ進みます。

### 3.2 導関数を直接評価しなくても済む

理論上はBDF法で高次導関数を仮定し $\mathrm{LTE}=O(h^{p+1})$ を見積もりますが、SPICE実装は「予測解 vs 修正解の差」からLTEを推定するため、回路がどれだけ非線形でも関係なく誤差管理ができるのです。
これによりMOSFETなどの大規模・複雑な回路解析にも実用的に対応できます。

---

## 4. 実際のステップ幅制御とオプション設定

### 4.1 誤差パラメータ (RELTOL, ABSTOL, VNTOL)

- **`.OPTIONS RELTOL=, ABSTOL=, VNTOL=`**
  - **RELTOL**: 相対誤差の許容度
  - **ABSTOL**: 微小電流の絶対誤差
  - **VNTOL**: 微小電圧の絶対誤差

設定が厳しすぎると極端にステップ幅が細かくなり計算が遅い、緩すぎると精度不足。回路に合わせてバランスを取るのが基本です。

### 4.2 BDF法の次数指定 (MAXORD など)

- **`.OPTIONS METHOD=GEAR MAXORD=k`**
  Gear法(BDF)で最大段数k（=最大次数）を制限（NGSPICE等）。
  - 理論的には k段 → k次精度 → $\mathrm{LTE}=O(h^{k+1})$
  - 非線形回路では2～3次が計算コストや安定性との兼ね合いで多く使われます。

### 4.3 METHOD=GEAR と METHOD=TRAP

- **Gear法(BDF)**: 数値ダンピングが大きく、Ringingしにくいが高周波成分をやや抑える。
- **台形法(Trapezoidal)**: 高周波を忠実に再現しやすいが振動しがち。
回路特性・解析目的に応じて使い分けます。

---

## 5. よくある質問 (FAQ)

### Q1. Gear法と台形法の使い分けは？
- **A:** 剛性問題や大規模回路ではGear法が安定しやすい。高速・高周波重視なら台形法が有利の場合もあり。両方試して比較するのが実務的。

### Q2. BDFの高次化で本当に誤差は減るの？
- **A:** 理論的には **LTE=$O(h^{p+1})$** で、次数 $p$（=段数 $k$）が上がるほど精度は良いはずです。しかし、非線形回路で理論の仮定が崩れたり、計算コストが大きくなる、また高次（特にBDF3以降）では安定領域が狭まる問題もあるため、**MAXORD=2～3** が多用されます。

### Q3. 誤差トレランス (RELTOL, ABSTOL, VNTOL) はどう調整すればいい？
- **A:** デフォルトから徐々に調整。厳しすぎると遅く、緩すぎると精度不足。回路特性・解析目的に合わせてバランスを取ります。

### Q4. Predictor-Corrector法の実際の近似解比較と理論的LTE評価はどう違う？
- **A:** 理論では「高次導関数が有界 → $\mathrm{LTE}=O(h^{p+1})$」。SPICE実装は“予測解 vs 修正解の差”からLTEを推定し、回路特有の非線形挙動を直接的に誤差制御に反映できます。

### Q5. Gear法はなぜ数値振動が起きにくい？
- **A:** 後退差分による陰的解法であり、数値的ダンピングが大きいからです。一方、台形法はエネルギー保存性が高く、ピークを過度に再現しがちでRingingを起こす場合があります。

### Q6. BDFと台形法を同時に使うことはできないの？
- **A:** 通常は解析1回につき1つのメイン手法 (Gear or 台形法) を `.OPTIONS METHOD=` で指定。内部切り替え可能な実装もありますが、多くのSPICE系は片方だけを選びます。

### Q7. 剛性が強い回路で台形法は危険？
- **A:** 危険というより、数値振動が大きくステップ幅が極端に小さくなる恐れがあります。Gear法はダンピングが効いて安定収束しやすいです。

### Q8. BDF法の4～6次はあまり使わない？
- **A:** 4次以上になると計算コストが増大するだけでなく、安定性の問題（特にBDF3以降はA-安定性を持たず、安定領域が狭まる）からメリットが見えづらくなります。2～3次が最も実務的とされています。

### Q9. BDF法の誤差評価を“数学的枠組み”と呼ぶのは正しい？
- **A:** 多段法理論で **k段→k次精度→$\mathrm{LTE}=O(h^{k+1})$** が厳密に示されるため、“数学的枠組み”としては正しいです。ただし現実の回路シミュレーションでこの理論通りの誤差挙動を示すかはケースバイケースです。

---

## 6. まとめ

1.  **Gear法(BDF)は理論的にも誤差が評価されている**
    - k段BDFで**ローカル離散誤差(LTE)=$O(h^{k+1})$** → グローバル誤差=$O(h^k)$.
    - 陰的(A-安定性に優れる)手法で剛性問題にも有効。

2.  **SPICEは予測子-修正子法でローカル誤差を推定**
    - 予測解と修正解の差からLTEを推定し、動的にステップ幅を制御しています。
    - 非線形・大規模回路でも実用的に回せる仕組み。

3.  **誤差パラメータと手法選択が要**
    - **RELTOL, ABSTOL, VNTOL**、BDF最大次数(MAXORD)の調整で精度と速度を両立。
    - Gear法(BDF)か台形法かは回路によって一長一短がある。

結果として、Gear法(BDF)の**数学的強み**（$\mathrm{LTE}=O(h^{p+1})$+陰的安定性）と、SPICEの**予測子-修正子**による**実用的な誤差制御**が合わさり、**大規模・非線形な回路を実用的な時間で精度よく**シミュレートできます。

---

## 参考情報: Gear法(BDF)の成り立ちとSPICE2での採用
Gear法(BDF)はなぜStiff問題に使えて、BDF2がA-安定で、$\mathrm{LTE} = O(h^{k+1})$（$k$は段数/次数）となるのか？
なぜそんなに都合よくできているのか? そんな疑問が沸きませんか?

多段法（Linear Multistep Method）の理論に基づき、**k段 → k次精度** になるようテイラー展開の整合性条件を満たし、**ローカル離散誤差(LTE)=$O(h^{k+1})$** を実現しているのがGear法です。また、後退差分を採用することで**剛性問題**にも対応しやすい陰的安定性を得られ、**BDF2** はA-安定であることが証明されています（BDF3以降はA-安定ではありませんが、実用上十分な安定性を持つことが多いです）。

**Gear法**は、まさにこうした仕様（Stiff問題にも安定＆$\mathrm{LTE}=O(h^{k+1})$の高次精度）を満たすよう**意図的に設計された**多段法と言えます。
1971年に**Gear氏**がこの成果を著書で発表し、1970年代中盤にかけて開発された**SPICE2**にも比較的早い段階で採用されたことで、剛性を含む回路解析に活かされることになりました。
つまり、先行研究（Dahlquistら）の多段法理論を剛性問題向けに応用し、**BDF法**を実装レベルで洗練したのが**Gear氏**の貢献なのです。

---

## 参考文献・関連リンク

- G. Dahlquist, *“A Special Stability Problem for Linear Multistep Methods”*, BIT Numerical Mathematics, 3(1), 1963.
- C. Gear, *“Numerical Initial Value Problems in Ordinary Differential Equations”*, Prentice-Hall, 1971.
- E. Hairer, S. P. Nørsett, G. W. Wanner, *“Solving Ordinary Differential Equations I: Nonstiff Problems”*, 2nd rev. ed., Springer, 1993.
- E. Hairer, G. Wanner, *“Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”*, 2nd rev. ed., Springer, 1996.
- K. S. Kundert, *“The Designer's Guide to SPICE and Spectre”*, Springer, 1995.
- [NGSPICE公式ドキュメント](https://ngspice.sourceforge.io/docs.html)

---

以上