SPICE 過渡解析：MNA・数値積分・ニュートン法・Stiff回路への対応

---

## SPICE 過渡解析：MNA・数値積分・ニュートン法・Stiff回路への対応

### はじめに

SPICE（Simulation Program with Integrated Circuit Emphasis）は、抵抗・コンデンサ・インダクタ・ダイオード・トランジスタなど多種多様な電子素子を含む回路をコンピュータ上でシミュレーションするための代表的なソフトウェアです。  
その中でも過渡解析は、時間の経過に沿って回路の各ノード電圧や素子電流がどのように変化していくかを数値的に解く機能です。

過渡解析では、以下のステップを踏むのが基本となります。

1. MNA（Modified Nodal Analysis）方程式で回路を定式化  
2. 数値積分（差分近似）により微分演算を時間離散化  
3. ニュートン法を用いて非線形連立方程式を解く（予測子修正子法やスパースソルバも活用）

さらに、実際の大規模回路では**スティフ（stiff）**な特性を示すことがあり、ギア法（BDF）など安定性に優れた積分法を選ぶことが多い点も重要です。本記事では、こうした観点に加え、ヤコビアン行列をどのように作り、各ニュートン反復でどう更新するのかといった実装イメージを詳しく解説します。

---

## 1. MNA（Modified Nodal Analysis）方程式で回路を定式化

### 1-1. 通常のノード解析との違い
- ノード解析そのものは、ノード電圧を未知数とし、抵抗や電流源をKCL（キルヒホッフの電流則）で表すには便利です。  
- しかし、電圧源やインダクタがあると、ノード電位だけでは未知数が不足したり行列が特異になったりします。

Modified Nodal Analysis (MNA) では、  
1. ノード電圧（グランド以外）  
2. 電圧源・インダクタなどの枝電流  

を未知数に追加し、電圧源やインダクタの振る舞いをスムーズに組み込めるようにします。

#### 勘違いしやすいポイント
「電圧源は電圧が既に決まっているのに、なぜ未知数？」  
→ 実際に未知なのは「電圧源に流れる電流」だからです。MNA ではこの枝電流を未知数に加え、(ノード差 - 電圧源 = 0) という代数方程式を連立に追加します。

### 1-2. 素子スタンピング（Matrix Stamping）
- **抵抗 (R)**: \$i = \frac{v}{R}\$。ノード \$i,\! j\$ 間なら、行列の該当位置に \$\pm \frac{1}{R}\$ を加算  
- **電流源**: KCL の右辺ベクトルに反映  
- **電圧源**: 枝電流を未知数に置き、(ノード電位差 - 電圧源値 = 0) を1行追加  
- **コンデンサ／インダクタ**: DC 解析時はオープン／短絡扱いだが、過渡解析時に後述の数値積分を使って取り込む  
- **非線形素子 (ダイオード, トランジスタ)**: \$i = f(v)\$ をKCLに組み込み、ニュートン法で偏微分を扱う

---

## 2. 過渡解析のための数値積分

### 2-1. 非線形微分代数方程式（DAE）

回路にコンデンサ \$i_C = C \frac{dv}{dt}\$ やインダクタ \$v_L = L \frac{di}{dt}\$ が含まれると、微分演算子＋代数方程式＋非線形が混じった**非線形 DAE (Differential-Algebraic Equation)** となります。  
SPICE の過渡解析では、これを時刻ステップごとに解き進めていきます。

### 2-2. スティフ回路と積分法

大規模回路では非常に異なるタイムコンスタントが共存し、数値的に不安定になりがち（スティフ）です。

- **ギア法 (BDF)**: スティフ系に強く、SPICE2などで主力  
- **トラペゾイド法**: 誤差特性が良いが、振動（ラングリング）を起こす場合あり  
- **後退オイラー法 (Backward Euler)**: 安定性は高いが数値拡散が大きい  

### 2-3. 時間離散化後の非線形方程式

これらの積分法を用いて微分を差分近似すると、時刻 \$t_{n+1}\$ の未知数 \$\mathbf{z}_{n+1}\$ を解く非線形連立方程式  

$$
\boldsymbol{\Phi}(\mathbf{z}_{n+1}) = 0
$$

が生まれます。

- 例: 後退オイラー法なら $i_{C,n+1} = C \frac{v_{C,n+1} - v_{C,n}}{\Delta t}$ など。

---

## 3. ニュートン法（予測子修正子法）とヤコビアン行列

### 3-1. なぜニュートン法？

回路にはダイオードやトランジスタなどの非線形素子が含まれるため、単純な線形代数方程式にはならず、反復解法が必要となります。

- ニュートン法は解の近傍で誤差が二乗で減る（2次収束）特性があり、大規模非線形方程式にも効果的。  
- 一部 SPICE では予測子修正子法を用い、予測ステップと修正ステップを分けて誤差判定やステップ幅制御を行います。

---

## 4. ヤコビアン行列はどう作られ、各反復でどう更新される？

### 4-1. ヤコビアン行列とは？

離散化後の非線形方程式 \$\boldsymbol{\Phi}(\mathbf{z}) = 0\$ に対し、\$ \boldsymbol{\Phi} \$ の各成分を各変数で偏微分した要素を並べた行列が**ヤコビアン行列** \$J(\mathbf{z})\$ です。  
ニュートン法では

$$ J\bigl(\mathbf{z}^{(k)}\bigr)\Delta \mathbf{z} = -\boldsymbol{\Phi}\bigl(\mathbf{z}^{(k)}\bigr) $$ 
  
という線形方程式を解き、\$\Delta \mathbf{z}\$ を得ます。

### 4-2. ヤコビアン行列の作り方（フロー）

1. **MNA方程式 + 数値積分公式** で \$\boldsymbol{\Phi}(\mathbf{z})\$ を定義  
2. **素子ごとに偏微分係数を計算**  
   - 抵抗: \$\frac{\partial i_R}{\partial v_R} = \frac{1}{R}\$  
   - コンデンサ（後退オイラー例）:  

$$
\frac{\partial}{\partial v_{C,n+1}}\Bigl(C \frac{v_{C,n+1}-v_{C,n}}{\Delta t}\Bigr) = \frac{C}{\Delta t}
$$  

   - ダイオード: \$\frac{\partial i_D}{\partial v_D} = f'(v_D)\$  
   - など…  
3. ノード/未知数のインデックスを参照しつつ**行列にスタンプ**  
   - 例: ノード \$i\$ の KCL = \$\Phi_i\$、枝電流 \$j = \mathbf{z}_j\$ なら、(i, j) 要素に \$\frac{\partial \Phi_i}{\partial z_j}\$ を加算。  
4. 推定値 \$\mathbf{z}^{(k)}\$ に依存  
   - ダイオードやコンデンサ等の偏微分係数は \$\mathbf{z}^{(k)}\$ によって変わるので、各反復で再評価し、行列の数値を更新。

### 4-3. 「各反復で一から作り直すわけではない」？

- 構造（行列の形やどの位置が非ゼロか）は、MNA で決まった時点で固定。  
- 各ニュートン反復では、「素子の導関数値」が変わるため、行列要素の数値を再スタンプする。  
- スパースソルバなどを使う場合、行列の**非ゼロパターン**は一定なので、まったくゼロから構造を作り直すわけではなく、数値だけを更新して再因数分解する、という実装が一般的です。

#### 結論
「フルニュートン法」は各反復でヤコビアン行列を「数値的には毎回再計算」するが、行列構造は変わらないので“完全な再構築”ではありません。

---

## 5. DCOP (DC 動作点解析) と最初のステップ

### 5-1. DCOP は「1ステップ目だけ」直接影響

SPICE では、過渡解析前に `.DCOP`（DC動作点解析）を行い、

- コンデンサをオープン・インダクタを短絡として「DC 条件」下での回路解（ノード電圧・枝電流）を求めます。  
- これを最初のステップ (1ステップ目) の初期値に用い、ニュートン法が大きくずれるリスクを抑えます。  
- 2ステップ目以降は前ステップの解を初期値にするため、DCOP は直接効きませんが、最初のステップが安定すると、後続ステップもスムーズに進みます。

### 5-2. ユーザによる初期条件指定も可能

ユーザが「コンデンサ電圧 = 5V」など直接指定する方法もありますが、もし現実と乖離する初期値だと 1ステップ目から発散しかねないため、`.DCOP` による初期動作点の自動算出が多用されます。

---

## 6. ステップサイズの自動制御 (Adaptive Time Step)

### 6-1. なぜステップ制御が必要か

- **大きいステップ**: 計算速度は上がるが急峻な波形に追従できず発散・大誤差の恐れ  
- **小さいステップ**: 精度は上がるが計算時間が増える、丸め誤差も蓄積  
- **Adaptive Time Step**: 収束状況や誤差を見ながらステップ幅を自動調整し、安定性と効率を両立

### 6-2. 予測子修正子法と収束判定

一部の SPICE 実装では、

- 予測子(predictor) で前ステップの解などから次ステップの値を仮推定  
- 修正子(corrector) でニュートン法を回し、誤差を補正  
- 予測値と修正値の差を見て、ステップを増減させる  

といった手段で過渡解析を安定させます。

---

## 7. スティフ回路への対応

### 7-1. スティフ（Stiff）とは

複数の大きく異なるタイムコンスタントが同時存在し、数値的に扱いづらい系を「スティフ（stiff）」と呼びます。

- 例: 超高速トランジスタ部と巨大RC定数部が混在  
- 通常の積分法（トラペゾイド）だと発散・振動を起こす場合があり、ギア法 (BDF) や後退オイラー法が選ばれることが多い。

### 7-2. Gear 法 (BDF) の特徴

- SPICE2 ではギア法が標準。スティフな系でも安定に動作しやすい。  
- ただし波形が丸まる（拡散的になる）傾向があり、ユーザはトラペゾイド法との切り替えを試すこともある。

---

## 8. 特殊ケース：複数解・発振回路・モデル進化

### 8-1. 複数解

非線形回路は複数解を持つ場合あり、ニュートン法は初期値付近の解に収束。別解を狙うには初期条件を変える・ホモトピー法を使うなどが必要。

### 8-2. 発振回路

DCOP では「0V 安定」に見える回路でも、実際には微小ノイズから振動が起こる場合がある。過渡解析を回せば発振波形が再現される。

### 8-3. モデル進化と普遍的フレームワーク

MOSFET やダイオードのモデル（BSIM, EKV など）が高度化しても、「MNA + 数値積分 + ニュートン法 + ヤコビアン構築」の基本枠組みは変わらない。導関数が増えても、行列にスタンプする実装形態は同じ。

---

## 9. まとめ

1. **MNA**: ノード電圧＋電圧源/インダクタの枝電流を未知数とする拡張ノード解析。  
2. **数値積分**: 後退オイラー／トラペゾイド／Gear 法などで微分項を差分近似し、時刻 \$t_{n+1}\$ での非線形連立方程式を作る。  
   - スティフ回路にはギア法(BDF)や後退オイラーが安定。  
3. **ニュートン法（予測子修正子法）**:  
   - 各ステップで生成された非線形方程式を局所線形化、ヤコビアン行列を用いて解を反復更新。  
   - 予測子修正子法の仕組みでステップ幅を調整しながら収束判定を行う。  
4. **ヤコビアン行列の再構築**:  
   - 構造（どこが非ゼロか）は固定だが、素子の偏微分係数は推定値 \$\mathbf{z}^{(k)}\$ に依存 → 毎回“数値”を再スタンプして更新。  
   - スパースソルバでは行列形状を流用し、要素値だけ変えて再因数分解。  
5. **.DCOP と最初のステップ**:  
   - 過渡解析開始前の DCOP で得られた動作点は 1ステップ目の初期値に反映され、大きなズレを防ぐ。  
   - 2ステップ目以降は前ステップの解が初期値となるため、DCOP は直接関与しないが、最初のステップが安定化することで全体の収束が向上。

このように、SPICE の過渡解析は  
「MNAで回路を定式化 → 数値積分で時刻ごとに非線形方程式 → ニュートン法で解を反復更新 → ヤコビアン行列は各素子の偏微分を再スタンプ」  
と進み、スティフな回路にはギア法などで安定性を確保します。  
最初のステップは `.DCOP` 結果による初期値を用いることで、ニュートン法が大きく外れずに収束しやすくなり、結果として全体の解析も円滑に進むのです。

以上