集成电路仿真器（SPICE）的实现原理 - TheWind

　　本文系统地介绍类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)∈E；
• 关联（Associative）：给定节点n，如果存在支路(s,t)∈E满足s=n或t=n，就称该支路与节点n关联；
• 回路（Circuit）：如果从节点n0出发，能沿与之关联的支路b0到达节点n1，再以n1为起点，重复上述过程并最终返回节点n0，并且形成的访问序列<n0,b0,n1,b1,⋯,n0>没有重复的节点，则该序列指出一个回路；

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

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

3.1 NA方程的推导#

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

　　将电路抽象为网络。首先为网络中各支路指定电流参考方向，使其成为有向网络，设其关联矩阵为：An×b=[ajk]={1,支路k参考电流从节点j流出−1,支路k参考电流从节点j流入0,支路k与节点j无关联

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

　　设支路电压向量为U=[u1,u2,⋯,ub]T；支路电流向量为I=[i1,i2,⋯,ib]T。

　　根据KCL定律[2]：

　　其中Is=[is1,is2,⋯,isn]T为注入电流，isj>0表示电流注入节点j，反之表示电流流出节点j。

　　设节点电压向量为Un=[un1,un2,⋯,unn]T。在网络中任意选定一个节点作为参考节点，对于余下(n−1)个节点，规定节点电压为该节点相对于参考节点的电位差。

　　根据KVL定律[2]：

　　设支路导纳矩阵为：

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

　　将式(2)代入式(3)，得到：

　　再将式(4)代入式(1)中，得到：

　　令Y=AGAT，最终得到节点分析方程：

　　其中Y称为节点导纳矩阵。

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

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

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

　　按照关联矩阵的定义，若支路k与节点i或j无关联，则aikajk=0。若支路k与节点i与j均有关联，则aik=±1且ajk=∓1，此时有aikajk={−1(i≠j)1(i=j)。

　　由此可知，b∑k=1gk⋅(aikajk)=yij (i≠j)的物理意义为：若节点i与节点j之间有直接关联的支路（i≠j），则yij就是两节点i和j之间关联的所有支路的导纳之和取负数，否则yij=0。把它定义为互导yij (i≠j)。
　　b∑k=1gk⋅(aikaik)=yii的物理意义为：yii=n∑j=1j≠iyij，即所有与节点i直接关联的支路的导纳之和，把它定义为自导yii。

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

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

　　按照自导和互导的物理意义，就能直接知道节点1的自导y11为R3和U1的电导之和，节点1到2的互导y12为U1的电导取负数，其它同理，可直接得出节点导纳矩阵如下：

　　需要注意y1,3=y3,1=0，y2,4=y4,2=0，因为节点1和3、2和4之间没有直接关联的支路。另外，因为U1是理想电压源，GU1=∞，因此这个电路用NA无法直接求解，这个例子揭示了NA的局限性。

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

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

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

　　MNA混合方程如下：

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

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

4.1 方程生成算法#

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

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

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

　　算法的关键是贡献值cY(e)、cB(e)、cC(e)、cD(e)、cI(e)、cE(e)。

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

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

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

4.2.1 阻抗元件#

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

　　其中ik为元件所在支路电流。

　　整理上式得：

　　对应MNA矩阵形式为：

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

即加上贡献值：

4.2.2 独立电压源（VS）#

　　设VS的电压值设定为es。支路约束方程为：

　　其中ik为VS所在支路的电流。

　　对应MNA矩阵形式为：

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

4.2.3 独立电流源（CS）#

　　设CS电流值设定为ik，支路约束方程为：

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

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

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

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

4.3.1 电压控制电压源（VCVS）#

　　设VCVS电压放大倍数为μ，控制电压为(up−uq)。支路约束方程为：

　　其中ik为VCVS受控端口所在支路的电流。

　　对应MNA矩阵形式为：

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

4.3.2 电压控制电流源（VCCS）#

　 设VCCS的转移电导为g，控制电压为(up−uq)。支路约束方程为：

　　对应MNA矩阵形式为：

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

4.3.3 电流控制电压源（CCVS）#

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

　　其中ik为CCVS受控端口所在支路电流。

　　对应MNA矩阵形式为：

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

4.3.4 电流控制电流源（CCCS）#

　　设CCCS的电流放大倍数为α，控制电流为ic。支路约束方程为：

　　对应矩阵形式为：

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

5 稀疏矩阵计算#

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

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

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

6 非线性元件的分析#

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

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

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

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

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

　　根据MNA方程(11)，设

　　根据式(12)：

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

　　利用雅可比矩阵将式(12)改写为：

　　再将式(13)代入，整理可得迭代格式：

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

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

6.2 理想PN结模型#

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

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

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

　　将PN结作为动态电流源注入支路s→t，根据独立电流源支路的结论(9)，有：

　　代入牛顿-拉夫逊迭格式(15)：

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

　　再考虑式(15)，右边式子可展开为：

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

　　至此式(15)左右两边都已确定，得到MNA方程组：

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

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

　　给定初值x(0)（这其中包含PN结两端的节点电压us(0)、ut(0)），代入式(16)和式(17)可得线性代数方程组(19)，解线性代数方程组可以得到x(1)，再将x(1)作为新的初值，如此迭代。当相邻两次迭代获得的解x(k+1)与x(k)满足‖x(k+1)−x(k)‖<ϵ时，就可认为迭代收敛，可以取近似解x(k+1)。

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

6.2.1 收敛性问题#

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

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

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

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

　　一种最简单的阈值限制算法是x′=max(x,10Uth)，将结电压限制在10Uth以下，但限制后的V-I曲线在Ud=10Uth处的导数不存在，大于10Uth后导数为0，造成不收敛。

　　SPICE中的限制算法[7]更合理，当Ud>Uth时，采用以电流Id为变量的迭代（解决反函数）；当Ud≤Uth时，采用正常的迭代格式。

6.3 收敛条件#

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

　　并且，当|J(k+1)b−J(k)b|≤ξr⋅Jb,max+ξ时，认为迭代收敛。

　　其中Un,max=max(|U(k+1)n|,|U(k)n|)；U(k+1)n,U(k)n为相邻两次迭代的结果。Jb同理。

7 直流扫描分析（DC Sweep）#

　　至此，我们搭建的框架可以实现SPICE的第一个应用——直流扫描分析，即遍历参数，输出各参数下电路的静态工作点。

7.1 特殊情形#

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

　　电容在直流下的容抗为limω→01jωC=∞，显然直流分析时电容应视为两端开路。

　　电感在直流下的感抗为limω→0jωL=0，显然直流分析时应视为两端短路。

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

7.2 直流分析的过程#

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

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

8 交流扫描分析（AC Sweep）#

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

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

8.1 电容的相量模型#

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

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

8.2 电感的相量模型#

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

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

8.3 交流分析的过程#

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

9 复频域分析（s域）#

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

　　s域的MNA混合方程变为：

9.1 Laplace变换#

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

9.2 Laplace逆变换#

　　对于上面得到的频域响应，可通过逆变换得到时域响应[L−1[U(s)]L−1[J(s)]]T。

10 瞬态分析（时域分析）#

　　上节给出了从复频域变换到时域的分析方法，下面直接在时域进行分析，这也是SPICE采用的方法。

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

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

10.1 线性多步：隐式GEAR法#

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

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

　　其中hn=tn+1−tn称为第n步的步长。

　　将式(20)写成迭代格式：

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

　　其中p是步数。ai、bi是常系数。

　　利用泰勒展开构造线性多步法[5]，ai、bi应满足：

　　选取一些特殊的p、ai、bi值，就能构造出不同的迭代方法。当b−1≠0时为隐式方法，当b−1=0时为显式方法。对于隐式GEAR：

　　给定阶数k=1,2,⋯，联立式(23)与式(22)可以解出系数ai、bi。从结果来看，一阶隐式GEAR就是隐式欧拉法（Implicit Euler's method）。4阶隐式GEAR迭代格式如下：

10.2 线性多步法中的迭代#

　　在线性多步法迭代格式(21)中，式子左侧为待求的量xn+1，而式子右侧也依赖于待求量xn+1，因此待求量无法直接计算。解决办法是解方程，设方程f(xn+1)=xn+1−∑pi=0aixn−i+hn∑pi=−1bif(xn−i,tn−i)=0，则只需通过迭代解出该方程的根xn+1即可。迭代格式为：

　　迭代到‖xm+1n+1−xmn+1‖小于给定误差限时，可以取xn+1=xm+1n+1。

10.3 预报-校正法#

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

　　这里采用显式欧拉法实现预报器，预报值x0n+1计算为

10.4 自适应步长控制算法#

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

　　由6.3节可知，hn的选取应使得第n步的局部截断误差ξ(n+1)L=|x(tn+1)−xn+1|满足：

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

　　自适应步长控制中，先设定一个足够小的初始步长h0，每进行一次迭代，就计算出本次迭代的局部截断误差ξ(n+1)L，再通过式(24)判定步长的好坏：

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

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

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

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

　　下图为瞬态仿真实例：

　　图中显示了步长自适应调整，时间轴是均匀的。

10.4.1 断点#

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

　　另一个问题出现在维持大步长时（如上图中的平稳部分）突然发生离散事件。在混合数字仿真中，波形存在大量间断点（电平跳变的瞬间）。步长过大会直接越过间断点。需要注意在事件驱动的数字仿真系统中，模拟仿真的误差估计并不会考虑离散事件造成的的跳变，这就是为什么自适应步长控制会失败。

　　因此，涉及数字-模拟混合仿真时，必须使用断点（simulink中称为过零检测）技术，在电平跳变之前插入断点。使模拟仿真器在时间到达断点前，强制减少积分步长，避免错过离散事件。对于模拟仿真，特别是涉及受控开关时，断点同样重要，可有效避免步长振荡。

　　稍微解释一下为什么数字仿真可以预测电平跳变。每个事件从进入队列到被调度执行，需间隔该事件指定的延时，因此事件的执行总是慢于入队。在事件入队时，就可预知断点的位置（入队时间+延迟），从而提前通知SPICE仿真器。

以后有机会可能补充离散-连续系统联合仿真的方法。

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

10.5.1 电容的时域模型#

　　电容的时域特性描述为：

　　利用线性多步法(21)得到迭代格式：

　　将上式展开并移项（b−1≠0），可得

　　geq从物理上可解释为等效电导，ieq从物理上可解释为等效电流源，于是得到了电容的等效线性化模型，如图4所示。

　　根据之前推出的阻抗元件支路对MNA各子矩阵的贡献（式6）和独立电流源支路对MNA各子矩阵贡献值（式9）可知对应MNA方程为：

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

　　同时应在程序中应设置标记，指出参数geq(n)和ieq(n)是需要迭代计算的。给定步长hn，初值u0=us(0)−ut(0)，根据上式生成MNA方程，可解出节点电压us(1)、ut(1)，以此类推。

10.5.2 电感的时域模型#

　　电感的时域特性描述为：

　　利用线性多步法(21)得到迭代格式：

　　整理并移项（b−1≠0），可得

　　从物理意义上看，电感可等效为动态电阻req与独立电压源ueq串联，如图5所示。

　　根据之前推出的独立电压源的贡献值（式8）和电流控制电压源支路对MNA各子矩阵贡献值（式10）可知对应MNA方程为：

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

　　同时应在程序中应设置标记，指出参数req(n)和ueq(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