在 EMT 仿真里，电感和电容通常会先被离散成 companion model，也就是当前时间步的等效导纳和历史电流源。这个写法非常适合进入节点电压方程，也是 Dommel 类算法里最常见的元件建模方式。
但如果 L L L C C C u = L d i / d t u=L\,di/dt u = L d i / d t i = C d u / d t i=C\,du/dt i = C d u / d t
ψ = L i , q = C u
\psi=L i,\qquad q=C u
ψ = L i , q = C u 这篇文章先从常见 LC 建模方式讲起，再说明为什么 EMT 仿真里有时会推崇磁链、电荷建模，最后从 d ( x y ) d t \frac{d(xy)}{dt} d t d ( x y ) θ \theta θ
常见 LC 建模方式
对普通电感和电容，最常见的起点是：
u L = L d i d t
u_L=L\frac{di}{dt}
u L = L d t d i i C = C d u d t
i_C=C\frac{du}{dt}
i C = C d t d u 这里默认 L L L C C C
i n + 1 = G u n + 1 + I h i s
i_{n+1}=G u_{n+1}+I_{\mathrm{his}}
i n + 1 = G u n + 1 + I his 其中，G G G I h i s I_{\mathrm{his}} I his
先看电容。对 i = C d u / d t i=C\,du/dt i = C d u / d t θ \theta θ
u n + 1 = u n + Δ t [ θ i n + 1 C + ( 1 − θ ) i n C ]
u_{n+1}=u_n+\Delta t\left[\theta \frac{i_{n+1}}{C}+(1-\theta)\frac{i_n}{C}\right]
u n + 1 = u n + Δ t [ θ C i n + 1 + ( 1 − θ ) C i n ] 整理得到：
i n + 1 = C θ Δ t u n + 1 − C θ Δ t u n − 1 − θ θ i n
i_{n+1}
=\frac{C}{\theta\Delta t}u_{n+1}
-\frac{C}{\theta\Delta t}u_n
-\frac{1-\theta}{\theta}i_n
i n + 1 = θ Δ t C u n + 1 − θ Δ t C u n − θ 1 − θ i n 因此可以写成：
i n + 1 = G C u n + 1 + I h i s , C
i_{n+1}=G_C u_{n+1}+I_{\mathrm{his},C}
i n + 1 = G C u n + 1 + I his , C 其中：
G C = C θ Δ t
G_C=\frac{C}{\theta\Delta t}
G C = θ Δ t C I h i s , C = − C θ Δ t u n − 1 − θ θ i n
I_{\mathrm{his},C}
=-\frac{C}{\theta\Delta t}u_n
-\frac{1-\theta}{\theta}i_n
I his , C = − θ Δ t C u n − θ 1 − θ i n 再看电感。对 u = L d i / d t u=L\,di/dt u = L d i / d t θ \theta θ
i n + 1 = i n + Δ t L [ θ u n + 1 + ( 1 − θ ) u n ]
i_{n+1}=i_n+\frac{\Delta t}{L}\left[\theta u_{n+1}+(1-\theta)u_n\right]
i n + 1 = i n + L Δ t [ θ u n + 1 + ( 1 − θ ) u n ] 整理得到：
i n + 1 = θ Δ t L u n + 1 + i n + ( 1 − θ ) Δ t L u n
i_{n+1}
=\frac{\theta\Delta t}{L}u_{n+1}
+i_n+\frac{(1-\theta)\Delta t}{L}u_n
i n + 1 = L θ Δ t u n + 1 + i n + L ( 1 − θ ) Δ t u n 同样写成：
i n + 1 = G L u n + 1 + I h i s , L
i_{n+1}=G_L u_{n+1}+I_{\mathrm{his},L}
i n + 1 = G L u n + 1 + I his , L 其中：
G L = θ Δ t L
G_L=\frac{\theta\Delta t}{L}
G L = L θ Δ t I h i s , L = i n + ( 1 − θ ) Δ t L u n
I_{\mathrm{his},L}
=i_n+\frac{(1-\theta)\Delta t}{L}u_n
I his , L = i n + L ( 1 − θ ) Δ t u n 这里的重点不是公式长什么样，而是结构：当前拍未知量只出现在 u n + 1 u_{n+1} u n + 1
这就是常见 LC 建模的基本形式。
（如果读者还不熟悉 θ \theta θ 《梯形积分为什么会数值震荡》 ；那篇文章从数值阻尼和高频误差的角度解释了为什么工程 EMT 里经常在梯形法和 θ \theta θ ）
磁链和电荷建模
从物理状态看，电感真正存储的变量不是电流本身，而是磁链：
ψ = L i
\psi=L i
ψ = L i 端口电压由磁链变化率决定：
u = d ψ d t
u=\frac{d\psi}{dt}
u = d t d ψ 电容真正存储的变量不是电压本身，而是电荷：
q = C u
q=C u
q = C u 端口电流由电荷变化率决定：
i = d q d t
i=\frac{dq}{dt}
i = d t d q 从这个角度看，电力电子里常说的“伏秒平衡”其实就是磁链平衡。因为：
Δ ψ = ∫ t 0 t 1 u ( t ) d t
\Delta \psi=\int_{t_0}^{t_1}u(t)\,dt
Δ ψ = ∫ t 0 t 1 u ( t ) d t 如果电感在一个开关周期或一个稳态周期内回到同一个磁链状态，就有：
∫ t 0 t 0 + T u ( t ) d t = 0
\int_{t_0}^{t_0+T}u(t)\,dt=0
∫ t 0 t 0 + T u ( t ) d t = 0 这就是伏秒平衡。它不是额外假设，而是 u = d ψ / d t u=d\psi/dt u = d ψ / d t
Δ q = ∫ t 0 t 1 i ( t ) d t
\Delta q=\int_{t_0}^{t_1}i(t)\,dt
Δ q = ∫ t 0 t 1 i ( t ) d t 周期稳态下电容电荷不漂移，因此一个周期内电容电流积分为零。
当 L L L C C C
u = d ( L i ) d t = L d i d t
u=\frac{d(Li)}{dt}=L\frac{di}{dt}
u = d t d ( L i ) = L d t d i i = d ( C u ) d t = C d u d t
i=\frac{d(Cu)}{dt}=C\frac{du}{dt}
i = d t d ( C u ) = C d t d u EMT 仿真有时推崇从磁链、电荷出发建模，原因主要有两个。
第一，状态变量更接近物理本质。对电感来说，连续性首先作用在磁链 ψ \psi ψ q q q L L L C C C
当 L L L C C C
u = d ( L i ) d t = L d i d t + i d L d t
u=\frac{d(Li)}{dt}
=L\frac{di}{dt}+i\frac{dL}{dt}
u = d t d ( L i ) = L d t d i + i d t d L i = d ( C u ) d t = C d u d t + u d C d t
i=\frac{d(Cu)}{dt}
=C\frac{du}{dt}+u\frac{dC}{dt}
i = d t d ( C u ) = C d t d u + u d t d C 这里的 i d L / d t i\,dL/dt i d L / d t u d C / d t u\,dC/dt u d C / d t
不过，这里要区分两类“可变”。
如果 L L L C C C L n L_{n} L n L n + 1 L_{n+1} L n + 1 C n C_{n} C n C n + 1 C_{n+1} C n + 1
但如果把 L L L C C C L = L ( i ) L=L(i) L = L ( i ) C = C ( u ) C=C(u) C = C ( u ) ψ = L ( i ) i \psi=L(i)i ψ = L ( i ) i q = C ( u ) u q=C(u)u q = C ( u ) u
所以，磁链、电荷建模并不自动意味着“非线性”。关键在于 L L L C C C
什么是非线性方程
前面几篇文章里，我们大多数时候都在处理线性电路。KCL 和 KVL 本身是线性的：节点处电流相加为零，回路中电压相加为零。这类方程的系数通常可以装配成一个矩阵，例如：
G u = i
G\mathbf u=\mathbf i
G u = i 如果元件关系也是线性的，比如：
i = G u
i=Gu
i = G u 或者离散后的电感、电容可以写成：
i n + 1 = G u n + 1 + I h i s
i_{n+1}=G u_{n+1}+I_{\mathrm{his}}
i n + 1 = G u n + 1 + I his 那么整个网络方程仍然是线性代数方程。求解器只需要组装矩阵、求解矩阵方程即可。
非线性通常不是 KCL、KVL 自己产生的，而是元件的端口关系产生的。比如二极管的电流电压关系：
i D = I S ( e u D n V T − 1 )
i_D=I_S\left(e^{\frac{u_D}{nV_T}}-1\right)
i D = I S ( e n V T u D − 1 ) 这不是 i = G u i=Gu i = G u
ψ = ψ ( i )
\psi=\psi(i)
ψ = ψ ( i ) 此时电感量可以理解为工作点相关的增量量：
L i n c = d ψ d i
L_{\mathrm{inc}}=\frac{d\psi}{di}
L inc = d i d ψ 再比如电压相关电容：
q = C ( u ) u
q=C(u)u
q = C ( u ) u 或者电流相关电感：
ψ = L ( i ) i
\psi=L(i)i
ψ = L ( i ) i 这些关系里，未知量和未知量相乘，或者未知量出现在指数、查表曲线、饱和函数里。离散以后，当前拍未知量不会只以一次项出现，因此不能直接整理成一个固定的 G G G
这类方程也不能简单使用 s → z s \rightarrow z s → z s → z s \rightarrow z s → z
工程上处理非线性方程，通常不是直接找一个全局的 s → z s \rightarrow z s → z
F ( x n + 1 ) = 0
F(\mathbf x_{n+1})=0
F ( x n + 1 ) = 0 然后在每个时间步内求这个方程。常见做法包括：
用 Newton-Raphson 方法在当前工作点线性化，反复更新雅可比矩阵和残差； 用割线法、固定点迭代或预测校正方法做近似迭代； 在 EMT 程序里，把非线性元件在某个迭代工作点局部线性化成增量导纳和等效注入源，再和网络方程一起求解。 所以，后面讨论 d ( x y ) d t \frac{d(xy)}{dt} d t d ( x y ) x x x y y y x x x y y y
从 Dxy/Dt 到 Theta 离散
现在考虑更一般的乘积状态：
z = x y
z=x y
z = x y 连续时间下有熟悉的乘积微分：
d ( x y ) d t = x d y d t + y d x d t
\frac{d(xy)}{dt}=x\frac{dy}{dt}+y\frac{dx}{dt}
d t d ( x y ) = x d t d y + y d t d x 但在离散时间里，问题会更细。因为从 n n n n + 1 n+1 n + 1
x n + 1 y n + 1 − x n y n Δ t = Δ ( x y ) Δ t
\frac{x_{n+1}y_{n+1}-x_n y_n}{\Delta t}
=\frac{\Delta(xy)}{\Delta t}
Δ t x n + 1 y n + 1 − x n y n = Δ t Δ ( x y ) 为了把它拆成“一个量乘以另一个量的增量”，可以有不止一种写法。
第一种拆法是加减 x n + 1 y n x_{n+1}y_n x n + 1 y n
x n + 1 y n + 1 − x n y n = x n + 1 y n + 1 − x n + 1 y n + x n + 1 y n − x n y n
x_{n+1}y_{n+1}-x_n y_n
=x_{n+1}y_{n+1}-x_{n+1}y_n+x_{n+1}y_n-x_n y_n
x n + 1 y n + 1 − x n y n = x n + 1 y n + 1 − x n + 1 y n + x n + 1 y n − x n y n 也就是：
x n + 1 y n + 1 − x n y n = x n + 1 ( y n + 1 − y n ) + y n ( x n + 1 − x n )
x_{n+1}y_{n+1}-x_n y_n
=x_{n+1}(y_{n+1}-y_n)+y_n(x_{n+1}-x_n)
x n + 1 y n + 1 − x n y n = x n + 1 ( y n + 1 − y n ) + y n ( x n + 1 − x n ) 因此：
x n + 1 y n + 1 − x n y n Δ t = x n + 1 y n + 1 − y n Δ t + y n x n + 1 − x n Δ t
\frac{x_{n+1}y_{n+1}-x_n y_n}{\Delta t}
=x_{n+1}\frac{y_{n+1}-y_n}{\Delta t}
+y_n\frac{x_{n+1}-x_n}{\Delta t}
Δ t x n + 1 y n + 1 − x n y n = x n + 1 Δ t y n + 1 − y n + y n Δ t x n + 1 − x n 第二种拆法是加减 x n y n + 1 x_n y_{n+1} x n y n + 1
x n + 1 y n + 1 − x n y n = x n + 1 y n + 1 − x n y n + 1 + x n y n + 1 − x n y n
x_{n+1}y_{n+1}-x_n y_n
=x_{n+1}y_{n+1}-x_n y_{n+1}+x_n y_{n+1}-x_n y_n
x n + 1 y n + 1 − x n y n = x n + 1 y n + 1 − x n y n + 1 + x n y n + 1 − x n y n 也就是：
x n + 1 y n + 1 − x n y n = x n ( y n + 1 − y n ) + y n + 1 ( x n + 1 − x n )
x_{n+1}y_{n+1}-x_n y_n
=x_n(y_{n+1}-y_n)+y_{n+1}(x_{n+1}-x_n)
x n + 1 y n + 1 − x n y n = x n ( y n + 1 − y n ) + y n + 1 ( x n + 1 − x n ) 因此：
x n + 1 y n + 1 − x n y n Δ t = x n y n + 1 − y n Δ t + y n + 1 x n + 1 − x n Δ t
\frac{x_{n+1}y_{n+1}-x_n y_n}{\Delta t}
=x_n\frac{y_{n+1}-y_n}{\Delta t}
+y_{n+1}\frac{x_{n+1}-x_n}{\Delta t}
Δ t x n + 1 y n + 1 − x n y n = x n Δ t y n + 1 − y n + y n + 1 Δ t x n + 1 − x n 这两种写法都不是近似，而是严格的代数恒等式。区别只是乘积中哪个因子取 n n n n + 1 n+1 n + 1
为了和 θ \theta θ
x n + θ = θ x n + 1 + ( 1 − θ ) x n
x_{n+\theta}=\theta x_{n+1}+(1-\theta)x_n
x n + θ = θ x n + 1 + ( 1 − θ ) x n y n + 1 − θ = ( 1 − θ ) y n + 1 + θ y n
y_{n+1-\theta}=(1-\theta)y_{n+1}+\theta y_n
y n + 1 − θ = ( 1 − θ ) y n + 1 + θ y n 直观上，这里确实有一种“平均”的思想：不是把乘积微分里的两个因子都粗暴地取在 n n n n + 1 n+1 n + 1 θ \theta θ θ = 0.5 \theta=0.5 θ = 0.5 θ \theta θ 1 1 1 n + 1 n+1 n + 1
那么可以得到严格恒等式：
x n + 1 y n + 1 − x n y n = x n + θ ( y n + 1 − y n ) + y n + 1 − θ ( x n + 1 − x n )
x_{n+1}y_{n+1}-x_n y_n
=x_{n+\theta}(y_{n+1}-y_n)
+y_{n+1-\theta}(x_{n+1}-x_n)
x n + 1 y n + 1 − x n y n = x n + θ ( y n + 1 − y n ) + y n + 1 − θ ( x n + 1 − x n ) 两边除以 Δ t \Delta t Δ t
x n + 1 y n + 1 − x n y n Δ t = x n + θ y n + 1 − y n Δ t + y n + 1 − θ x n + 1 − x n Δ t
\frac{x_{n+1}y_{n+1}-x_n y_n}{\Delta t}
=x_{n+\theta}\frac{y_{n+1}-y_n}{\Delta t}
+y_{n+1-\theta}\frac{x_{n+1}-x_n}{\Delta t}
Δ t x n + 1 y n + 1 − x n y n = x n + θ Δ t y n + 1 − y n + y n + 1 − θ Δ t x n + 1 − x n 这个式子很重要：它看起来像连续时间的乘积微分离散形式，但到目前为止仍然是精确代数恒等式，还没有引入数值近似。
把它放回电感。令：
x = L , y = i , x y = ψ
x=L,\qquad y=i,\qquad xy=\psi
x = L , y = i , x y = ψ 则：
ψ n + 1 − ψ n Δ t = L n + θ i n + 1 − i n Δ t + i n + 1 − θ L n + 1 − L n Δ t
\frac{\psi_{n+1}-\psi_n}{\Delta t}
=L_{n+\theta}\frac{i_{n+1}-i_n}{\Delta t}
+i_{n+1-\theta}\frac{L_{n+1}-L_n}{\Delta t}
Δ t ψ n + 1 − ψ n = L n + θ Δ t i n + 1 − i n + i n + 1 − θ Δ t L n + 1 − L n 而电感端口满足：
u n + θ = ψ n + 1 − ψ n Δ t
u_{n+\theta}=\frac{\psi_{n+1}-\psi_n}{\Delta t}
u n + θ = Δ t ψ n + 1 − ψ n 所以：
u n + θ = L n + θ i n + 1 − i n Δ t + i n + 1 − θ L n + 1 − L n Δ t
u_{n+\theta}
=L_{n+\theta}\frac{i_{n+1}-i_n}{\Delta t}
+i_{n+1-\theta}\frac{L_{n+1}-L_n}{\Delta t}
u n + θ = L n + θ Δ t i n + 1 − i n + i n + 1 − θ Δ t L n + 1 − L n 其中：
u n + θ = ( 1 − θ ) u n + θ u n + 1
u_{n+\theta}=(1-\theta)u_n+\theta u_{n+1}
u n + θ = ( 1 − θ ) u n + θ u n + 1 这个式子就是把非线性乘积状态 ψ = L i \psi=L i ψ = L i θ \theta θ
类似地，对电容令：
x = C , y = u , x y = q
x=C,\qquad y=u,\qquad xy=q
x = C , y = u , x y = q 可得：
i n + θ = C n + θ u n + 1 − u n Δ t + u n + 1 − θ C n + 1 − C n Δ t
i_{n+\theta}
=C_{n+\theta}\frac{u_{n+1}-u_n}{\Delta t}
+u_{n+1-\theta}\frac{C_{n+1}-C_n}{\Delta t}
i n + θ = C n + θ Δ t u n + 1 − u n + u n + 1 − θ Δ t C n + 1 − C n 其中：
i n + θ = ( 1 − θ ) i n + θ i n + 1
i_{n+\theta}=(1-\theta)i_n+\theta i_{n+1}
i n + θ = ( 1 − θ ) i n + θ i n + 1 如果 L n + 1 L_{n+1} L n + 1 C n + 1 C_{n+1} C n + 1 L n + 1 L_{n+1} L n + 1 C n + 1 C_{n+1} C n + 1 L n + θ i n + 1 L_{n+\theta}i_{n+1} L n + θ i n + 1 C n + θ u n + 1 C_{n+\theta}u_{n+1} C n + θ u n + 1
这就是可变 LC 建模里最需要先说清楚的边界：我们到底是在离散一个已知参数调度下的线性元件，还是在求解一个参数也随当前未知状态变化的非线性元件。
磁链和电荷写法的 Dommel 形式
有了前面的乘积差分公式，就可以直接写出可变电感、电容进入节点方程时的 companion model。这里先假设 L L L C C C R R R
先看串联 R L RL R L
d d t ( L ( t ) i ( t ) ) = u ( t ) − R i ( t )
\frac{d}{dt}\left(L(t)i(t)\right)=u(t)-R i(t)
d t d ( L ( t ) i ( t ) ) = u ( t ) − R i ( t ) 使用前面的 θ \theta θ
L n + θ ( i n + 1 − i n ) + i n + 1 − θ ( L n + 1 − L n ) = Δ t ( u n + θ − R i n + θ )
L_{n+\theta}(i_{n+1}-i_n)
+i_{n+1-\theta}(L_{n+1}-L_n)
=\Delta t\left(u_{n+\theta}-R i_{n+\theta}\right)
L n + θ ( i n + 1 − i n ) + i n + 1 − θ ( L n + 1 − L n ) = Δ t ( u n + θ − R i n + θ ) 如果 R R R
i n + 1 = L n − R Δ t ( 1 − θ ) L n + 1 + R Δ t θ i n + Δ t L n + 1 + R Δ t θ u n + θ
i_{n+1}
=\frac{L_n-R\Delta t(1-\theta)}
{L_{n+1}+R\Delta t\theta}i_n
+\frac{\Delta t}
{L_{n+1}+R\Delta t\theta}u_{n+\theta}
i n + 1 = L n + 1 + R Δ tθ L n − R Δ t ( 1 − θ ) i n + L n + 1 + R Δ tθ Δ t u n + θ 而：
u n + θ = ( 1 − θ ) u n + θ u n + 1
u_{n+\theta}=(1-\theta)u_n+\theta u_{n+1}
u n + θ = ( 1 − θ ) u n + θ u n + 1 所以它可以写成当前拍端口电压加历史项：
i n + 1 = G R L , ψ u n + 1 + I h i s , R L , ψ
i_{n+1}=G_{RL,\psi}u_{n+1}+I_{\mathrm{his},RL,\psi}
i n + 1 = G R L , ψ u n + 1 + I his , R L , ψ 其中：
G R L , ψ = θ Δ t L n + 1 + R Δ t θ
G_{RL,\psi}
=\frac{\theta\Delta t}
{L_{n+1}+R\Delta t\theta}
G R L , ψ = L n + 1 + R Δ tθ θ Δ t I h i s , R L , ψ = L n − R Δ t ( 1 − θ ) L n + 1 + R Δ t θ i n + ( 1 − θ ) Δ t L n + 1 + R Δ t θ u n
I_{\mathrm{his},RL,\psi}
=\frac{L_n-R\Delta t(1-\theta)}
{L_{n+1}+R\Delta t\theta}i_n
+\frac{(1-\theta)\Delta t}
{L_{n+1}+R\Delta t\theta}u_n
I his , R L , ψ = L n + 1 + R Δ tθ L n − R Δ t ( 1 − θ ) i n + L n + 1 + R Δ tθ ( 1 − θ ) Δ t u n 如果还考虑 R R R
R n + θ = ( 1 − θ ) R n + θ R n + 1
R_{n+\theta}=(1-\theta)R_n+\theta R_{n+1}
R n + θ = ( 1 − θ ) R n + θ R n + 1 则可写成：
i n + 1 = L n − R n + θ Δ t ( 1 − θ ) L n + 1 + R n + θ Δ t θ i n + Δ t L n + 1 + R n + θ Δ t θ u n + θ
i_{n+1}
=\frac{L_n-R_{n+\theta}\Delta t(1-\theta)}
{L_{n+1}+R_{n+\theta}\Delta t\theta}i_n
+\frac{\Delta t}
{L_{n+1}+R_{n+\theta}\Delta t\theta}u_{n+\theta}
i n + 1 = L n + 1 + R n + θ Δ tθ L n − R n + θ Δ t ( 1 − θ ) i n + L n + 1 + R n + θ Δ tθ Δ t u n + θ 再看串联 R C RC RC
d d t ( C ( t ) u C ( t ) ) = i ( t )
\frac{d}{dt}\left(C(t)u_C(t)\right)=i(t)
d t d ( C ( t ) u C ( t ) ) = i ( t ) 和支路电压关系：
u C ( t ) = v ( t ) − R i ( t )
u_C(t)=v(t)-R i(t)
u C ( t ) = v ( t ) − R i ( t ) 出发，其中 v v v R C RC RC θ \theta θ
i n + 1 = C n + 1 C n + 1 R n + 1 + Δ t θ v n + 1 + C n R n − Δ t ( 1 − θ ) C n + 1 R n + 1 + Δ t θ i n − C n C n + 1 R n + 1 + Δ t θ v n
i_{n+1}
=\frac{C_{n+1}}
{C_{n+1}R_{n+1}+\Delta t\theta}v_{n+1}
+\frac{C_nR_n-\Delta t(1-\theta)}
{C_{n+1}R_{n+1}+\Delta t\theta}i_n
-\frac{C_n}
{C_{n+1}R_{n+1}+\Delta t\theta}v_n
i n + 1 = C n + 1 R n + 1 + Δ tθ C n + 1 v n + 1 + C n + 1 R n + 1 + Δ tθ C n R n − Δ t ( 1 − θ ) i n − C n + 1 R n + 1 + Δ tθ C n v n 因此它的 Dommel 形式为：
i n + 1 = G R C , q v n + 1 + I h i s , R C , q
i_{n+1}=G_{RC,q}v_{n+1}+I_{\mathrm{his},RC,q}
i n + 1 = G RC , q v n + 1 + I his , RC , q 其中：
G R C , q = C n + 1 C n + 1 R n + 1 + Δ t θ
G_{RC,q}
=\frac{C_{n+1}}
{C_{n+1}R_{n+1}+\Delta t\theta}
G RC , q = C n + 1 R n + 1 + Δ tθ C n + 1 I h i s , R C , q = C n R n − Δ t ( 1 − θ ) C n + 1 R n + 1 + Δ t θ i n − C n C n + 1 R n + 1 + Δ t θ v n
I_{\mathrm{his},RC,q}
=\frac{C_nR_n-\Delta t(1-\theta)}
{C_{n+1}R_{n+1}+\Delta t\theta}i_n
-\frac{C_n}
{C_{n+1}R_{n+1}+\Delta t\theta}v_n
I his , RC , q = C n + 1 R n + 1 + Δ tθ C n R n − Δ t ( 1 − θ ) i n − C n + 1 R n + 1 + Δ tθ C n v n 这里可以看到，磁链写法的 R L RL R L R C RC RC L n , L n + 1 L_n,L_{n+1} L n , L n + 1 C n , C n + 1 C_n,C_{n+1} C n , C n + 1
磁链和电荷建模的意义
到这里可以回头回答一个更本质的问题：为什么要引入磁链和电荷？它最重要的意义，并不是简单地“允许电压、电流突变”，而是明确区分开关事件中真正保存系统记忆的状态量，以及可以随代数关系发生跳变的端口量。
电容和电感的基本状态方程是：
d q d t = i
\frac{dq}{dt}=i
d t d q = i d ψ d t = u
\frac{d\psi}{dt}=u
d t d ψ = u 在一个开关瞬间，对这两个方程积分：
q + − q − = ∫ t − t + i ( t ) d t
q^+-q^-=\int_{t^-}^{t^+}i(t)\,dt
q + − q − = ∫ t − t + i ( t ) d t ψ + − ψ − = ∫ t − t + u ( t ) d t
\psi^+-\psi^-=\int_{t^-}^{t^+}u(t)\,dt
ψ + − ψ − = ∫ t − t + u ( t ) d t 如果开关瞬间没有冲激电流和冲激电压，那么右侧积分为零，于是：
q + = q −
q^+=q^-
q + = q − ψ + = ψ −
\psi^+=\psi^-
ψ + = ψ − 也就是说，通常真正连续的是电荷和磁链。电压、电流是否连续，还要看参数和本构关系是否改变。
对于固定电容：
q = C v
q=Cv
q = C v 若 q q q C C C
v + = v −
v^+=v^-
v + = v − 所以普通开关事件不会让固定电容电压无条件突变，但电容电流可以突变。对于时变电容，如果电容从 C − C^- C − C + C^+ C +
q + = q −
q^+=q^-
q + = q − 于是：
v + = q − C + = C − C + v −
v^+=\frac{q^-}{C^+}
=\frac{C^-}{C^+}v^-
v + = C + q − = C + C − v − 这说明电荷连续并不意味着电压一定连续。当电容参数突变时，电压会根据新的 q = C v q=Cv q = C v
电感完全对偶。若：
ψ = L i
\psi=Li
ψ = L i 且 L L L L − L^- L − L + L^+ L +
i + = ψ − L + = L − L + i −
i^+=\frac{\psi^-}{L^+}
=\frac{L^-}{L^+}i^-
i + = L + ψ − = L + L − i − 所以，磁链连续时，时变电感的电流也可能突变。
更准确地说，磁链和电荷建模不是人为放宽电压、电流的连续性，而是先保存真正的物理状态，再通过事件后的本构关系计算新的电压或电流。它避免了把连续性错误地施加在不应该连续的变量上。
下面用一个孤立时变电容的小例子说明这个差别。
设一个电容与外部电路隔离，端口电流始终为零：
i t e r m i n a l = 0
i_{\mathrm{terminal}}=0
i terminal = 0 事件前：
C − = 10 μ F , v − = 10 V
C^-=10\,\mu\mathrm{F},\qquad v^-=10\,\mathrm{V}
C − = 10 μ F , v − = 10 V 事件后：
C + = 5 μ F
C^+=5\,\mu\mathrm{F}
C + = 5 μ F 初始电荷为：
q − = C − v − = 10 μ F × 10 V = 100 μ C
q^-=C^-v^-
=10\,\mu\mathrm{F}\times 10\,\mathrm{V}
=100\,\mu\mathrm{C}
q − = C − v − = 10 μ F × 10 V = 100 μ C 如果使用完整的电荷平衡模型：
q n − q n − 1 = h i n
q_n-q_{n-1}=h i_n
q n − q n − 1 = h i n 又因为：
q n = C n v n
q_n=C_n v_n
q n = C n v n 所以：
C n v n − C n − 1 v n − 1 = h i n
C_n v_n-C_{n-1}v_{n-1}=h i_n
C n v n − C n − 1 v n − 1 = h i n 本例中 i n = 0 i_n=0 i n = 0
q n = q n − 1 = 100 μ C
q_n=q_{n-1}=100\,\mu\mathrm{C}
q n = q n − 1 = 100 μ C 事件后电压应为：
v + = q − C + = 100 μ C 5 μ F = 20 V
v^+=\frac{q^-}{C^+}
=\frac{100\,\mu\mathrm{C}}{5\,\mu\mathrm{F}}
=20\,\mathrm{V}
v + = C + q − = 5 μ F 100 μ C = 20 V 也就是说，电容从 10 μ F 10\,\mu\mathrm{F} 10 μ F 5 μ F 5\,\mu\mathrm{F} 5 μ F 100 μ C 100\,\mu\mathrm{C} 100 μ C 10 V 10\,\mathrm{V} 10 V 20 V 20\,\mathrm{V} 20 V
如果错误地只写成：
C n ( v n − v n − 1 ) = h i n
C_n(v_n-v_{n-1})=h i_n
C n ( v n − v n − 1 ) = h i n 那么在 i n = 0 i_n=0 i n = 0
v n = v n − 1
v_n=v_{n-1}
v n = v n − 1 也就是电压仍保持 10 V 10\,\mathrm{V} 10 V
q + = C + v + = 5 μ F × 10 V = 50 μ C
q^+=C^+v^+
=5\,\mu\mathrm{F}\times 10\,\mathrm{V}
=50\,\mu\mathrm{C}
q + = C + v + = 5 μ F × 10 V = 50 μ C 这意味着一个端口电流始终为零的孤立电容，凭空损失了 50 μ C 50\,\mu\mathrm{C} 50 μ C
C d v d t
C\frac{dv}{dt}
C d t d v 却遗漏了完整关系中的参数变化项：
i = d ( C v ) d t = C d v d t + v d C d t
i=\frac{d(Cv)}{dt}
=C\frac{dv}{dt}+v\frac{dC}{dt}
i = d t d ( C v ) = C d t d v + v d t d C 这里比较的不是“所有电压建模”和“电荷建模”的优劣。即使网络未知量仍然选择电压，只要离散方程使用完整的 Δ ( C v ) \Delta(Cv) Δ ( C v ) q = C v q=Cv q = C v C Δ v C\Delta v C Δ v
这个例子还有一个能量问题。电容储能为：
W = 1 2 C v 2 = q 2 2 C
W=\frac12Cv^2=\frac{q^2}{2C}
W = 2 1 C v 2 = 2 C q 2 事件前：
W − = 1 2 ( 10 μ F ) ( 10 V ) 2 = 0.5 m J
W^-=\frac12(10\,\mu\mathrm{F})(10\,\mathrm{V})^2
=0.5\,\mathrm{mJ}
W − = 2 1 ( 10 μ F ) ( 10 V ) 2 = 0.5 mJ 事件后：
W + = 1 2 ( 5 μ F ) ( 20 V ) 2 = 1 m J
W^+=\frac12(5\,\mu\mathrm{F})(20\,\mathrm{V})^2
=1\,\mathrm{mJ}
W + = 2 1 ( 5 μ F ) ( 20 V ) 2 = 1 mJ 增加的能量不是由电气端口提供的，因为端口电流为零；它来自改变电容参数的外部机械端口、控制端口或其他参数端口。这一点也提醒我们：时变参数元件不是普通固定元件，它隐含了一个可以与电磁系统交换能量的参数通道。
下面这张波形图把两种写法的结果放在一起。电荷建模会保持完整状态 q = C v q=Cv q = C v C C C
孤立时变电容：电荷建模保持 q = C v，固定 C 写法会错误地丢失电荷。 下面这段代码可以用来复现实验中的电压和电荷变化。它的重点不是说“选择 q q q v v v q = C v q=Cv q = C v
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
import numpy as np
t_event = 1.0
t = np . linspace ( 0.0 , 2.0 , 9 )
# C changes from 10 uF to 5 uF at t = 1 s.
C = np . where ( t < t_event , 10e-6 , 5e-6 )
v0 = 10.0
q0 = C [ 0 ] * v0
# Complete charge-balance model: q[n] - q[n-1] = h*i[n].
# The isolated capacitor has zero terminal current, so q is constant.
q_balance = np . full_like ( t , q0 )
v_balance = q_balance / C
# Incomplete model: C[n]*(v[n] - v[n-1]) = h*i[n].
# With i = 0, this incorrectly keeps voltage constant.
v_incomplete = np . full_like ( t , v0 )
q_incomplete = C * v_incomplete
print ( "charge balance voltage:" , v_balance )
print ( "incomplete model voltage:" , v_incomplete )
print ( "charge balance charge/uC:" , q_balance * 1e6 )
print ( "incomplete model charge/uC:" , q_incomplete * 1e6 )
支付宝
微信