本文系统地介绍类SPICE集成电路仿真器的实现原理，包括改进节点分析（MNA）、非线性器件建模、DC/AC分析、时域/（复）频域仿真以及涉及的数值方法。
基于介绍的原理，实现了SPICE-like仿真器：https://github.com/cassuto/CSIM

1 理论基础

　　任何集总参数电路都能依据基尔霍夫电流定律（KCL）、基尔霍夫电压定律（KVL）和支路约束方程建立模型并通过解析法或数值法求解，进而实现计算机辅助电路分析（CACA，Computer Aided Circuit Analysis）。

　　CACA核心：

　　1. 建立电路数学模型：

• a）拓扑约束：利用图论分析电路，建立KCL、KVL方程。常见方法包括割集电压法、节点分析（NA, Nodal Analysis）法和C.W. Ho\(^{[1]}\)提出的改进节点分析（MNA, Modified Nodal Analysis）法等。UC Berkeley开发的SPICE即采用MNA方法\(^{[7]}\)。
• b）元件约束：根据元件的物理特性和分析的目标建立VCR约束。（复）频域分析可使用s域模型或相量（正弦稳态）模型建立VCR约束；时域分析则可用微分代数方程（DAE）建立VCR约束，或先使用s域模型求解最后进行Laplace逆变换。

　　2. 求解数学模型。主要有解析法和数值法两种。并非所有模型都容易找到解析解。常采用数值方法，对于线性方程组，可用LU分解法等；对于非线性方程，需通过迭代法将其线性化。对于微分方程，常使用数值积分。

　　3. 分析模型。主要进行误差和灵敏度分析。

　　CACA把上述过程总结为一套算法，让计算机自动完成。

2 术语约定

在不同参考资料中，相关术语的定义各有差别。本文统一采用如下规定：

• 网络（Network）：描述电路拓扑的无向图（带参考方向时为有向图），用\(G(V,E)\)表示；
• 网表（Netlist）：网络的一个实例；
• 端子（Pin）：元件的接线处；
• 端口（Port）：一对端子构成一个端口；
• 节点（Node）：连接两个或两个以上端子的交汇点；
• 支路（Branch）：对于一个二端元件，若两个端子分别接在节点\(s\)和\(t\)上，则该二端元件构成一条支路\((s,t) \in E\)；
• 关联（Associative）：给定节点\(n\)，如果存在支路\((s,t) \in E\)满足\(s=n\)或\(t=n\)，就称该支路与节点\(n\)关联；
• 回路（Circuit）：如果从节点\(n_0\)出发，能沿与之关联的支路\(b_0\)到达节点\(n_1\)，再以\(n_1\)为起点，重复上述过程并最终返回节点\(n_0\)，并且形成的访问序列\(<n_0,b_0,n_1,b_1,\cdots ,n_0>\)没有重复的节点，则该序列指出一个回路；

3 改进节点分析（MNA）前置内容

　　MNA是节点分析（NA）的增广，这里先介绍NA，如果您对此熟悉可以跳过本节。

3.1 NA方程的推导

　　我们从拓扑学角度推导NA方程。

　　将电路抽象为网络。首先为网络中各支路指定电流参考方向，使其成为有向网络，设其关联矩阵为：\(\mit{A}_{n \times b} = \begin{bmatrix} a_{jk} \end{bmatrix} = \begin{cases} 1, & \text{支路k参考电流从节点j流出}\\ -1, & \text{支路k参考电流从节点j流入} \\ 0, & \text{支路k与节点j无关联} \end{cases}\)

　　其中\(n\)为节点个数，\(b\)为支路个数。

　　设支路电压向量为\(\mit{U}=\begin{bmatrix} u_1, u_2, \cdots ,u_b \end{bmatrix} ^\mathrm{T}\)；支路电流向量为\(\mit{I}=\begin{bmatrix} i_1, i_2, \cdots ,i_b \end{bmatrix} ^\mathrm{T}\)。

　　根据KCL定律\(^{[2]}\)：

　　其中\(\mit{I}_s=\begin{bmatrix} i_{s_1}, i_{s_2}, \cdots ,i_{s_n} \end{bmatrix} ^\mathrm{T}\)为注入电流，\(i_{s_j}>0\)表示电流注入节点\(j\)，反之表示电流流出节点\(j\)。

　　设节点电压向量为\(\mit{U}_n=\begin{bmatrix} u_{n_1}, u_{n_2}, \cdots ,u_{n_n} \end{bmatrix} ^\mathrm{T}\)。在网络中任意选定一个节点作为参考节点，对于余下\((n-1)\)个节点，规定节点电压为该节点相对于参考节点的电位差。

　　根据KVL定律\(^{[2]}\)：

　　设支路导纳矩阵为：

　　通过导纳描述支路约束方程：

　　将式\((\ref{NA_KVL})\)代入式\((\ref{NA_BRAN})\)，得到：

　　再将式\((\ref{NA_KVL_BRAN})\)代入式\((\ref{NA_KCL})\)中，得到：

　　令\(\mit{Y} = \mit{A}\mit{G}\mit{A}^\mathrm{T}\)，最终得到节点分析方程：

　　其中\(\mit{Y}\)称为节点导纳矩阵。

　　给定网表可计算节点导纳矩阵\(\mit{Y}\)，也可由电流源支路确定\(\mit{I}_s\)（若无，置0）。最后只需解上述线性代数方程组，即可得出节点电压\(\mit{U_n}\)。

3.2 节点导纳矩阵的物理意义

　　上节已经推出节点导纳矩阵的表达式：

　　按照关联矩阵的定义，若支路\(k\)与节点\(i\)或\(j\)无关联，则\(a_{ik}a_{jk}=0\)。若支路\(k\)与节点\(i\)与\(j\)均有关联，则\(a_{ik}=\pm 1\)且\(a_{jk}=\mp 1\)，此时有\(a_{ik}a_{jk}=\begin{cases} -1 & (i \ne j) \\ 1 & (i = j) \end{cases}\)。

　　由此可知，\(\sum\limits_{k=1}^{b} g_k \cdot (a_{ik}a_{jk}) =y_{ij} \space (i \ne j)\)的物理意义为：若节点\(i\)与节点\(j\)之间有直接关联的支路（\(i \ne j\)），则\(y_{ij}\)就是两节点\(i\)和\(j\)之间关联的所有支路的导纳之和取负数，否则\(y_{ij}=0\)。把它定义为互导\(y_{ij} \space (i \ne j)\)。
　　\(\sum\limits_{k=1}^{b} g_k \cdot (a_{ik}a_{ik})=y_{ii}\)的物理意义为：\(y_{ii}=\sum\limits_{j=1 \\ j \neq i}^{n} y_{ij}\)，即所有与节点\(i\)直接关联的支路的导纳之和，把它定义为自导\(y_{ii}\)。

　　于是节点导纳矩阵就能简单表示为：

　　下面以一个实例来说明这种表示的好处：

　　按照自导和互导的物理意义，就能直接知道节点1的自导\(y_{11}\)为\(R_3\)和\(U_1\)的电导之和，节点1到2的互导\(y_{12}\)为\(U_1\)的电导取负数，其它同理，可直接得出节点导纳矩阵如下：

　　需要注意\(y_{1,3}=y_{3,1}=0\)，\(y_{2,4}=y_{4,2}=0\)，因为节点1和3、2和4之间没有直接关联的支路。另外，因为U1是理想电压源，\(G_{U_1}=\infty\)，因此这个电路用NA无法直接求解，这个例子揭示了NA的局限性。

4 改进节点分析（MNA）的原理

　　节点分析（NA）直接处理含无伴电压源（内阻为0）的支路时会遇到困难，因为这些支路的导纳为无穷大，方程无法求数值解。上文图1中的\(U_1\)就属于这种情况。

　　MNA解决NA局限性的思路很简单：对NA方程组进行增广。NA只分析节点电压，而MNA能同时分析支路电流，将两种状态变量混合在一起求解。

　　MNA混合方程如下：

　　其中\(\mit{Y}\)是节点导纳矩阵；\(\mit{U}\)是节点电压向量；\(\mit{I}\)是节点电流向量，对应方程\(\mit{Y}\mit{U}=\mit{I}\)与NA无异。此外，\(\mit{J}\)是支路电流向量，\(\mit{B}\)、\(\mit{C}\)、\(\mit{D}\)、\(\mit{E}\)是新增加的方程的系数矩阵。

　　这些系数矩阵用“橡皮图章”法（Rubber Stamps）机械地生成\(^{[1]}\)，实现电路分析的程序化。

4.1 方程生成算法

　　C.W. Ho提出的元件橡皮图章算法（Element Rubber Stamps）\(^{[1]}\)可以根据网表直接生成各子矩阵的值。算法初始时\(\mit{Y}=\mit{B}=\mit{C}=\mit{D}=\mit{I}=\mit{E}=0\)。遍历网表中所有元件，每遇到一个元件\(e\)，就将该类型元件对应的贡献值\(c(e)\)加到相应的子矩阵，遍历结束时各子矩阵就有了正确的值，从而直接产生方程组。因为每种元件的贡献值是常量，可以将贡献值编制成表格，也称为“表格法”。

　　在稍后的推导中可以看到，算法“将贡献值加到对应子矩阵”的的理论依据是线性电路的叠加性。MNA本来只能处理线性电路，但是稍后可以看到如何利用线性化处理非线性电路。

　　方程生成算法可以描述为：

　　算法的关键是贡献值\(c_Y(e)\)、\(c_B(e)\)、\(c_C(e)\)、\(c_D(e)\)、\(c_I(e)\)、\(c_E(e)\)。

　　下面推导了一些常见一端口和二端口线性元件的贡献值。

4.2 常见一端口线性元件对Y、B、C、D、I、E的贡献

　　假设一端口元件所在支路为\(k\)，电压与电流取关联参考方向\(s \rightarrow t\)。

4.2.1 阻抗元件

　　设元件导纳为\(y\)，支路约束方程为：

　　其中\(i_k\)为元件所在支路电流。

　　整理上式得：

　　对应MNA矩阵形式为：

　　如果多个阻抗元件并入支路，则支路的导纳就是它们之和。因此每当一个阻抗元件并入支路\(s \rightarrow t\)时，只需在MNA方程的分块矩阵Y中加上导纳：

即加上贡献值：

4.2.2 独立电压源（VS）

　　设VS的电压值设定为\(e_s\)。支路约束方程为：

　　其中\(i_k\)为VS所在支路的电流。

　　对应MNA矩阵形式为：

　　如果多个电压源串在同一支路，则支路电压为它们之和，因此每当一个VS串入支路\(s \rightarrow t\)时，只需加上贡献值：

4.2.3 独立电流源（CS）

　　设CS电流值设定为\(i_k\)，支路约束方程为：

　　如果多个电流源并入支路，则支路总电流就是它们之和，因此每当一个CS并入支路\(s \rightarrow t\)时，只需加上贡献值：

　　上述这些例子其实反映了线性电路的叠加性。这就是为什么任何元件加入电路时都只需在MNA对应矩阵加上贡献值。

4.3 常见二端口线性元件对Y、B、C、D、I、E的贡献

　　假设元件的端口1所在支路为\(k\)，电压与电流取关联参考方向\(s \rightarrow t\)；端口2所在支路为\(c\)，电压与电流取关联参考方向\(p \rightarrow q\)。

4.3.1 电压控制电压源（VCVS）

　　设VCVS电压放大倍数为\(\mu\)，控制电压为\((u_p - u_q)\)。支路约束方程为：

　　其中\(i_k\)为VCVS受控端口所在支路的电流。

　　对应MNA矩阵形式为：

　　可知每当一个VCVS加入支路\(s \rightarrow t\)、\(p \rightarrow q\)时，贡献为：

4.3.2 电压控制电流源（VCCS）

　 设VCCS的转移电导为\(g\)，控制电压为\((u_p - u_q)\)。支路约束方程为：

　　对应MNA矩阵形式为：

　　可知每当一个VCCS加入支路\(s \rightarrow t\)、\(p \rightarrow q\)时，贡献为：

4.3.3 电流控制电压源（CCVS）

　　设CCVS的转移电阻为\(r\)，控制电流为\(i_c\)，支路约束方程为：

　　其中\(i_k\)为CCVS受控端口所在支路电流。

　　对应MNA矩阵形式为：

　　可知每当一个CCVS加入支路\(s \rightarrow t\)、\(p \rightarrow q\)时，贡献为：

4.3.4 电流控制电流源（CCCS）

　　设CCCS的电流放大倍数为\(\alpha\)，控制电流为\(i_c\)。支路约束方程为：

　　对应矩阵形式为：

　　可知每当一个CCCS加入支路\(s \rightarrow t\)、\(p \rightarrow q\)时，贡献为：

5 稀疏矩阵计算

　　对于线性电路，按照上述方法可以建立MNA方程：

　　直接采用高斯列主元素消元法、LU分解法、雅各比迭代法等都能求解上述方程，具体参考教科书[3]。

　　但是，求解稠密矩阵方程需要\(O(n^3)\)的时间。如果观察到电路对应的系数矩阵\(\mit{A}\)是稀疏矩阵，就可以使用更优化的算法，因为而稀疏矩阵中存在大量零元素，利用稀疏矩阵算法，存储和计算时都可以跳过大量零元素，从而使算法所需的时间和空间大幅减少。

6 非线性元件的分析

　　MNA可以方便地分析线性电路，但无法直接处理非线性电路。

　　幸运的是，如果利用迭代法将非线性元件线性化，使每一步迭代都能用等效的线性元件替代，就能利用MNA分析和求解非线性电路。SPICE就采用这种方法\(^{[7]}\)。

6.1 牛顿-拉夫逊法求解非线性MNA方程

　　牛顿-拉夫逊法（Newton-Ralfsnn's method）是最常用求解非线性方程近似根的算法。

　　牛顿-拉夫逊法通过如下迭代格式计算非线性方程\(f(\mit{x}) = \mit{0}\)的近似根\(\mit{x}\)：

　　根据MNA方程\((\ref{MNA})\)，设

　　根据式\((\ref{equ:newtown_fx})\)：

　　这其中涉及的雅克比矩阵有：

　　利用雅可比矩阵将式\((\ref{equ:newtown_fx})\)改写为：

　　再将式\((\ref{equ:fMNA})\)代入，整理可得迭代格式：

　　即\(\mit{A}'(\mit{x}^{(k)}) \cdot \mit{x}^{(k+1)}=\mit{z}'(\mit{x}^{(k)})\)，可见迭代格式与线性代数方程组形式保持一致。给定初值\(\mit{x}^{(0)}\)，计算\(\mit{A}'(\mit{x}^{(0)})\)、\(\mit{z}'(\mit{x}^{(0)})\)，然后求解线性代数方程组，得到\(\mit{x}^{(1)}\)，再将\(\mit{x}^{(1)}\)作为新的初值，重复上述过程，生成迭代序列\(\{\mit{x}^{(k)}\}\)，直到\(\|\mit{x}^{(k+1)} - \mit{x}^{(k)}\|\)小于设定的误差限时，可认为迭代收敛，取近似根\(\tilde{\mit{x}} = \mit{x}^{(k+1)}\)。

　　牛顿-拉夫逊迭代的本质是将非线性问题分成若干线性问题，这就启发我们用该方法实现元件的线性化，从而允许用MNA分析非线性元件。

6.2 理想PN结模型

　　理论上任何非线性元件都可以用动态电压源或电流源等效代替。为了便于应用MNA，这里使用动态电流源代替。

　　设理想PN结位于支路\(s \rightarrow t\)上，理想PN结电流的近似方程为：

　　其中\(I_0\)为反向饱和电流，\(U_T\)为温度电压当量，这两个参数都是由物理特性决定的。

　　将PN结作为动态电流源注入支路\(s \rightarrow t\)，根据独立电流源支路的结论\((\ref{equ:cs_contrib})\)，有：

　　代入牛顿-拉夫逊迭格式\((\ref{equ:mna_newtown_format})\)：

　　对比之前推出的阻抗元件支路的结论（式\(\ref{equ:mna_y_comp}\)），可以认为\((\ref{equ:pn_newtown_left})\)描述的是等效电阻，其电导\({g_d}^{(k)}\)随迭代次数\(k\)动态变化，然而在本轮迭代内是不变的，可视作线性元件：

　　再考虑式\((\ref{equ:mna_newtown_format})\)，右边式子可展开为：

　　其中：

　　从物理意义上看，式\((\ref{equ:PN_i_d})\)描述的是动态电流源\(i_d\)与动态电阻\(g_{eq}\)并联，如图2所示，这样就实现了元件的线性化，每次迭代都可以用MNA分析了。

　　至此式\((\ref{equ:mna_newtown_format})\)左右两边都已确定，得到MNA方程组：

　　每当一个理想PN结加入支路\(s \rightarrow t\)时，只需对子矩阵作如下更新：

　　值得注意的是整个求解过程是迭代进行的，每轮迭代都要重新计算等效电路的参数并重新求解MNA，即：

　　给定初值\(\mit{x}^{(0)}\)（这其中包含PN结两端的节点电压\({u_s}^{(0)}\)、\({u_t}^{(0)}\)），代入式\((\ref{equ:pn_newtown_left})\)和式\((\ref{equ:pn_newtown_right})\)可得线性代数方程组\((\ref{equ:new_rn_mna})\)，解线性代数方程组可以得到\(\mit{x}^{(1)}\)，再将\(\mit{x}^{(1)}\)作为新的初值，如此迭代。当相邻两次迭代获得的解\(\mit{x}^{(k+1)}\)与\(\mit{x}^{(k)}\)满足\(\| \mit{x}^{(k+1)}-\mit{x}^{(k)} \| \lt \epsilon\)时，就可认为迭代收敛，可以取近似解\(\mit{x}^{(k+1)}\)。

　　下面给出了非线性电路分析的算法流程。

6.2.1 收敛性问题

　　在PN结特性方程的牛顿-拉夫逊迭代中，存在如下图所示的异常情况：

　　图中第\(k+1\)步迭代时，由于指数函数的迅猛增长，\(i_d^{(k+1)}\)超出机器数所能表示的范围，产生上溢，使得迭代无法进行下去。

　　考虑到实际电路中不可能出现如此大的电流（双精度浮点数最大值约为\(10^{308}\)）；另外在结电流方程中，当\(y\)急剧增长时，\(x\)的变化范围却很小。因此可以将PN结电压限制在较小范围内，以避免数值溢出\(^{[6]}\)。

　　PN结临界电压是V-I曲线中曲率半径最小的点，当PN结电压大于临界电压时，结电流开始急剧增加，因此可用PN结临界电压\(U_{th}\)作为阈值的参考。

　　结电压限制算法一般基于经验来构造。一种最简单的阈值限制算法是\(U'_d = max(U_d, 10U_{th})\)。

　　这种算法将结电压限制在\(10U_{th}\)以下，但限制后的V-I曲线在\(U_d=10U_{th}\)处的导数不存在，大于\(10U_{th}\)后导数为0，由此造成不收敛。

　　SPICE中的限制算法\(^{[7]}\)要复杂得多，当\(U_d > U_{th}\)时，采用以电流\(I_d\)为变量的迭代；当\(U_d \le U_{th}\)时，采用正常的迭代格式。

6.3 收敛条件

　　对于MNA方程中的节点电压\(U\)和支路电流\(J\)，SPICE采用独立的收敛条件。设\(RELTOL\)为相对误差限，\(ABSTOL\)为绝对误差限。当：

　　并且，当\(|J^{(k+1)} _b- J^{(k)}_b| \le RELTOL \cdot J_{b,max} + ABSTOL\)时，认为迭代收敛。

　　其中\(U_{n,max}=max(|U^{(k+1)}_n|, |U^{(k)}_n|)\)；\(U^{(k+1)}, U^{(k)}\)为相邻两次迭代的结果。\(J_b\)同理。

7 直流扫描分析（DC Sweep）

　　直流分析用于求静态工作点，而直流扫描分析遍历参数，输出各参数下电路的静态工作点。

7.1 特殊情形

　　直流分析反映的是输入为直流（即频率\(\omega=0\)）时的状态，需要特殊处理动态元件。

　　电容在直流下的容抗为$ \lim \limits_{\omega \rightarrow 0} \frac{1}{j \omega C} = \infty $，显然直流分析时电容应视为两端开路。

　　电感在直流下的感抗为$ \lim \limits_{\omega \rightarrow 0} j \omega L = 0 $，显然直流分析时应视为两端短路。

　　此外所有作为信号源的电压源视为短路、作为信号源的电流源视为开路。

7.2 直流分析的过程

　　设目标参数\(p\)，扫描范围\([p_{min}, p_{max}]\)，扫描步长\(s\)。线性扫描共需要\(\frac{p_{max}-p_{min}}{s}\)次方程求解，每次将目标参数设定为\(p(n)=p_{min}+ns\)，通过改进节点分析（MNA）建立的方程解出对应的节点电压\(U(n)\)。这样\(U(n)\)就形成了直流扫描分析的结果。

　　实用中，有时需要使用对数步进来扫描。

8 交流扫描分析（AC Sweep）

　　AC Sweep分析是正弦稳态电路在频域上的小信号分析。输入变量是正弦频率\(\omega\)，输出变量是电路中各节点电压的频率响应（幅度和相位）。

　　交流分析采用相量法，电压电流都采用相量表示，仍然利用MNA求解，只不过MNA中各矩阵都定义在复数域上。

8.1 电容的相量模型

　　理想电容在正弦稳态电路中的VCR表示为

　　每当一个理想电容加入支路\(s \rightarrow t\)时，只需对MNA的子矩阵的值作如下更新：

8.2 电感的相量模型

　　类似地，理想电感在正弦稳态电路中的VCR表示为

　　每当一个电感加入支路\(s \rightarrow t\)时，只需对MNA的子矩阵的值作如下更新：

8.3 交流分析的过程

　　交流分析非常重要的假设是小信号。在小信号模型中，非线性元件可以在静态工作点处线性化，例如PN结可通过动态电阻\(g_d(u)\)等效。因此在进行交流分析之前，先进行直流分析，确定电路静态工作点。静态工作点确定，式\((\ref{equ:mna_newtown_format})\)中所有雅克比矩阵的值也都确定。这样在交流分析时，不必迭代，而是将其视作线性方程来处理。

9 复频域分析（s域）

　　对电路的微分方程进行Laplace变换可以得到s域上的代数方程，这些代数方程可以用与上一节AC分析相同的方法建立和求解。事实上，上节所述的相量模型可以看作\(s = j\omega\)的特殊情况，这也反映了频域和复频域的关系。

　　s域的MNA混合方程变为：

9.1 Laplace变换

　　给定时域激励信号f(t)，可通过Laplace变换得到复频域上的激励\(F(s)=\mathscr{L}[f(t)]\)，将F(s)代入激励源模型中，求解MNA方程即可得到节点电压和分支电流的频域响应\(\begin{bmatrix} \mit{U(s)} & \mit{J(s)} \end{bmatrix}^{T}\)。

9.2 Laplace逆变换

　　对于上面得到的频域响应，可通过逆变换得到时域响应\(\begin{bmatrix}\mathscr{L}^{-1}[\mit{U(s)}] & \mathscr{L}^{-1}[\mit{J(s)}]\end{bmatrix}^{T}\)。

10 瞬态分析（时域分析）

　　上节已经给出了从复频域变换到时域的分析方法，下面介绍直接在时域进行分析的方法。

　　时域上理想电容和电感的VCR为常微分方程：

　　计算机无法直接处理连续时间系统。可以将时间离散化，然后利用数值方法求解。

10.1 线性多步：隐式GEAR法

　　对于电容或电感特性方程中所出现的形如\(f(x,t) = \frac{dx}{dt}\)的常微分方程，可通过积分\(x = \int f(x,t) dt\)来求解。

　　根据黎曼积分的定义，连续时间域上的积分\(x = \int f(x,t) dt\)可通过离散时间域上的累积来近似：

　　其中\(h_n=t_{n+1}-t_n\)称为第\(n\)步的步长。

　　将式\((\ref{equ:num_int})\)写成迭代格式：

　　这是一阶显式单步法。单步法的收敛性可参考资料\(^{[3]90-92}\)。为了获得更精确的解，这里采用线性多步法。线性多步法是前\(p\)步解的线性组合，单步法正是线性多步法的特例。

　　其中\(p\)是步数。\(a_i\)、\(b_i\)是常系数。

　　利用泰勒展开构造线性多步法\(^{[5]}\)，\(a_i\)、\(b_i\)应满足：

　　选取一些特殊的\(p\)、\(a_i\)、\(b_i\)值，就能构造出不同的迭代方法。当\(b_{-1} \ne 0\)时为隐式方法，当\(b_{-1} = 0\)时为显式方法。对于隐式GEAR：

　　给定阶数\(k=1,2,\cdots\)，联立式\((\ref{equ:gear_pb})\)与式\((\ref{equ:lin_multi_step_ab})\)可以解出系数\(a_i\)、\(b_i\)。从结果来看，一阶隐式GEAR就是隐式欧拉法（Implicit Euler's method）。4阶隐式GEAR迭代格式如下：

10.2 线性多步法中的迭代

　　在线性多步法迭代格式\((\ref{equ:multi_step_format})\)中，式子左侧为待求的量\(x_{n+1}\)，而式子右侧也依赖于待求量\(x_{n+1}\)，因此待求量无法直接计算。解决办法是解方程，设方程\(f(x_{n+1}) = x_{n+1} - \sum_{i=0}^{p}a_i x_{n-i} + h_n\sum_{i=-1}^{p}b_i f(x_{n-i},t_{n-i}) = 0\)，则只需通过迭代解出该方程的根\(x_{n+1}\)即可。迭代格式为：

　　迭代到\(\|x_{n+1}^{m+1}-x_{n+1}^{m}\|\)小于给定误差限时，可以取\(x_{n+1}=x_{n+1}^{m+1}\)。

10.3 预报-校正法

　　从上节可以看出，每步线性迭代中又包含若干\(x_{n+1}^{m+1}=g(x_{n+1}^{m})\)这样的迭代。初值\(x_{n+1}^0\)的选取将直接影响到迭代次数，因此初值的选取十分重要。相比于随机给定一个初值，通过预报器预测的初值可能更接近真实值，再进行线性多步迭代时，所需迭代次数将减少。

　　这里采用显式欧拉法实现预报器，预报值\(x_{n+1}^{0}\)计算为

10.4 自适应步长控制算法

　　瞬态分析中，如果时间步长\(h_n\)选取过大会造成局部截断误差偏大，甚至得出完全错误的结果；而如果步长选取太小则会使得计算量增加。手工选取步长非常繁琐，而且固定步长在某些区间内往往不是最优，因此一般采用变步长算法。

　　\(h_n\)的选取应使得第\(n\)步的局部截断误差\(\xi_{L}^{(n+1)} = |x(t_{n+1}) - x_{n+1}|\)满足：

　　其中\(\xi\)是设定的绝对误差限；\(\xi_r\)是设定的相对误差限。一般来说局部截断误差无法精确计算，只能通过Milne公式估计。

　　自适应步长控制中，先设定一个足够小的初始步长\(h_0\)，每进行一次迭代，就计算出本次迭代的局部截断误差\(\xi_{L}^{(n+1)}\)，再通过式\((\ref{equ:adpt_step})\)判定步长的好坏：

• 若\(q \gt 1\)，则说明局部截断误差大于设定的误差限，步长偏大；
• 若\(q \lt 1\)，则说明局部截断误差已经小于设定误差限，若q远小于1则说明步长过小。

　　根据\(q\)对步长进行动态调整：

　　其中\(k\)是线性多步法的阶数。

　　假设当前仿真时间为\(t_{n+1}\)，首先利用线性多步法求出解\(x_{n+1}\)，再估计局部截断误差\(\xi_{L}^{(n+1)}\)，按上式计算新步长\(h'_{n+1}\)，若\(h'_{n+1}<h_{n+1}\)，则说明局部截断误差过大，此时仿真时间不向前推进，而是按新步长重新计算当前时间的解\(x_{n+1}\)，重复上述过程；反之则说明局部截断误差已经满足要求，仿真时间可以推进到\(t_{n+2}\)，令\(h_{n+2}=h'_{n+1}\)，结束本次步长调整。

　　下图为瞬态仿真实例：

　　图中显示了步长的自适应调整，时间轴是均匀的，当信号变化接近线性时数据点比较稀疏，反之则比较稠密。

　　算法能保证每步迭代局部截断误差的估计值不超出规定的误差限，然而，每次调整步长都要重新计算线性方程组，对于信号快速变化的电路（例如开关电路），步长可能会频繁振荡，从而使仿真器将大量时间花费在寻找合适的步长上。因此，数字部分需要不同于模拟部分的专用算法来仿真。

10.5 常见储能元件的时域模型

10.5.1 电容的时域模型

　　电容的时域特性描述为：

　　利用线性多步法\((\ref{equ:multi_step_format})\)得到迭代格式：

　　将上式展开并移项（\(b_{-1} \ne 0\)），可得

　　其中：

　　\(g_{eq}\)从物理上可解释为等效电导，\(i_{eq}\)从物理上可解释为等效电流源，于是得到了电容的等效线性化模型，如图4所示。

　　根据之前推出的阻抗元件支路对MNA各子矩阵的贡献（式\(\ref{equ:admit_elem_contrib}\)）和独立电流源支路对MNA各子矩阵贡献值（式\(\ref{equ:cs_contrib}\)）可知对应MNA方程为：

　　因此可知每当一个电容加入支路\(s \rightarrow t\)时，只需对MNA各子矩阵作如下更新：

　　同时应在程序中应设置标记，指出参数\({g_{eq}}^{(n)}\)和\({i_{eq}}^{(n)}\)是需要迭代计算的。给定步长\(h_n\)，初值\(u_0={u_s}^{(0)}-{u_t}^{(0)}\)，根据上式生成MNA方程，可解出节点电压\({u_s}^{(1)}\)、\({u_t}^{(1)}\)，以此类推。

10.5.2 电感的时域模型

　　电感的时域特性描述为：

　　利用线性多步法\((\ref{equ:multi_step_format})\)得到迭代格式：

　　整理并移项（\(b_{-1} \ne 0\)），可得

　　其中：

　　从物理意义上看，电感可等效为动态电阻\(r_{eq}\)与独立电压源\(u_{eq}\)串联，如图5所示。

　　根据之前推出的独立电压源的贡献值（式\(\ref{equ:vs_contrib}\)）和电流控制电压源支路对MNA各子矩阵贡献值（式\(\ref{equ:ccvs_contrib}\)）可知对应MNA方程为：

　　因此每当一个电感加入支路\(s \rightarrow t\)时，只需对MNA各子矩阵作如下更新：

　　同时应在程序中应设置标记，指出参数\({r_{eq}}^{(n)}\)和\({u_{eq}}^{(n)}\)是需要迭代计算的。

10.6 时域分析总流程

• 1 初始化：建立MNA方程。设置当前时间\(t=0\)，设置积分步\(s\)、阶数\(ord\)为初值
• 2 用\(s\)、\(ord\)计算GEAR系数
• 3 解MNA方程（流程图见6.2）
• 4 更新时间\(t'=t+s\)
• 5 动态调整积分步\(s\)、阶数\(ord\)
• 6 判断终止条件，不终止则循环执行 2

11 程序实现

　　详见文章开头给出的github链接。

参考资料

[1] C.W. Ho; Ruehli, A.; Brennan, P. The modified nodal approach to network analysis[J]. IEEE, doi:10.1109/tcs.1975.1084079, 1975: 0–509.

[2] 邱关源. 电路[M]. 第5版, 高等教育出版社, 2006: 391-392.

[3] 李建良等. 计算机数值方法[M]. 东南大学出版社, 2009.

[5] Timothy Sauer. Numerical Analysis[M]. 2nd Edition, George Masonry University, 2011: 336-339.

[6] Thomas L. Quarles. Analysis of performance and convergence issues for circuit simulation[R], University of California, Berkeley Technical Report No. UCB/ERL M89/42, 1989: 30-31.

[7] L. W. Nagel. SPICE2: A Computer Program to Simulate Semiconductor Circuits[R]. University of California, Berkeley Technical Report No. UCB/ERL M520, 1975: 138-142