Python实现轻量级SPICE内核：直流/交流/瞬态电路仿真

简介这是一份用纯Python实现的轻量级SPICE风格电路仿真工具面向电子工程学习者、电路教学人员以及需要在Python项目中嵌入仿真能力的开发者。它无需依赖外部EDA软件或C编译环境完整覆盖网表解析、电路拓扑建模、节点电压法矩阵构建并支持工作点、直流扫描、交流小信号与时域瞬态四类分析同时可处理二极管等非线性器件。资源包共20个文件约370KB以6个Python源码模块为核心分别承担语法解析、器件模型定义、分析算法封装与求解流程协调另含图形界面入口与启动主程序配套10个PNG图标用于打开、保存、仿真、撤销等按钮以及说明文档和缓存文件。目前已有87人学习下载。读者可借此理解SPICE核心机制用于教学演示、电路原理验证或作为二次开发与功能嵌入的参考实现。1. 用 Python 重写 SPICE 内核一个能跑直流/交流/瞬态分析的轻量级电路仿真工具做硬件验证的同行大概都有过这种体验板子焊完想快速验证一个分压偏置或 RC 滤波的频响打开商业 EDA 要等半天授权写个脚本又只能算线性网络碰到二极管、三极管就歇菜。我最近拆的这个 Python 电路仿真工具就是冲着这个空档来的——它用纯 Python 实现了一个轻量级 SPICE 内核支持直流工作点、交流小信号扫频和瞬态分析并且把非线性器件的牛顿迭代也做了进去。整个包不依赖商业求解器装完 numpy 和 scipy 就能跑源码结构清晰到可以当教学参考。适合两类人一是想搞懂 SPICE 内部到底怎么解方程的嵌入式/硬件工程师二是需要把电路仿真嵌进自己 Python 工具链、又不想背 ngspice 那套 C 接口的开发者。下面我按“它怎么算 → 怎么用 → 哪里会翻车”的顺序把这份资源拆开讲。2. 从网表到矩阵直流工作点与牛顿迭代的实现逻辑2.1 改进节点分析法MNA是怎么搭出方程组的SPICE 类工具的核心不是“画电路”而是把电路翻译成一个线性方程组G·x z然后交给数值库解。这个工具用的是改进节点分析法Modified Nodal AnalysisMNA比标准节点法多了一列“支路电流”未知量专门用来处理电压源和受控源——因为理想电压源的电流没法直接用节点电压表示。具体做法是先给每个节点编号地节点固定为 0然后遍历网表里的每个器件往一个(NM) × (NM)的矩阵里“盖章”stamp。电阻盖章到节点导纳位置电压源额外增加一行一列记录支路电流电容电感在直流分析里分别当开路和短路处理。这一步是整个仿真的地基盖错了后面全错。我一般会先拿一个纯电阻分压网络验证 MNA 矩阵对不对因为它的解析解手算就能得到。比如下面这个 3 节点网络import numpy as np def build_mna_resistors(nodes, resistors, vsources): nodes: 节点总数不含地 resistors: [(n1, n2, R), ...] n 为节点号0 表示地 vsources: [(n_plus, n_minus, V), ...] 返回 G 矩阵和 z 向量 n nodes m len(vsources) size n m G np.zeros((size, size)) z np.zeros(size) # 电阻盖章导纳 g 1/R 加到对应节点 for n1, n2, R in resistors: g 1.0 / R if n1 0: G[n1-1, n1-1] g if n2 0: G[n2-1, n2-1] g if n1 0 and n2 0: G[n1-1, n2-1] - g G[n2-1, n1-1] - g # 电压源盖章增加支路电流未知量 for idx, (np_, nm, V) in enumerate(vsources): row n idx if np_ 0: G[np_-1, row] 1 G[row, np_-1] 1 if nm 0: G[nm-1, row] - 1 G[row, nm-1] - 1 z[row] V return G, z # 验证12V 电源经 R11k 和 R22k 分压 G, z build_mna_resistors(2, [(1, 0, 1000), (1, 2, 2000)], [(1, 0, 12)]) x np.linalg.solve(G, z) print(节点1电压:, x[0], 节点2电压:, x[1])这段代码里nodes是不含地的节点数resistors用元组列表描述连接关系vsources每项多出一个电流未知量。跑出来节点 2 的电压应该是 8V12 × 2k/(1k2k)如果对不上先查节点编号有没有从 1 开始、地有没有写成 0。这个“盖章”思路和商业 SPICE 完全一致只是后者用稀疏矩阵和更复杂的排序优化。2.2 非线性器件牛顿-拉夫逊迭代怎么收敛二极管、BJT 这些器件的电流是电压的非线性函数没法一次解出来。工具用的是牛顿-拉夫逊迭代在当前猜测电压附近把非线性 I-V 曲线线性化求切线转成等效电导和电流源塞进 MNA 矩阵解一次再用新电压重新线性化直到电压变化小于收敛容差。二极管模型通常用肖克利方程I Is·(exp(V/(n·Vt)) - 1)其中Vt kT/q ≈ 25.85mV室温n是理想因子。线性化后等效电导gd Is/(n·Vt)·exp(V/(n·Vt))等效电流源Id I - gd·V。迭代时有个经典坑指数项在电压稍大时就会溢出所以工程实现里都会加电压钳位比如限制单步电压变化不超过 0.5V和exp参数上限。import numpy as np def diode_stamp(Vd, Is1e-14, n1.0, Vt0.02585): 返回二极管的等效电导 gd 和等效电流源 Id # 钳位防止 exp 溢出这是 SPICE 里的标准做法 Vcrit n * Vt * np.log(n * Vt / (np.sqrt(2) * Is)) if Vd Vcrit: Vd Vcrit exp_term np.exp(min(Vd / (n * Vt), 40)) # 上限 40 防溢出 Id Is * (exp_term - 1) gd Is / (n * Vt) * exp_term Ieq Id - gd * Vd return gd, Ieq def newton_solve_diode(Vsupply, R, tol1e-6, max_iter100): 简单电路电源 - 电阻 - 二极管到地求工作点 Vd 0.6 # 初始猜测太离谱会不收敛 for i in range(max_iter): gd, Ieq diode_stamp(Vd) # 节点方程(Vsupply - Vd)/R gd*Vd Ieq G 1.0/R gd z Vsupply/R - Ieq Vd_new z / G if abs(Vd_new - Vd) tol: print(f收敛于第 {i1} 次迭代Vd {Vd_new:.4f} V) return Vd_new Vd Vd_new raise RuntimeError(牛顿迭代未收敛检查初值或增加迭代次数) newton_solve_diode(5.0, 1000)Vcrit是临界电压超过它之后 exp 增长过快必须钳位max_iter给 100 是保险正常三极管电路十几次就收敛。如果报“未收敛”八成是初值给得太离谱或者电路本身没有直流解比如把二极管反接还串了个理想电流源。这个工具在非线性部分还支持简单的 BJT Ebers-Moll 模型思路一样只是多两个结要同时迭代。3. 交流小信号与瞬态分析从扫频到时间步进3.1 交流分析复数矩阵与自动频率扫描交流小信号分析的前提是先把非线性器件在工作点处线性化得到小信号等效电路然后用复数表示电抗——电容导纳jωC电感阻抗jωL。工具的做法是先跑一次直流工作点把每个非线性器件的跨导、输出电阻算出来替换成线性小信号模型再对每个频点组装复数 MNA 矩阵求解。频率扫描一般用对数间隔因为工程上关心的是十倍频程内的趋势。下面这段演示了 RC 低通在 10Hz 到 1MHz 的幅频响应计算import numpy as np def ac_sweep(R, C, f_start10, f_stop1e6, points50): 一阶 RC 低通节点1经R到节点2节点2经C到地节点1为输入 freqs np.logspace(np.log10(f_start), np.log10(f_stop), points) H [] for f in freqs: w 2 * np.pi * f # 复数 MNA节点2 的 KCL(V2-V1)/R jwC*V2 0 # 输入 V1 1单位激励解 V2 Y 1.0/R 1j * w * C V2 (1.0/R) / Y H.append(V2) H np.array(H) # 找 -3dB 点 mag_db 20 * np.log10(np.abs(H)) idx_3db np.argmin(np.abs(mag_db 3)) print(f截止频率约 {freqs[idx_3db]:.1f} Hz理论值 {1/(2*np.pi*R*C):.1f} Hz) return freqs, H freqs, H ac_sweep(1000, 100e-9)points决定扫频密度50 个点画波特图够用要精确定位谐振峰就加到 200 以上。注意这里输入固定为 1V算出来的是传递函数如果电路有多个独立源得逐个激励求叠加。工具内部对每个频点重新组装矩阵虽然比稀疏矩阵慢但胜在代码直观几千个节点以内完全能接受。3.2 瞬态分析时间步进与数值稳定性瞬态分析要解微分方程工具用的是后向欧拉或梯形法配合自适应步长。电容和电感在每个时间步被替换成等效电导和电流源伴随模型这样每个时间步又变成一个线性方程组。后向欧拉无条件稳定但精度一阶梯形法二阶精度但可能在陡变处振荡工具默认用梯形法遇到不连续点自动切后向欧拉。时间步长的选择是瞬态分析里最玄学的一环。步长太大高频成分被抹掉太小仿真时间爆炸。常见做法是取最小时间常数的 1/10 到 1/20同时限制单步电压变化。下面是一个 RC 充电的瞬态求解骨架import numpy as np def transient_rc(Vsupply, R, C, t_stop, dt): 后向欧拉求解 RC 充电返回时间轴和电容电压 n_steps int(t_stop / dt) t np.linspace(0, t_stop, n_steps 1) Vc np.zeros(n_steps 1) # 后向欧拉C*(Vc[n]-Vc[n-1])/dt (Vs - Vc[n])/R # 整理得 Vc[n] (Vc[n-1]*C/dt Vs/R) / (C/dt 1/R) a C / dt b 1.0 / R for n in range(1, n_steps 1): Vc[n] (Vc[n-1] * a Vsupply * b) / (a b) tau R * C print(f时间常数 {tau:.3e} sdt/tau {dt/tau:.3f}) print(f终值 {Vc[-1]:.4f} V理论 {Vsupply:.4f} V) return t, Vc t, Vc transient_rc(5.0, 1000, 1e-6, 5e-3, 1e-5)dt取tau/10时终值误差在千分之一量级取tau/100更准但步数翻十倍。工具内部会根据相邻步的电压变化自动调整dt变化大就缩小平稳就放大这个策略叫局部截断误差控制。如果你自己写循环记得每步都要重新组装矩阵——非线性器件的工作点随时间是变的。4. 避坑与排查仿真跑不通时先看这几处4.1 直流不收敛现象、原因与解决现象直流分析报“牛顿迭代未收敛”或结果全是 NaN。原因最常见的是初值给得太差比如所有节点电压初始化为 0而电路里有正反馈或高增益级其次是模型参数超出物理范围比如二极管Is设成 1e-3 这种离谱值。解决给关键节点加.ic初始条件工具支持在网表里写初始电压或者先断开反馈环跑一次开环工作点再接入闭环。我一般还会把max_iter临时调到 500看它到底卡在哪个节点。4.2 交流结果和手算对不上先查单位再查参考方向现象截止频率差了 1000 倍或者相位反了 180 度。原因九成是单位问题——电容写成100以为是 100nF实际按 100F 算了剩下一成是电流参考方向定义反了导致传递函数符号翻转。解决网表里强制用科学计数法写值100n而不是100并在交流分析后打印几个关键节点的幅值和理论分压比对照。这个工具对后缀p/n/u/m/k/meg的解析是大小写敏感的M是毫MEG才是兆写错一个字母结果差 1e9。4.3 瞬态波形振荡梯形法的锅还是步长太大现象方波激励下输出出现高频振铃但电路里根本没有电感。原因梯形法在激励突变时会产生数值振荡这是算法固有特性不是电路问题。解决在激励源边沿附近强制切后向欧拉或者把最大步长限制到上升时间的 1/10。工具里有个method参数可以全局切换积分方法调试阶段我通常先用后向欧拉确认趋势再用梯形法出最终波形。4.4 非线性器件模型参数缺失默认值不是万能药现象换了个型号的二极管仿真结果和手册曲线对不上。原因工具内置的默认模型参数Is1e-14, n1只适合通用小信号二极管功率二极管或肖特基的Is和Rs差好几个数量级。解决从器件手册里抄Is、n、Rs三个关键参数至少把串联电阻Rs加上否则大电流下的压降会偏小。这个工具支持在网表里用.model语句覆盖默认值格式和标准 SPICE 兼容。4.5 节点编号与地节点一个字母引发的血案现象电路明明有直流通路却报“奇异矩阵”。原因网表里地节点写成了GND而不是0工具没识别导致整个网络浮空。解决统一用0表示地所有需要接地的器件明确连到0。另外悬空节点也会导致矩阵奇异检查每个节点是否至少连接了两个器件。这个坑我踩过不止一次后来养成习惯网表写完后先跑一个“节点连通性检查”把只连了一个器件的节点全列出来。5. 把仿真嵌进自己的工具链参数化扫描与结果验证这份资源最实用的地方是它把求解器做成了可调用的 Python 类而不是一个只能跑网表的黑匣子。你可以直接在脚本里改器件参数、循环扫描、拿结果做优化。我常用的一个套路是参数化扫描固定电路拓扑让某个电阻或电容在范围内变化批量跑交流分析提取截止频率或谐振点画成曲线。import numpy as np from my_circuit_sim import Circuit, ACSweep # 假设的模块名按实际包结构调整 def sweep_cutoff(R_values, C_fixed100e-9): 扫描电阻提取 RC 低通的 -3dB 频率 results [] for R in R_values: ckt Circuit() ckt.add_resistor(1, 2, R) ckt.add_capacitor(2, 0, C_fixed) ckt.add_vsource(1, 0, 1.0) ac ACSweep(ckt, f_start10, f_stop1e6, points100) freqs, H ac.run() mag_db 20 * np.log10(np.abs(H)) idx np.argmin(np.abs(mag_db 3)) results.append((R, freqs[idx])) return results for R, fc in sweep_cutoff([500, 1000, 2000, 5000]): print(fR{R}Ω, fc{fc:.1f}Hz, 理论{1/(2*np.pi*R*100e-9):.1f}Hz)这段代码的关键是Circuit对象可以反复实例化每次改一个参数重新跑互不干扰。ACSweep的points在扫描时不用给太高50 到 100 足够定位 -3dB 点要精确值可以在粗扫后局部加密。验证方法很简单把扫描结果和解析公式fc 1/(2πRC)对比误差应该在 1% 以内如果偏差大先查扫频范围是否覆盖了截止点再查points是否太稀导致漏掉拐点。另一个进阶用法是把瞬态分析结果导出成 CSV喂给 matplotlib 或 pandas 做后处理。工具的输出接口返回的是 numpy 数组时间轴和节点电压一一对应直接np.savetxt就能存。我习惯在每次仿真后强制走一遍“能量守恒检查”对纯阻性电路输入功率应该等于各电阻耗散功率之和偏差超过 1e-6 就说明矩阵组装有问题。这个习惯帮我逮住过好几次盖章时符号写反的低级错误。从那以后我每次拿到一个新的电路仿真内核都会先用分压网络和 RC 充电这两个解析解已知的电路各跑一遍确认直流、交流、瞬态三条路径都正确再往复杂拓扑上堆。这份 Python 实现的轻量级工具胜在代码透明、改起来没有心理负担适合当自己工具链里的“计算器”用。希望帮到你。本文还有配套的精品资源点击获取