从支路电压到节点电压法

2026-06-27 2026-07-10 约 14300 字预计阅读 29 分钟 - 次阅读

上一篇文章里，我们一直把一个支路两端的电压写成 $u_k$ 。

这里的下标 $k$ 表示离散时刻，不是节点编号。换句话说， $u_k$ 的意思是“第 $k$ 个时间步的支路电压”。这一点很容易和节点电压法里的节点编号混在一起，所以从这篇文章开始，我们把两个概念分开写。

把支路两端电压先记成 $u$ ，是一种很适合入门的写法。读者可以先把注意力放在“一个元件或一个支路怎么离散”，不用一上来就面对完整电路网络里的节点编号、矩阵装配和右端向量符号约定。

但是实际电路仿真程序通常不会真的把每条支路的端口电压都当成独立未知量。更常见的做法是节点电压法：先选定参考地，再把每个非参考节点的电压作为未知量。支路两端电压不是新的未知量，而是两个节点电压之差。

也就是说，前面写的：

u

到了网络求解器里，通常要换成：

u_{ab}=u_a-u_b

其中 $a$ 和 $b$ 是节点编号， $u_a$ 和 $u_b$ 是当前时刻这两个节点对参考地的节点电压。若要显式写离散时刻，可以写成：

u_{ab}^{(k)}=u_a^{(k)}-u_b^{(k)}

这里上标 $(k)$ 表示第 $k$ 个离散时刻。为了让公式更轻，后文默认都在同一个当前时刻讨论，于是省略时间上标，只写：

u_{ab}=u_a-u_b

从一个 R 支路开始

先看最简单的电阻支路。设电阻连接在节点 $a$ 和节点 $b$ 之间，支路电流参考方向从 $a$ 流向 $b$ 。

电阻支路关系是：

i=G(u_a-u_b)

其中 $G=1/R$ 。

这句话已经很接近节点电压法了，因为右边只剩节点电压。对节点 $a$ 来说，这条支路从节点 $a$ 流出去的电流是：

i_a=G(u_a-u_b)

对节点 $b$ 来说，同一条支路从节点 $b$ 流出去的电流是反方向：

i_b=G(u_b-u_a)

把这两行放在一起：

\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G&-G\\-G&G\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}

这就是最普通的电阻支路对节点导纳矩阵的贡献。

它也解释了为什么电路程序喜欢节点电压法。只要每条支路电流都能写成节点电压的线性组合，程序就可以把每条支路贡献的 $2\times2$ 小矩阵装配到全局矩阵里。支路连接哪两个节点，就把这四个数加到哪两个节点对应的行列上。

L 支路

电阻支路只有当前节点电压，没有历史项。接下来换成电感支路，就会看到 Dommel 形式真正进入节点电压法的样子。

上一篇文章已经推过电感的梯形离散形式：

i=G_Lu+I_{L,\mathrm{his}}

其中：

G_L=\frac{T}{2L}

I_{L,\mathrm{his}}=i_{k-1}+G_Lu_{k-1}

这里的 $k-1$ 仍然表示上一拍离散时刻。现在把这条电感支路接在节点 $a$ 和节点 $b$ 之间，并规定支路电流参考方向从 $a$ 流向 $b$ ，支路电压参考方向也取：

u=u_a-u_b

于是支路电流就是：

i=G_L(u_a-u_b)+I_{L,\mathrm{his}}

现在要规定清楚正负号。本文采用下面这套约定：

支路电流参考方向从节点 $a$ 流向节点 $b$ 。
对每个节点写方程时，从该节点流出的支路电流记为正。
所以对节点 $a$ ，这条支路电流是流出节点，贡献为 $+i$ 。
对节点 $b$ ，同一条支路电流是流入节点，按“流出为正”的约定，贡献为 $-i$ 。

因此：

i_a=G_L(u_a-u_b)+I_{L,\mathrm{his}}

i_b=-i=G_L(u_b-u_a)-I_{L,\mathrm{his}}

把两行放到一起，就得到节点电压法里最常见的形式：

\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G_L&-G_L\\-G_L&G_L\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}+\begin{bmatrix}I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}\end{bmatrix}

这就是截图里那种结构。把电感换成任意已经整理好的 Dommel 支路，只要它能写成：

i=G_{\mathrm{eq}}u+I_{\mathrm{his}}

并且支路方向仍然从 $a$ 指向 $b$ ，就有：

\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G_{\mathrm{eq}}&-G_{\mathrm{eq}}\\-G_{\mathrm{eq}}&G_{\mathrm{eq}}\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}+\begin{bmatrix}I_{\mathrm{his}}\\-I_{\mathrm{his}}\end{bmatrix}

左边的两个量 $i_a,i_b$ ，表示这条支路分别对节点 $a$ 、节点 $b$ 的“流出电流贡献”。如果程序采用“流入节点为正”或者把所有已知项移到方程右端，历史项向量的符号会跟着变，但物理方向和支路方程本身没有变。

R+L+C 串联支路

现在回到上一篇文章最后的 RLC 串联支路。

支路层面已经得到：

i=G_{RLC}u+I_{RLC,\mathrm{his}}

其中：

G_{RLC}=\frac{1}{R+\frac{2L}{T}+\frac{T}{2C}}

I_{RLC,\mathrm{his}}=-G_{RLC}\left[u_{C,k-1}-u_{L,k-1}+\left(\frac{T}{2C}-\frac{2L}{T}\right)i_{k-1}\right]

注意这里的 $k-1$ 是离散时刻，不是节点编号。本文统一用 $a,b$ 表示节点编号，用 $u_a,u_b$ 表示节点电压，避免把“时间步下标”和“节点编号下标”混在一起。

如果 RLC 串联支路连接节点 $a$ 和节点 $b$ ，并规定支路电流从 $a$ 流向 $b$ ，那么：

u=u_a-u_b

代入支路 Dommel 形式：

i=G_{RLC}(u_a-u_b)+I_{RLC,\mathrm{his}}

于是它对两个节点的贡献就是：

\begin{bmatrix}i_a\\i_b\end{bmatrix}=\begin{bmatrix}G_{RLC}&-G_{RLC}\\-G_{RLC}&G_{RLC}\end{bmatrix}\begin{bmatrix}u_a\\u_b\end{bmatrix}+\begin{bmatrix}I_{RLC,\mathrm{his}}\\-I_{RLC,\mathrm{his}}\end{bmatrix}

这和电阻支路的矩阵结构完全一样。不同的只是：

电阻支路的导纳 $G$ 是常数。
RLC 串联支路的等效导纳 $G_{RLC}$ 来自梯形离散后的伴随模型。
RLC 串联支路多了历史项 $I_{RLC,\mathrm{his}}$ ，它要进入当前步的右端向量。

所以从网络求解器的角度看，Dommel 形式的好处就在这里：动态支路虽然内部有电感、电容和历史状态，但在当前时刻装配矩阵时，它仍然像一个“等效导纳加已知源”的二端支路。

展开 RLC 的内部节点

上面的 RLC 写法，是把整个 R+L+C 串联支路压缩成一个二端支路，只看外部两个节点 $a$ 和 $b$ 。这一步对网络求解器很友好，因为它只需要给节点 $a$ 和节点 $b$ 装配一个 $2\times2$ 小矩阵。

现在把难度加一点：不再把 RLC 压缩成一个整体，而是把 R、L、C 三个元件各自看成一条单独支路。这样串联支路内部会多出两个节点。设节点顺序从左到右为：

a\rightarrow c\rightarrow d\rightarrow b

其中：

R 支路连接节点 $a$ 和节点 $c$ 。
L 支路连接节点 $c$ 和节点 $d$ 。
C 支路连接节点 $d$ 和节点 $b$ 。

三个支路的电压分别是：

u_R=u_a-u_c

u_L=u_c-u_d

u_C=u_d-u_b

支路方向都取从左到右，于是三条支路可以写成：

i_R=G_R(u_a-u_c)

i_L=G_L(u_c-u_d)+I_{L,\mathrm{his}}

i_C=G_C(u_d-u_b)+I_{C,\mathrm{his}}

这里 $G_R=1/R$ ，而 $G_L$ 、 $G_C$ 来自电感和电容各自的 Dommel 伴随模型。历史项的正负号仍然沿用前面的约定：支路方向从左到右，流出节点为正。

接下来装配矩阵时，先固定一个全局节点顺序：

\mathbf u=\begin{bmatrix}u_a\\u_c\\u_d\\u_b\end{bmatrix}

后面的矩阵行、列都按照这个顺序来排。也就是说，第 1 行和第 1 列对应节点 $a$ ，第 2 行和第 2 列对应节点 $c$ ，第 3 行和第 3 列对应节点 $d$ ，第 4 行和第 4 列对应节点 $b$ 。

先看 R 支路。它连接节点 $a$ 和节点 $c$ ，所以它原本只有一个 $2\times2$ 小矩阵：

\begin{bmatrix}G_R&-G_R\\-G_R&G_R\end{bmatrix}

嵌入到全局节点顺序 $a,c,d,b$ 里，就变成：

G^{(R)}=\begin{bmatrix}G_R&-G_R&0&0\\-G_R&G_R&0&0\\0&0&0&0\\0&0&0&0\end{bmatrix}

意思是：R 只连接 $a$ 和 $c$ ，所以只会影响 $a,c$ 对应的四个位置，和 $d,b$ 没有关系。

再看 L 支路。它连接节点 $c$ 和节点 $d$ ，等效导纳贡献嵌入到全局矩阵以后是：

G^{(L)}=\begin{bmatrix}0&0&0&0\\0&G_L&-G_L&0\\0&-G_L&G_L&0\\0&0&0&0\end{bmatrix}

同时，电感 Dommel 形式里还有历史电流源。因为 L 支路方向取从 $c$ 到 $d$ ，按“流出节点为正”的约定，它对节点 $c$ 是 $+I_{L,\mathrm{his}}$ ，对节点 $d$ 是 $-I_{L,\mathrm{his}}$ ：

\mathbf I_{\mathrm{his}}^{(L)}=\begin{bmatrix}0\\I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}\\0\end{bmatrix}

最后看 C 支路。它连接节点 $d$ 和节点 $b$ ，所以等效导纳贡献是：

G^{(C)}=\begin{bmatrix}0&0&0&0\\0&0&0&0\\0&0&G_C&-G_C\\0&0&-G_C&G_C\end{bmatrix}

C 支路方向取从 $d$ 到 $b$ ，所以历史电流源对节点 $d$ 是 $+I_{C,\mathrm{his}}$ ，对节点 $b$ 是 $-I_{C,\mathrm{his}}$ ：

\mathbf I_{\mathrm{his}}^{(C)}=\begin{bmatrix}0\\0\\I_{C,\mathrm{his}}\\-I_{C,\mathrm{his}}\end{bmatrix}

三条支路相加，才得到总的节点导纳矩阵：

G=G^{(R)}+G^{(L)}+G^{(C)}=\begin{bmatrix}G_R&-G_R&0&0\\-G_R&G_R+G_L&-G_L&0\\0&-G_L&G_L+G_C&-G_C\\0&0&-G_C&G_C\end{bmatrix}

历史项也一样相加：

\mathbf I_{\mathrm{his}}=\mathbf I_{\mathrm{his}}^{(L)}+\mathbf I_{\mathrm{his}}^{(C)}=\begin{bmatrix}0\\I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}\\-I_{C,\mathrm{his}}\end{bmatrix}

所以完整的装配结果是：

\begin{bmatrix}i_a\\i_c\\i_d\\i_b\end{bmatrix}=\begin{bmatrix}G_R&-G_R&0&0\\-G_R&G_R+G_L&-G_L&0\\0&-G_L&G_L+G_C&-G_C\\0&0&-G_C&G_C\end{bmatrix}\begin{bmatrix}u_a\\u_c\\u_d\\u_b\end{bmatrix}+\begin{bmatrix}0\\I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}\\-I_{C,\mathrm{his}}\end{bmatrix}

这样看，中间节点 $d$ 的历史项为什么是：

-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}

就比较清楚了。节点 $d$ 左边接着 L，右边接着 C。L 支路方向是从 $c$ 到 $d$ ，所以 L 的历史电流对 $d$ 来说是流入，按“流出为正”记成负号；C 支路方向是从 $d$ 到 $b$ ，所以 C 的历史电流对 $d$ 来说是流出，记成正号。

这也解释了为什么程序里常说“装配”。所谓装配，不是重新推一遍全电路方程，而是每条支路只关心自己两端节点，把自己的小矩阵和历史项加到全局矩阵、全局向量对应的位置上。记住这个规则就够了：一条连接节点 $p,q$ 的二端支路，会把 $+G$ 加到 $(p,p)$ 和 $(q,q)$ ，把 $-G$ 加到 $(p,q)$ 和 $(q,p)$ ；如果它还有历史源，就按支路方向把 $+I_{\mathrm{his}}$ 加到起点节点，把 $-I_{\mathrm{his}}$ 加到终点节点。

实战：正弦电压源驱动 RLC 支路

最后做一个很小的实战。为了不把问题一下子推进到改进节点电压法，这里先采用最简单的电路：理想正弦电压源一端接地，另一端接 RLC 串联支路。

这时右端节点是参考地。读者可以思考一下：为什么这里必须加上参考地？

u_b=0

左端节点由电压源强制给定：

u_a=v_s(t)

所以 RLC 支路两端电压不是未知量，而是已知量：

u=u_a-u_b=v_s(t)

这一步很重要。这个例子不是在求一个含浮接理想电压源的完整节点矩阵，而是在已知端口电压的条件下，计算 RLC 串联支路的动态响应。也就是说，节点电压法里的“参考地”和“已知节点电压”已经先处理好了。

为了和上一篇文章衔接，先保留一种压缩支路写法：把整个 RLC 串联支路看成一个二端口，直接用已知的端口电压 $u=u_a-u_b$ 递推出支路电流。这个写法更短，但它还没有真正体现本章前面讲的“把 R、L、C 分别装配进节点导纳矩阵”。

取：

v_s(t)=100\sin(2\pi 50t)

示例参数为：

R=10\ \Omega,\qquad L=20\ \mathrm{mH},\qquad C=100\ \mu\mathrm{F},\qquad T=20\ \mu\mathrm{s}

计算流程可以按下面理解：

初始时刻默认 $u_C(0)=0$ 、 $i(0)=0$ ，因此第一拍的历史量来自零初值。
每一个当前步先由电压源得到 $u_k=v_s(t_k)$ 。
再把上一拍保存的 $u_{L,k-1}$ 、 $u_{C,k-1}$ 、 $i_{k-1}$ 合成历史项，求出当前步 $i_k$ 。
最后用当前步 $i_k$ 更新 $u_{L,k}$ 和 $u_{C,k}$ 。下一拍再计算历史项时，用的就是这组刚刚更新过的状态。

所以程序里看起来是在更新 $u_L$ 、 $u_C$ ，本质上是在保存电感和电容的内部状态。 $I_{\mathrm{his}}$ 不是一个必须单独长期保存的神秘变量；它可以在每一拍开始时，由上一拍状态重新算出来。

压缩支路写法的完整示例代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
import numpy as np
import matplotlib.pyplot as plt

R = 10.0
L = 20e-3
C = 100e-6
f = 50.0
U_peak = 100.0
T = 20e-6
t_end = 0.12

t = np.arange(0.0, t_end + T, T)
u = U_peak * np.sin(2.0 * np.pi * f * t)
i = np.zeros_like(t)
u_L = np.zeros_like(t)
u_C = np.zeros_like(t)

G_rlc = 1.0 / (R + 2.0 * L / T + T / (2.0 * C))

for k in range(1, len(t)):
    history = u_C[k - 1] - u_L[k - 1] + (T / (2.0 * C) - 2.0 * L / T) * i[k - 1]
    i[k] = G_rlc * (u[k] - history)
    u_L[k] = 2.0 * L / T * (i[k] - i[k - 1]) - u_L[k - 1]
    u_C[k] = u_C[k - 1] + T / (2.0 * C) * (i[k] + i[k - 1])

u_R = R * i

fig, (ax1, ax3) = plt.subplots(1, 2, figsize=(11.0, 3.8), sharex=True)

ax1.plot(t * 1000.0, u, color="#2563eb", linewidth=1.8, label="source voltage")
ax1.set_xlabel("time / ms")
ax1.set_ylabel("voltage / V", color="#2563eb")
ax1.tick_params(axis="y", labelcolor="#2563eb")
ax1.grid(True, color="#e2e8f0", linewidth=0.8)

ax2 = ax1.twinx()
ax2.plot(t * 1000.0, i, color="#dc2626", linewidth=1.8, label="branch current")
ax2.set_ylabel("current / A", color="#dc2626")
ax2.tick_params(axis="y", labelcolor="#dc2626")

lines = ax1.get_lines() + ax2.get_lines()
labels = [line.get_label() for line in lines]
ax1.legend(lines, labels, loc="upper right", frameon=False)

ax3.plot(t * 1000.0, u_R, color="#16a34a", linewidth=1.5, label="u_R")
ax3.plot(t * 1000.0, u_L, color="#9333ea", linewidth=1.5, label="u_L")
ax3.plot(t * 1000.0, u_C, color="#f97316", linewidth=1.5, label="u_C")
ax3.set_xlabel("time / ms")
ax3.set_ylabel("component voltage / V")
ax3.grid(True, color="#e2e8f0", linewidth=0.8)
ax3.legend(loc="upper right", frameon=False)

fig.tight_layout()
fig.savefig("rlc-voltage-source-result.svg", format="svg")

仿真结果如下：

左图中蓝色曲线是外加电压源，红色曲线是 RLC 串联支路电流；右图画出同一时刻 R、L、C 三个元件各自承担的电压。由于这里的电容、电感初值都取零，刚开始会有一段暂态；经过几个周波以后，响应逐渐进入 50 Hz 正弦稳态。

不过，从本章标题来看，这段代码还差一步：它仍然是在“支路电压已知”的层面算，而不是在程序里显式写出节点矩阵。下面把同一个电路改成节点电压法的写法。

用节点电压矩阵再写一遍

仍然采用节点顺序：

\mathbf u=\begin{bmatrix}u_a\\u_c\\u_d\\u_b\end{bmatrix}

为了和前面的 $u_a,u_c,u_d,u_b$ 区别开，下面把节点电压向量统一写成大写 $\mathbf V$ ：

\mathbf V=\begin{bmatrix}u_a\\u_c\\u_d\\u_b\end{bmatrix}

第一步，先写出完整的节点方程。

按照前面装配出来的矩阵，每一个时间步仍然保持这种形式：

\mathbf I=G\mathbf V+\mathbf I_{\mathrm{his}}

其中：

G=\begin{bmatrix}G_R&-G_R&0&0\\-G_R&G_R+G_L&-G_L&0\\0&-G_L&G_L+G_C&-G_C\\0&0&-G_C&G_C\end{bmatrix}

\mathbf I_{\mathrm{his}}=\begin{bmatrix}0\\I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}\\-I_{C,\mathrm{his}}\end{bmatrix}

这里 $\mathbf I$ 表示各节点的外部注入电流， $\mathbf I_{\mathrm{his}}$ 表示电感、电容历史项折算出来的已知电流项。

这里暂时只关心一个和初值有关的问题。第一拍计算时，历史项来自初始状态：

I_{L,\mathrm{his}}=i_{L,k-1}+G_Lu_{L,k-1}

I_{C,\mathrm{his}}=-G_Cu_{C,k-1}-i_{C,k-1}

若默认零初值，就是：

i_L(0)=0,\qquad u_C(0)=0

若用户给了非零初值，例如：

i_L(0)=2\ \mathrm A,\qquad u_C(0)=50\ \mathrm V

那么第一拍的 $I_{L,\mathrm{his}}$ 和 $I_{C,\mathrm{his}}$ 就会随之改变，仿真曲线一开始的暂态也会不同。

这里还有一个小边界：这三个元件是串联的，所以不能随便令 $i_L(0)$ 和 $i_C(0)$ 不相等。更自然的初值是指定同一个支路电流 $i_0$ 和电容电压 $u_C(0)$ ；电感初始电压再由这一拍的 KVL 关系确定。

第二步，分清哪些量已知，哪些量未知。

这个小电路里，左端节点由理想电压源直接给定，右端节点是参考地：

u_a=v_s(t),\qquad u_b=0

所以已知节点电压是：

\mathbf V_{\mathrm{known}}=\begin{bmatrix}u_a\\u_b\end{bmatrix}

真正要求的是内部节点：

\mathbf V_{\mathrm{unknown}}=\begin{bmatrix}u_c\\u_d\end{bmatrix}

内部节点 $c,d$ 没有额外外部电流注入，所以它们对应的节点电流为零：

\mathbf I_{\mathrm{unknown}}=\begin{bmatrix}0\\0\end{bmatrix}

第三步，只取内部节点对应的两行方程。

从完整 $4\times4$ 矩阵里取出第 2、3 行，也就是节点 $c,d$ 对应的方程：

\begin{bmatrix}0\\0\end{bmatrix}= \begin{bmatrix} -G_R&G_R+G_L&-G_L&0\\ 0&-G_L&G_L+G_C&-G_C \end{bmatrix} \begin{bmatrix}u_a\\u_c\\u_d\\u_b\end{bmatrix} +\begin{bmatrix}I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}\end{bmatrix}

再把这一行按“未知节点 $c,d$ ”和“已知节点 $a,b$ ”拆开：

\begin{bmatrix}0\\0\end{bmatrix}= \underbrace{\begin{bmatrix} G_R+G_L&-G_L\\ -G_L&G_L+G_C \end{bmatrix}}_{G_{uu}} \underbrace{\begin{bmatrix}u_c\\u_d\end{bmatrix}}_{\mathbf V_{\mathrm{unknown}}} + \underbrace{\begin{bmatrix} -G_R&0\\ 0&-G_C \end{bmatrix}}_{G_{uk}} \underbrace{\begin{bmatrix}u_a\\u_b\end{bmatrix}}_{\mathbf V_{\mathrm{known}}} + \underbrace{\begin{bmatrix}I_{L,\mathrm{his}}\\-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}\end{bmatrix}}_{\mathbf I_{\mathrm{his},u}}

第四步，把已知项移到右边。

于是内部节点电压满足：

G_{uu}\mathbf V_{\mathrm{unknown}}=-\left(G_{uk}\mathbf V_{\mathrm{known}}+\mathbf I_{\mathrm{his},u}\right)

程序里不需要真的写矩阵逆，而是用np.linalg.solve求解这个小线性方程组。这样写也更接近真实程序：左边是导纳矩阵，右边是已知电压和历史电流项整理出来的右端向量。完整代码如下：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
import numpy as np
import matplotlib.pyplot as plt


R = 10.0
L = 20e-3
C = 100e-6
f = 50.0
U_peak = 100.0
T = 20e-6
t_end = 0.12

t = np.arange(0.0, t_end + T, T)
source = U_peak * np.sin(2.0 * np.pi * f * t)

G_R = 1.0 / R
G_L = T / (2.0 * L)
G_C = 2.0 * C / T

G = np.array(
    [
        [G_R, -G_R, 0.0, 0.0],
        [-G_R, G_R + G_L, -G_L, 0.0],
        [0.0, -G_L, G_L + G_C, -G_C],
        [0.0, 0.0, -G_C, G_C],
    ]
)

unknown_nodes = np.array([1, 2])  # u_c, u_d
known_nodes = np.array([0, 3])  # u_a, u_b
G_uu = G[np.ix_(unknown_nodes, unknown_nodes)]
G_uk = G[np.ix_(unknown_nodes, known_nodes)]


def simulate(i0=0.0, u_C0=0.0):
    u_a = source.copy()
    u_b = np.zeros_like(t)
    u_c = np.zeros_like(t)
    u_d = np.zeros_like(t)

    u_R = np.zeros_like(t)
    u_L = np.zeros_like(t)
    u_C = np.zeros_like(t)
    i_R = np.zeros_like(t)
    i_L = np.zeros_like(t)
    i_C = np.zeros_like(t)

    i_R[0] = i0
    i_L[0] = i0
    i_C[0] = i0
    u_R[0] = R * i0
    u_C[0] = u_C0
    u_L[0] = u_a[0] - u_b[0] - u_R[0] - u_C[0]
    u_d[0] = u_b[0] + u_C[0]
    u_c[0] = u_d[0] + u_L[0]

    for k in range(1, len(t)):
        I_L_his = i_L[k - 1] + G_L * u_L[k - 1]
        I_C_his = -G_C * u_C[k - 1] - i_C[k - 1]

        I_his = np.array([0.0, I_L_his, -I_L_his + I_C_his, -I_C_his])
        V_known = np.array([u_a[k], u_b[k]])

        rhs = -(G_uk @ V_known + I_his[unknown_nodes])
        V_unknown = np.linalg.solve(G_uu, rhs)

        u_c[k], u_d[k] = V_unknown

        u_R[k] = u_a[k] - u_c[k]
        u_L[k] = u_c[k] - u_d[k]
        u_C[k] = u_d[k] - u_b[k]

        i_R[k] = G_R * u_R[k]
        i_L[k] = G_L * u_L[k] + I_L_his
        i_C[k] = G_C * u_C[k] + I_C_his

    return {
        "u_a": u_a,
        "u_b": u_b,
        "u_c": u_c,
        "u_d": u_d,
        "u_R": u_R,
        "u_L": u_L,
        "u_C": u_C,
        "i_R": i_R,
        "i_L": i_L,
        "i_C": i_C,
    }


zero = simulate()
nonzero = simulate(i0=2.0, u_C0=50.0)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(13.0, 3.8), sharex=True)

ax1.plot(t * 1000.0, zero["u_a"], color="#2563eb", linewidth=1.6, label="u_a")
ax1.plot(t * 1000.0, zero["u_c"], color="#16a34a", linewidth=1.4, label="u_c")
ax1.plot(t * 1000.0, zero["u_d"], color="#f97316", linewidth=1.4, label="u_d")
ax1.plot(t * 1000.0, zero["u_b"], color="#64748b", linewidth=1.2, label="u_b")
ax1.set_xlabel("time / ms")
ax1.set_ylabel("node voltage / V")
ax1.grid(True, color="#e2e8f0", linewidth=0.8)
ax1.legend(loc="upper right", frameon=False)

ax2.plot(t * 1000.0, zero["i_R"], color="#2563eb", linewidth=1.5, label="i_R")
ax2.plot(t * 1000.0, zero["i_L"], color="#16a34a", linewidth=1.5, linestyle="--", label="i_L")
ax2.plot(t * 1000.0, zero["i_C"], color="#dc2626", linewidth=1.2, linestyle=":", label="i_C")
ax2.set_xlabel("time / ms")
ax2.set_ylabel("branch current / A")
ax2.grid(True, color="#e2e8f0", linewidth=0.8)
ax2.legend(loc="upper right", frameon=False)

ax3.plot(t * 1000.0, zero["u_R"], color="#16a34a", linewidth=1.4, label="u_R")
ax3.plot(t * 1000.0, zero["u_L"], color="#9333ea", linewidth=1.4, label="u_L")
ax3.plot(t * 1000.0, zero["u_C"], color="#f97316", linewidth=1.4, label="u_C")
ax3.set_xlabel("time / ms")
ax3.set_ylabel("component voltage / V")
ax3.grid(True, color="#e2e8f0", linewidth=0.8)
ax3.legend(loc="upper right", frameon=False)

fig.tight_layout()
fig.savefig("rlc-nodal-voltage-result.svg", format="svg")

fig2, (bx1, bx2) = plt.subplots(1, 2, figsize=(11.0, 3.8), sharex=True)

bx1.plot(t * 1000.0, zero["i_L"], color="#2563eb", linewidth=1.5, label="zero initial")
bx1.plot(
    t * 1000.0,
    nonzero["i_L"],
    color="#dc2626",
    linewidth=1.5,
    label="i_L(0)=2 A, u_C(0)=50 V",
)
bx1.set_xlabel("time / ms")
bx1.set_ylabel("branch current / A")
bx1.grid(True, color="#e2e8f0", linewidth=0.8)
bx1.legend(loc="upper right", frameon=False)

bx2.plot(t * 1000.0, zero["u_C"], color="#2563eb", linewidth=1.5, label="zero initial")
bx2.plot(
    t * 1000.0,
    nonzero["u_C"],
    color="#dc2626",
    linewidth=1.5,
    label="i_L(0)=2 A, u_C(0)=50 V",
)
bx2.set_xlabel("time / ms")
bx2.set_ylabel("capacitor voltage / V")
bx2.grid(True, color="#e2e8f0", linewidth=0.8)
bx2.legend(loc="upper right", frameon=False)

fig2.tight_layout()
fig2.savefig("rlc-nodal-initial-comparison.svg", format="svg")

这个版本和前一个版本算的是同一个电路，但写法更贴近节点电压法。它每一拍都做了几件事：

先由上一拍状态算出 $I_{L,\mathrm{his}}$ 和 $I_{C,\mathrm{his}}$ 。
组装历史项向量 $\mathbf I_{\mathrm{his}}$ 。
把已知节点 $u_a,u_b$ 移到右端。
用节点矩阵求出内部节点 $u_c,u_d$ 。
再由节点电压差回算 $u_R,u_L,u_C$ 和三条支路电流。

仿真结果如下：

节点电压法求解正弦电压源驱动 RLC 支路 — 用节点电压矩阵求解同一个 RLC 串联支路

中间图里 $i_R$ 、 $i_L$ 、 $i_C$ 基本重合，因为这三个元件串联在同一条支路上。这个图反过来也可以当成一个小检查：如果三条支路电流明显对不上，大概率就是历史项符号、节点顺序或矩阵装配位置弄错了。

如果把初值改成：

i_L(0)=2\ \mathrm A,\qquad u_C(0)=50\ \mathrm V

同一个程序会得到下面这种对比：

不同初值下的 RLC 节点电压法仿真结果比较 — 改变电感电流初值和电容电压初值后，暂态响应会明显改变

蓝色曲线是默认零初值，红色曲线是非零初值。可以看到，稳态频率仍然由外加正弦源和电路参数决定，但刚开始那段暂态明显不一样。换句话说，初值不是“画图时随便给个起点”这么简单；在离散 Dommel 形式里，它会进入第一拍的历史项，然后通过后续递推继续影响一段时间。

读者的思考

很多资料讲 Dommel 形式，讲到：

i=Gv+I_{\mathrm{his}}

就停下来了。这个式子当然重要，但真正写程序时，很快会遇到一堆它没有直接回答的问题。

比如，刚刚的例子，求解的电路为什么一定要有参考地？如果一个网络完全浮起来，不接地，节点电压矩阵会发生什么？再比如，理想电压源如果不是一端接地，而是浮接在两个非参考节点之间，还能不能直接塞进普通节点导纳矩阵？

还有一个很容易被初学者忽略的问题：串联支路里的初值不能随便乱给。理想串联 R-L-C 只有一条电流路径，所以同一时刻应当有：

i_R(0)=i_L(0)=i_C(0)

但电感真正的状态量是电流，电容真正的状态量是电压。更自然的初值通常是：

i_L(0)=i_0

u_C(0)=u_{C0}

而不是再额外指定一个和 $i_L(0)$ 不一致的“电容电流初值”。电容电流由电路方程决定：

i_C=C\frac{\mathrm du_C}{\mathrm dt}

如果在程序里硬给 $i_L(0)\ne i_C(0)$ ，就相当于一开始在内部节点制造了一个不满足 KCL 的矛盾。放到本文的节点 $d$ 上，就是历史项：

-I_{L,\mathrm{his}}+I_{C,\mathrm{his}}

突然包含了一个不该有的净注入。工程上通常会令：

i_L(0)=i_C(0)=i_0

再给出电容电压初值 $u_C(0)=u_{C0}$ ，并由 KVL 算出初始电感电压：

u_L(0)=u_a(0)-u_b(0)-Ri_0-u_C(0)

也就是说，串联支路可以自由指定电感电流初值和电容电压初值，但不能再独立指定一个不一致的电容电流初值。

这篇文章先不深究这些问题。本文只做一件事：把单条支路的 Dommel 形式，接到节点电压差和节点矩阵装配上。至于参考地、浮接电压源、改进节点电压法、理想源约束方程、矩阵奇异这些问题，后面单独开文章再慢慢拆。