はじめに

電気回路シミュレーションでは、ダイオードやトランジスタなどの非線形要素を含む回路を時間領域で解析する必要があります。
本記事では、オープンソースの数値ライブラリ SUNDIALS と、スパース行列専用のLUソルバ KLU を組み合わせて、ダイオードを含むような非線形微分代数方程式(DAE)をどのように解いていくかを解説します。

SUNDIALS: 常微分方程式(ODE)や微分代数方程式(DAE)、非線形方程式を解くためのライブラリ群。
- 例: IDA, CVODE, KINSOLなど。
KLU: SuiteSparseプロジェクトの一部で、電気回路シミュレーションに特化したスパースLU分解ソルバ。
- SUNDIALSと連携させることで、大規模かつスパースなヤコビ行列を効率よく分解可能。

背景: 非線形素子と微分代数方程式(DAE)

電気回路シミュレーションでは、修正ノード解析(MNA) を行うことで、ノード電圧や電流の連立方程式を作成します。
しかし、ダイオードやトランジスタなどの非線形素子が存在すると、その電流電圧関係は指数関数など線形では表せない形式をとります。結果として、以下のような 微分代数方程式 (DAE) を解く必要が出てきます。

\mathbf{f}\bigl(\mathbf{y}, \mathbf{y}'\bigr) = 0

SUNDIALSのIDAライブラリは、このような非線形DAEを数値的に解くための強力なツールで、回路方程式を時間領域で解く際に有効です。

SUNDIALS + KLUの仕組み

SUNDIALS（特にIDAやCVODEなど）は、非線形方程式をニュートン法で解く際に、ヤコビ行列を用いた線形システムを反復解法する仕組みを持っています。
行列が大規模かつスパースである場合、通常の密行列ソルバはメモリ・計算時間ともに非効率になりがちです。一方で、KLUのようなスパースLUソルバを併用することで、行列のスパース性を活かし、大規模問題でも高速に解けるようになります。

非線形残差 ($ \mathbf{r}(\mathbf{y}, \mathbf{y}', t) = 0 $) をユーザーが定義する
ヤコビ行列 ($ \frac{\partial \mathbf{r}}{\partial \mathbf{y}} $) を差分近似または解析的に計算
IDA が内部ニュートン法でステップごとにヤコビ行列を構築し、KLUでスパースLU分解
解が収束するまで連立線形方程式を解くループを回す

こうした流れにより、ダイオードのような非線形要素を含んだ回路でも効率的に数値解を求められます。

ダイオード特性の例

ダイオードのI-V特性を簡単に書くと：

I_D = I_S \bigl( e^{(V_D / (n V_T))} - 1 \bigr)

(I_S): 飽和電流
(V_D): ダイオードの電圧差
(n): ideality factor
(V_T): thermal voltage

ダイオードが接続されたノード電圧 (vD) に対する電流 (iD) が、指数関数的に変化します。これをMNA方程式やDAEに組み込むと、非線形方程式が登場します。

ヤコビ行列の組み立て

ダイオードなどの指数関数的な素子を含むと、ヤコビ行列の該当部分は指数関数の微分が含まれます。
解析的ヤコビであれば、指数関数の微分を直接書けますし、数値的な精度・収束性の面で有利です。一方、差分近似は実装が簡単ですが、($ \Delta v $) の設定や数値誤差で収束が遅くなる可能性があるため、大規模回路では解析的ヤコビが望まれます。

サンプルコード (IDA + KLU, C++20)

以下に、SUNDIALS(IDA) と KLU を併用した最小限のサンプルコードを示します。

ダイオード特性自体は簡略化しています。
C++20 の std::format を使用するため、GCC 12以降・Clang 13以降などを推奨。

// filename: example_sundials_klu_diode.cpp

#include <cmath>
#include <format>
#include <iostream>

// =========== SUNDIALS関連のヘッダ ==========
#include <idas/idas.h>                   // IDACreate, IDAInitなど
#include <nvector/nvector_serial.h>      // N_VNew_Serialなど
#include <sunmatrix/sunmatrix_sparse.h>  // SUNSparseMatrix
#include <sunlinsol/sunlinsol_klu.h>     // SUNLinSol_KLU
#include <idas/idas_ls.h>               // IDASetLinearSolver

// ダイオードパラメータ (例: 定数を仮定)
constexpr double IS = 1e-12;    // 飽和電流
constexpr double VT = 25.85e-3; // thermal voltage (約25.85mV, 室温)
constexpr double n  = 1.0;      // ideality factor

// ---------- 2次元のDAE(例) ------------
// y[0]: ダイオード電圧 (vD)
// y[1]: 他の状態量 (dummy)
// r_1 = vD' + I_D(vD) = 0
// r_2 = y[1] - sin(t) = 0
static int resfn(double t,
                 N_Vector yy,
                 N_Vector yp,
                 N_Vector rr,
                 void* user_data)
{
    double* y       = NV_DATA_S(yy);   // {vD, dummy}
    double* yprime  = NV_DATA_S(yp);   // {vD', dummy'}
    double* r       = NV_DATA_S(rr);   // {r1, r2}

    double vD = y[0];
    double iD = IS * (std::exp(vD/(n*VT)) - 1.0);  // ダイオード電流

    // r_1: vD' + iD = 0
    r[0] = yprime[0] + iD;

    // r_2: y[1] - sin(t) = 0
    r[1] = y[1] - std::sin(t);

    return 0; // 正常終了
}

int main()
{
    // ---- 変数数 N=2 ----
    int N = 2;

    // ---- IDAメモリ作成 ----
    void* ida_mem = IDACreate();
    if (!ida_mem) {
        std::cerr << "Error: IDACreate failed.\n";
        return -1;
    }

    // ---- 初期状態 (y, y') ----
    N_Vector y  = N_VNew_Serial(N);
    N_Vector yp = N_VNew_Serial(N);
    if (!y || !yp) {
        std::cerr << "Error: N_VNew_Serial failed.\n";
        return -1;
    }

    // y[0] = 0 (ダイオード電圧初期値)
    // y'[0] = 0
    // y[1] = 0
    // y'[1] = d/dt(sin(t))のt=0 => cos(0)=1 (仮定)
    NV_Ith_S(y,  0) = 0.0;
    NV_Ith_S(yp, 0) = 0.0;
    NV_Ith_S(y,  1) = 0.0;
    NV_Ith_S(yp, 1) = 1.0;

    // ---- IDA初期化 ----
    // (resfn)で定義した残差関数を登録
    int flag = IDAInit(ida_mem, resfn, 0.0, y, yp);
    if (flag != IDA_SUCCESS) {
        std::cerr << "Error: IDAInit failed.\n";
        return -1;
    }

    // ---- 誤差許容値 ----
    // 相対誤差1e-6, 絶対誤差1e-6に設定
    flag = IDASStolerances(ida_mem, 1e-6, 1e-6);
    if (flag != IDA_SUCCESS) {
        std::cerr << "Error: IDASStolerances failed.\n";
        return -1;
    }

    // ---- KLU用の行列 & ソルバ ----
    // N=2, 最大4要素のスパース行列(CSC形式)
    SUNMatrix A = SUNSparseMatrix(N, N, 4, CSC_MAT);
    if (!A) {
        std::cerr << "Error: SUNSparseMatrix creation failed.\n";
        return -1;
    }

    // 線形ソルバLS (KLU)
    SUNLinearSolver LS = SUNLinSol_KLU(y, A);
    if (!LS) {
        std::cerr << "Error: SUNLinSol_KLU creation failed.\n";
        return -1;
    }

    // IDAにスパースソルバLSをセット
    flag = IDASetLinearSolver(ida_mem, LS, A);
    if (flag != IDA_SUCCESS) {
        std::cerr << "Error: IDASetLinearSolver failed.\n";
        return -1;
    }

    // ---- 時間方向に解く (例: t=0 ~ 1.0) ----
    double t  = 0.0;
    double tout = 0.1;
    while (tout <= 1.0) {
        // IDASolveで時刻toutまで積分
        flag = IDASolve(ida_mem, tout, &t, y, yp, IDA_NORMAL);
        if (flag < 0) {
            std::cerr << "IDA error at t=" << t << "\n";
            break;
        }
        double vD  = NV_Ith_S(y, 0);
        double y2  = NV_Ith_S(y, 1);
        // C++20 std::formatでの出力
        std::cout << std::format("t={:.2f}, vD={:.6f}, y2={:.6f}\n",
                                 t, vD, y2);
        tout += 0.1;
    }

    // ---- 後処理 ----
    N_VDestroy(y);
    N_VDestroy(yp);
    SUNLinSolFree(LS);
    SUNMatDestroy(A);
    IDAFree(&ida_mem);

    return 0;
}

ビルド例 (Ubuntu)

g++ -std=c++20 example_sundials_klu_diode.cpp -o example_sundials_klu_diode \
    -I/usr/include \
    -L/usr/lib/x86_64-linux-gnu \
    -lsundials_ida -lsundials_nvecserial -lsundials_sunlinsollklu -lsundials_sunmatrixsparse \
    -lklu -lamd -lcolamd -lbtf -lsuitesparseconfig \
    -lm

フローチャート

サンプルコード解説

IDAメモリ作成:
IDACreate() でIDAのコンテキストを生成します。
初期状態設定:
N_VNew_Serial で未知数ベクトル（y）とその微分（yp）を用意し、初期値を設定。
残差関数resfn登録:
IDAInit で、微分代数方程式(DAE)を定義する関数 resfn をIDAに伝えます。
誤差許容値設定:
IDASStolerances で相対・絶対誤差を設定し、安定かつ精度を確保。
KLUソルバの生成・設定:
SUNSparseMatrix でスパース行列を確保し、SUNLinSol_KLU でKLUソルバを作成。
IDASetLinearSolver でIDAに登録します。
時間ループ:
IDASolve を呼び出し、少しずつ時間を進めながら（ここではtout += 0.1）非線形DAEを積分。各ステップで vD と y2 を出力。
後処理:
使ったメモリやソルバリソースを解放します。

数式表記とダイオードの残差

ダイオードの指数関数特性を含む非線形微分代数方程式は、例えば以下のように書くことができます。

\begin{cases}
r_1 = \frac{dv_D}{dt} + I_S\bigl(e^{v_D/(nV_T)} - 1\bigr) = 0 \\
r_2 = y_2 - \sin(t) = 0
\end{cases}

SUNDIALS(IDA) はこの $ \mathbf{r}(t, \mathbf{y}, \mathbf{y}') = 0 $ をニュートン法で反復し、ヤコビ行列を組み立てるときにKLUでスパースLU分解を行って解を求めます。
この仕組みにより、ダイオードのような非線形素子が含まれていても高速に解を収束させることが可能です。

よくある質問(FAQ)

Q1. 非線形要素が多い大規模回路でも動作しますか？

A. はい。SUNDIALS+KLUは、大規模かつスパースな系に最適化されています。非線形素子が増えてもヤコビ行列のスパース構造を活かせれば、KLUが効率的にLU分解を行い、計算時間とメモリを節約できます。

Q2. 解析的ヤコビと差分近似、どちらが良いですか？

A. 解析的ヤコビ（手動または自動微分ツールを用いる）が理想的です。差分近似だと、($ \Delta v $) の設定や数値誤差で収束が遅くなる場合があります。ただし、実装コストが高い場合は差分近似から始めるのも手です。

Q3. ダイオード以外にもMOSFETやBJTは扱えますか？

A. 同じ仕組みで扱えます。MOSFETやBJTなどより複雑なI-V特性でも、残差関数に適切にモデルを組み込めばニュートン法で解くことができます。

Q4. 収束しない場合の対策はありますか？

A. 初期値を工夫する、相対/絶対誤差許容値を緩める、タイムステップを小さくする、ヤコビ行列を解析的に計算する等が有効です。ダイオードなど急峻な特性の場合、初期ステップが適切でないと発散してしまうことがあります。

Q5. SUNDIALSをGPU上で動かせますか？

A. 2023年現在、KLU自体はCPU向け実装が中心で、SUNDIALSのKLU連携もCPU計算が前提です。GPU対応のスパースソルバやSUNDIALSのGPUサポートも一部検討されていますが、KLUはまだGPUをネイティブにはサポートしていません。

Q6. 実際のSpiceライクなシミュレータと同等の精度が得られますか？

A. SUNDIALS+KLUの組み合わせは、高い汎用性と安定性を持ちますが、商用Spice等では専用最適化や独自モデルを用いている場合があります。理論上は同等の精度を得られますが、大規模回路でのパフォーマンスや周辺機能の充実度は最適化次第です。

Q7. ダイオードのオン・オフ切り替えなど、イベント的な振る舞いはどのように扱いますか？

A. 実際には指数関数モデルで連続的にオン・オフが表現されるため、ハードスイッチング的なギャップは存在しません。もし「ある電圧を境に一気にモードが変わる」ような挙動を取り入れたいなら、SUNDIALSのイベント検知機能（root-findingなど）を利用して分岐タイミングを検出し、モデルを切り替える方法もあります。

Q8. SUNDIALSでの時間ステップ制御はどうすればいいですか？

A. IDAやCVODEは自動ステップ制御を行い、局所誤差を推定しながらステップサイズを調整します。IDASStolerances 等で相対・絶対誤差許容値を設定すれば、通常は十分安定に積分されます。さらに、最大/最小ステップや初期ステップ設定、ニュートン反復回数などを調整すれば、より細かな制御が可能です。

Q9. 行列がランク落ち（特異行列）する場合は？

A. 回路にフローティングノードがある、接地や参照がない、またはモデルが不適切なときに行列が特異化することがあります。KLUはランク欠損自体を検出しませんが、ニュートン法で「Singular matrix」エラーが返る場合があるので、回路構成をチェックしてください。

Q10. ダイオードやMOSFETが多数ある巨大回路で、KLUのメモリ使用量が大きくなる場合は？

A. 大規模回路ではスパース性を最大限活かしつつ、それでも要素数が膨大になることがあります。

順序付け (Ordering) でフィルインを削減（AMD, COLAMDなど）
klu_refactor でシンボリック解析の再利用
サブ回路に分割し、境界条件を受け渡す形で解く分割解法などを検討

Q11. SUNDIALSのIDAではなくCVODEでも同じことができますか？

A. 常微分方程式(ODE)寄りの形式ならCVODEでも可能です。ただ、電気回路のKCL/KVLをそのまま表現するとDAEになることが多いため、IDAの利用が自然となります。

Q12. SpiceのAC解析やDC解析のような機能はありますか？

A. DC解析（定常解）は、SUNDIALSのKINSOL（非線形方程式ソルバ）で可能です。AC解析（線形化解析）は周波数領域での手法なので、SUNDIALS単体では提供していません。時間領域シミュレーションからフーリエ変換等で周波数特性を推定する方法はありますが、Spiceが持つ小信号解析などとはアプローチが異なります。

まとめ

ダイオードのような非線形素子を含む回路でも、SUNDIALS+KLUの組み合わせで解ける
非線形DAE ($ \mathbf{f}(\mathbf{y}, \mathbf{y}')=0 $) をIDAがニュートン法で解き、スパースヤコビ行列はKLUで高速に分解
ダイオード以外にも、トランジスタやオペアンプなど多数の非線形素子に同様の方法を適用できる
数値的な安定性や初期値、差分ヤコビ vs 解析ヤコビの選択など考慮点はいくつかありますが、回路の大きさが増えるほどKLUの恩恵が大きいと言えます

本記事では最小限の例を取り上げましたが、より大規模・複雑な回路でも同じ仕組みで拡張が可能です。

以上