はじめに

本記事は、前回の 「SPICEの過渡解析における誤差推定とステップ制御」 の続編です。前回は主に 台形法（Trapezoidal） と BDF（Gear）法 に触れましたが、じつは 両者を1ステップ内で組み合わせた"TR‑BDF2法" という選択肢もあります。そこで今回は、その TR‑BDF2法 についてまとめてみます。
（なお、MATLAB/Simulink では ode23tbというソルバー名称で実装されています）

TR‑BDF2法の特徴は次のとおりです。

台形とBDF2の"いいとこ取り"ができる
1ステップを 台形ステージ → BDF2ステージ の2段で進めることで、台形の2次精度と、BDF系の強い安定性（$L$安定）を両立できます。特に $ \gamma = 2-\sqrt{2} \approx 0.5858 $ を用いる設計が実装上では扱いやすいです。
SPICE流の実装に載せやすい
両ステージとも残差を
$$
R(\mathbf{x})=\mathbf{f}(\mathbf{x},t_\star)+\alpha\bigl(\mathbf{q}(\mathbf{x})-\hat{\mathbf{q}}\bigr)=\mathbf{0}
$$
の形にでき、ニュートン反復のヤコビアンは $ \mathbf{J}=\mathbf{G}+\alpha\mathbf{C} $ になります。さらに $ \gamma=2-\sqrt{2} $ では 両ステージの $\alpha$ が一致し、同じ行列 $ \mathbf{G}+\alpha\mathbf{C} $ のLU分解を再利用できるため、計算コストや収束の安定性でメリットが出やすいです。
"台形の数値リンギング"と"Gearの過減衰"の折衷として有効
スイッチングや半導体非線形が強い回路で出やすい数値リンギング（台形）と、位相・波形が重くなりがちな過減衰（Gear）の間を、現実的にバランスできます。
TCADでは実績があるが、SPICE系では未普及
デバイス方程式のような強剛性問題では TR‑BDF2が実務採用されています。一方で、一般的なSPICE系メジャーツールでは trap/gear中心でTR‑BDF2は未採用で、このギャップの背景（イベント処理や既存資産との親和性など）を整理する価値がありそうです。

本記事では、まず、フローチャートで、TR‑BDF2法の実装フローを可視化します。さらに、TCADでは使われているがSPICE系では一般的でない理由、$ \gamma $ の選び方と 同一行列再利用、誤差推定・ステップ制御 の勘所まで、前回の文脈に沿って実装者目線で解説します。

sukimaspice での実装状況（2025-11-03 更新）

✅ TR-BDF2 完全版実装済み（Phase 8.8）

sukimaspice では、本記事で解説する γ = 2 - √2 による LU分解再利用最適化を完全に実装しました：

シンボリック分解（analyze）を1回のみ実行: 両ステージで同じ反復行列 $ \mathbf{G}+\alpha\mathbf{C} $ を使用

Stage 1（台形）: 最初のNewton反復で analyze_and_factorize()、以降は refactorize()

Stage 2（BDF2）: 全Newton反復で refactorize()（Stage 1のシンボリック分解を再利用）

1. TR‑BDF2は「台形」と「BDF2」を両方解く

TR‑BDF2（Trapezoidal Rule–BDF2）は 1ステップの内部で

台形（Crank–Nicolson, TR）の前段ステージと、
BDF2 の後段ステージ

を 連続して 解く 二段一体 の方法です。$ \gamma \in (0,1) $ を段間位置として、時刻 $t_n \to t_{n+\gamma} \to t_{n+1}$ の順に進めます。実装上は 両段で同じ反復行列（$G+\alpha C$）を使い回せる設計にでき、強剛性でも安定に動作します。

イメージ図

時間 t --->

t_n                           t_{n+γ}                     t_{n+1}
 |<------------------------------|-------------------------->|
 |<-----------  γ h  ----------->|<-------  (1−γ) h  ------->|
 |              TRZ              |           BDF2            |

要点：区間幅 h を「前半 γh（TRZ）」「後半 (1−γ)h（BDF2）」に分割します。さらに

$$
\frac{2}{\gamma} = \frac{2-\gamma}{1-\gamma}
$$

となるように $\gamma$ を選ぶと（解は $\gamma = 2-\sqrt{2} \approx 0.5858$）、両段の係数が揃って 同じ反復行列 $G+\alpha C$ のLU分解を再利用できます。

2. 数式でざっくり（MNAのDAEに対して）

回路のMNAはしばしば
$$
\frac{d\mathbf{q}(\mathbf{x})}{dt} + \mathbf{f}(\mathbf{x},t) = \mathbf{s}(t)
$$
のようなDAEで表せます（$\mathbf{x}$は状態，$\mathbf{q}$は蓄積量，$\mathbf{f}$は抵抗性素子や非線形素子の電流写像）です。ステップ幅を $h$，$ \gamma = 2-\sqrt{2}$ とします。

TRステージ（$t_n \rightarrow t_{n+\gamma}$）
台形則を長さ $\gamma h$ で適用し、非線形方程式

$$
R^{\text{TR}}(\mathbf{x}{n+\gamma})= \mathbf{f}(\mathbf{x}{n+\gamma}, t_{n+\gamma})+ \alpha \bigl(\mathbf{q}(\mathbf{x}_{n+\gamma}) - \widehat{\mathbf{q}}1\bigr)- \mathbf{s}(t{n+\gamma}) = \mathbf{0},
$$
( ここ、なぜかうまく数式が表示されない ... 原因不明 )
をニュートンで解きます。ここで $ \alpha = \dfrac{2}{\gamma h}$，$\widehat{\mathbf{q}}_1$ は既知量から組む仮想右辺です。

BDF2ステージ（$t_{n+\gamma} \rightarrow t_{n+1}$）
長さ $(1-\gamma)h$ のBDF2を適用し，
$$
R^{\text{BDF2}}(\mathbf{x}{n+1})
= \mathbf{f}(\mathbf{x}{n+1}, t_{n+1})
- \alpha' \bigl(\mathbf{q}(\mathbf{x}_{n+1}) - \widehat{\mathbf{q}}_2\bigr)
- \mathbf{s}(t_{n+1}) = \mathbf{0},
  $$
  ( ここも、なぜかうまく数式が表示されない ... 原因不明 )
  を解きます。ここで $ \alpha' = \dfrac{2-\gamma}{(1-\gamma)h}$ です。

$ \gamma = 2-\sqrt{2} $ とすると
$$
\alpha = \alpha' = \frac{2+\sqrt{2}}{h}
$$
となり、両ステージで同一の係数行列 $G + \alpha C$ を用いる実装（同一LUの再利用）が可能です。ステップ受理・棄却は、TR‑BDF2本体（2次）と埋め込み高次近似（3次）との差を用いた 埋め込み誤差推定 が定番です。

3. 実装フロー

TR‑BDF2 を 回路DAE に適用したときの、1ステップあたりの代表的なフローです（ブレークポイント処理・ニュートン反復・LU再利用を含めて図示します）。

実務上のポイント

同一係数行列の再利用：$ \gamma = 2-\sqrt{2} $ なら TR 段と BDF2 段で $G+\alpha C$ が一致し、1回のLU分解で両段に使い回しできます。
イベント（ブレークポイント）：区間内にステップ境界が存在する場合、ステップ短縮で内部ステージ $t_{n+\gamma}$ が境界を跨がないように調整します。
誤差推定とステップ制御：$p=2$ として $h_{\text{new}}=\mathrm{safety}\cdot h\cdot E^{-1/(p+1)}=\mathrm{safety}\cdot h\cdot E^{-1/3}$ を基本にします（上下限ガードを入れると安定します）。

4. どのツールが TR‑BDF2 を採用しているか

TCAD（半導体デバイスシミュレータ）
公開資料や論文で、Synopsys Sentaurus Device や Silvaco ATLAS などに TR‑BDF2 の採用記述が見られます。強剛性の連立PDEを安定に積分する目的で、既定あるいは主要選択肢として用いられる構成が一般的です。
一般のSPICE系回路シミュレータ
Cadence Spectre／Xyce／ngspice／LTspice などの主流では、積分法は 台形（trap）／Gear（BDF）中心 です。ヘルプやユーザーガイドの選択肢に TR‑BDF2 は無く、公式サポートはされていないようです。

研究用・自作系では、MATLAB/Simulink の ode23tb（TR‑BDF2実装） や Julia の DifferentialEquations.jl TRBDF2() を MNA に適用して試すことが現実的です。

5. なぜ TCAD では使われ、SPICE系では一般的でないのか（勝手な考察）

問題構造の違い（剛性とPDE/DAEの規模）
TCAD は連立PDE（連続の式＋ポアソン等）由来の強剛性で、$L$安定かつ 1ステップ内2段で誤差と安定性を両立するTR‑BDF2のメリットが大きいと考えられます。
行列分解の再利用が効きやすい
$ \gamma=2-\sqrt{2} $ で 両段の反復行列が同一になるため、1回のLU分解を2段で再利用できます。大規模3Dメッシュではこの特性が効率に直結します。
SPICEのイベント密集・離散切替との相性
スイッチングや理想素子のオン/オフでイベント境界が多い回路では、内部ステージ $t_{n+\gamma}$ が境界を跨がないよう細やかなステップ分割が必要です。TR‑BDF2は各ステップで2段解くため、イベント処理の分岐と履歴管理が複雑になりがちです。長年最適化された trap/gear＋ブレークポイント制御 の枠組みの方が実装資産を活かせます。
既存フロー資産とユーザ設定の簡潔さ
現場のフローやモデル・パラメータは trap/gear 前提 に最適化されていることが多く、「選択肢を増やすコスト＞得られる利得」になりがちです。LTspice の modified‑trap のように、台形の欠点（数値リンギング）を低コストで抑える改良も普及しました。

6. 係数の要点（ $ \gamma = 2-\sqrt{2} $ の効果）

代表的パラメータ選択は
$$
\gamma = 2-\sqrt{2} \approx 0.5858
$$
です。これにより

$L$安定（高周波成分の数値的減衰が強く、剛性問題に有利）
TR段/BDF2段で同一反復行列（$G+\alpha C$）
$$
\alpha = \frac{2}{\gamma h} = \frac{2-\gamma}{(1-\gamma)h} = \frac{2+\sqrt{2}}{h}
$$
が成り立ち、1回のLU分解を2段で再利用できます。

7. 参考：TR‑BDF2 の段方程式（常微分形のイメージ）

常微分形 $ \dot{u}=f(u,t) $ への適用例（段方程式の一例）です。

TR段（$t_n\to t_{n+\gamma}$）
$$
u_{n+\gamma} - \frac{\gamma h}{2} f(u_{n+\gamma},t_{n+\gamma})
= u_n + \frac{\gamma h}{2} f(u_n,t_n)
$$
BDF2段（$t_{n+\gamma}\to t_{n+1}$）
$$
u_{n+1} - \frac{(1-\gamma)h}{2-\gamma} f(u_{n+1},t_{n+1})
= \frac{1}{\gamma(2-\gamma)}u_{n+\gamma} - \frac{(1-\gamma)^2}{\gamma(2-\gamma)} u_n
$$

この形からも、$ \gamma=2-\sqrt{2} $ で係数が揃い、同じ係数行列を使えることがわかります。

8. まとめ

TR‑BDF2 は台形(TR)とBDF2を 1ステップ内で 連続して解く二段法 です。$ \gamma = 2-\sqrt{2} $ を用いると $L$安定 で、同一行列の再利用により 強剛性かつ大規模 な問題で実装効率が高いです。
TCAD（Sentaurus/ATLAS など）では採用例がある一方、一般SPICE系（Spectre, Xyce, ngspice, LTspice 等）では trap / gear（＋modified‑trap）中心で、TR‑BDF2は一般的ではありません。
実務Tips：スイッチング多め・イベント密集ならまず trap / gear / modified‑trap と誤差制御を見直し、デバイスレベルの強剛性 では TR‑BDF2 を検討する、という棲み分けが現場で扱いやすそうです。

Q&A

Q. TR‑BDF2はGear(BDF2)の"改良版"ですか？
A. 別物です。台形→BDF2 の 二段合成 で 二次精度かつ$L$安定 を狙った方法（one‑step, two‑stage）です。BDF2単体は多段法（履歴2点が必要）で、開始やイベント処理の取り回しがやや異なります。
Q. SPICEでもTR‑BDF2を試せますか？
A. 私が調べた範囲ではメジャーツールの公式選択肢にはないです。ただし、MATLAB/Simulink の ode23tb（TR‑BDF2実装）をMNAへ適用すれば検証できそうです。

参考文献・リンク（抜粋）

R. C. Y. Chin, L. F. Shampine, "Analysis and Implementation of TR-BDF2"（Hosea–Shampine, 1996 などのTR‑BDF2解説論文の系譜）
MATLAB Documentation: ode23tb — TR-BDF2 solver（TR‑BDF2と同一iteration matrixの解説）
Hairer & Wanner, Solving Ordinary Differential Equations II（BDF/SDIRK/Rosenbrock の背景）
Synopsys Sentaurus Device User Guide / Training資料（過渡のデフォルトがTR‑BDF2である記述がある版があります）
Silvaco ATLAS User's Manual / Release Notes（TR‑BDF2の離散化に関する記述）
Cadence Spectre User Guide（過渡 method=trap/gear2）
Xyce User Guide（過渡はTrapezoid/Gear中心）
ngspice Manual / 開発者ML（過渡はTrapezoid/Gearの2種が基本）
LTspice ヘルプ／解説（Trap/Gear/Modified‑Trap）

実装メモ(自分用) 　

ニュートンのヤコビアンは常に $ \mathbf{J}=\mathbf{G}+\alpha\mathbf{C} $ の形にする。
$ \gamma=2-\sqrt{2} $ 固定で開始し、TR段のLU分解をBDF2段でも再利用する。
誤差推定は埋め込み（2次/3次差）で $E$ を算出、$h_{\text{new}}=\mathrm{safety}\cdot h\cdot E^{-1/3}$。
ブレークポイント境界を跨がないよう、必要なら $h$ を短縮してから2段を実行する。