梯形积分为什么会数值震荡

RLC侠收录于

2026-06-17 2026-07-01 约 6000 字预计阅读 12 分钟 - 次阅读

前面几篇文章一直在给梯形积分铺路：从 $1/s$ 积分算子，到 $s$ 到 $z$ 离散，再到电感、电容的 Dommel 伴随模型。梯形积分确实很好用，它有二阶精度，对线性电路也非常适合整理成“等效导纳 + 历史源”的形式。

但工程里还有另一句话：梯形法容易在开关突变后出现数值震荡。

这句话听起来有点矛盾。梯形法不是稳定性很好吗？为什么又会震荡？这一篇就专门拆这个问题。

梯形法好在哪里

先看最普通的积分问题：

\frac{\mathrm{d}y}{\mathrm{d}t}=f(t,y)

梯形法写成：

y_k=y_{k-1}+\frac{T}{2}(f_k+f_{k-1})

其中 $T$ 是仿真步长。它的直觉很简单：一步内的积分面积不用左端点，也不用右端点，而是用两端平均值形成的梯形面积。

对平滑波形来说，这个做法很漂亮：

精度比一阶欧拉高；
对很多线性系统稳定性很好；
在电磁暂态仿真中，电感、电容可以自然整理成 Dommel 形式。

问题恰恰也藏在这里：梯形法太“公平”了，它把上一步和当前步各取一半。对平滑波形，这是优点；对开关突变后的高频误差，它就可能缺少足够的数值阻尼。

数值震荡从哪里来

考虑一个最简单的衰减方程：

\frac{\mathrm{d}y}{\mathrm{d}t}=-ay

其中 $a>0$ 。精确解会单调衰减：

y(t)=y(0)e^{-at}

对它使用梯形法：

y_k=y_{k-1}+\frac{T}{2}(-ay_k-ay_{k-1})

整理得到：

y_k=\frac{1-\frac{aT}{2}}{1+\frac{aT}{2}}y_{k-1}

令 $q=aT$ ，放大因子为：

A_{\mathrm{trap}}=\frac{1-\frac{q}{2}}{1+\frac{q}{2}}

如果 $q$ 很小， $A_{\mathrm{trap}}$ 是一个接近 1 的正数，数值解平滑衰减。

但如果 $q>2$ ，分子变成负数， $A_{\mathrm{trap}}$ 也变成负数。于是 $y_k$ 会一步正、一步负地交替变化。更极端地，当 $q$ 非常大时：

A_{\mathrm{trap}}\rightarrow -1

这就意味着：真实系统可能已经在一个步长内衰减得几乎没有了，但梯形离散结果却可能留下一个接近 $(-1)^k$ 的交替分量。

这就是很多电磁暂态仿真里说的 numerical oscillation、numerical chatter，或者中文里常说的数值震荡。

频率和步长有什么关系

一步正、一步负的序列可以写成：

1,-1,1,-1,\cdots

也就是：

(-1)^k

在离散系统里，这对应每两个采样点完成一个周期。若步长为 $T$ ，震荡周期就是：

T_{\mathrm{osc}}=2T

对应频率为：

f_{\mathrm{osc}}=\frac{1}{2T}

这正是采样系统里的 Nyquist 频率。也就是说，梯形法在开关突变后最容易留下来的那种交替震荡，往往是离散系统能表示的最高频附近的数值分量。

这里容易说错。不能说“震荡频率是步长的 2 倍频”。更准确的说法是：这个交替分量两拍完成一个周期，周期为 $2T$ ，频率为 $1/(2T)$ 。如果把采样频率记为 $f_s=1/T$ ，那么：

f_{\mathrm{osc}}=\frac{f_s}{2}

也就是采样频率的一半，而不是 $2/T$ 。

所以步长越小，这个震荡频率越高：

$T=50\,\mu s$ 时，交替震荡频率约为 $10\,kHz$ ；
$T=10\,\mu s$ 时，交替震荡频率约为 $50\,kHz$ 。

但要注意，减小步长并不等于自动消除震荡。只要开关突变使局部时间常数仍然远小于步长，梯形法仍然可能把高频误差留在计算结果里。

一个开关突变的 RC 例子

现在看一个很小的电路例子。电容初始电压为 $1\,V$ 。在 $t_s$ 之前，开关断开，电容电压保持不变；在 $t_s$ 之后，开关闭合，电容通过一个很小的电阻放电。

闭合后的方程是：

\frac{\mathrm{d}u_C}{\mathrm{d}t}=-\frac{1}{RC}u_C

令时间常数：

\tau=RC

如果仿真步长 $T=50\,\mu s$ ，而开关闭合后的 $\tau=2\,\mu s$ ，则：

q=\frac{T}{\tau}=25

这时梯形法的放大因子为：

A_{\mathrm{trap}}=\frac{1-\frac{25}{2}}{1+\frac{25}{2}}\approx -0.852

也就是说，电容电压每一步都会翻符号，并且幅值只慢慢衰减。这不是物理振荡，而是数值离散造成的交替分量。

下面这段 Python 代码复现了这个现象，同时给出后向欧拉的结果：

 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
import numpy as np
import matplotlib.pyplot as plt

T = 50e-6
t_end = 2.5e-3
t_switch = 0.5e-3
tau_closed = 2e-6

t = np.arange(0.0, t_end + T, T)
v_trap = np.ones_like(t)
v_be = np.ones_like(t)

q = T / tau_closed
a_trap = (1.0 - q / 2.0) / (1.0 + q / 2.0)
a_be = 1.0 / (1.0 + q)

for k in range(1, len(t)):
    if t[k] >= t_switch:
        v_trap[k] = a_trap * v_trap[k - 1]
        v_be[k] = a_be * v_be[k - 1]
    else:
        v_trap[k] = v_trap[k - 1]
        v_be[k] = v_be[k - 1]

plt.plot(t * 1e3, v_trap, "o-", label="Trapezoidal")
plt.plot(t * 1e3, v_be, "s-", label="Backward Euler")
plt.axvline(t_switch * 1e3, color="k", linestyle="--", linewidth=1)
plt.xlabel("Time / ms")
plt.ylabel("Capacitor voltage / V")
plt.legend()
plt.grid(True)
plt.show()

图里可以看到：开关闭合后，后向欧拉很快把电压压到接近 0；梯形法却留下了明显的一步正、一步负交替震荡。

这个例子故意选了比较极端的参数，是为了把现象放大。真实电磁暂态程序里，开关、二极管、电力电子器件、避雷器、理想电感电容组合，都可能制造类似的局部突变。

这里还可以顺手澄清一个问题：并不是只有状态突变时，梯形法才“拥有”数值震荡模式。更准确地说，梯形法本来就缺少对最高频数值分量的阻尼；只是平滑运行时，这类高频误差很小，不容易被看见。开关突变、拓扑突变或参数突变会突然改变微分方程的局部形态，制造很强的高频误差，其中就可能包含接近 $(-1)^k$ 的交替分量。梯形法又刚好不太衰减这个分量，于是它就残留下来，表现为一步正、一步负的数值震荡。

所以，突变事件不是数值震荡模式存在的唯一原因，而是最容易把这个模式激发出来的“敲击”。

从初值不协调看数值震荡

还有一种常见解释，是从“初值不协调”或者“历史项不协调”的角度看这个问题。

开关动作前后，电路方程往往不是同一个方程。开关断开时，电容可能几乎没有放电通路；开关闭合后，电容突然接上一个很小的电阻，新的局部时间常数可能非常小。物理上，电容电压不能突变，电感电流也不能突变：

u_C^+=u_C^-

i_L^+=i_L^-

但是，电路的代数约束、支路电流、导数关系已经变了。也就是说，状态量本身连续，不代表它们和新拓扑下的导数、支路电流、历史项立刻协调。

梯形法尤其容易暴露这个问题，因为它不是只看当前点，而是把上一拍和当前拍平均起来：

y_k=y_{k-1}+\frac{T}{2}(f_k+f_{k-1})

这里的 $f_{k-1}$ 来自开关前的旧电路，而 $f_k$ 来自开关后的新电路。如果开关瞬间直接沿用旧历史项，就相当于让积分器拿着旧电路的记忆去接新电路的方程。这个交接处只要留下一个小误差，就可能被梯形法的交替模态接住。

仍然看测试方程：

y'=-ay

梯形法的误差递推因子为：

A_{\mathrm{trap}}=\frac{1-\frac{aT}{2}}{1+\frac{aT}{2}}

如果开关后的局部时间常数很小，也就是 $aT$ 很大，则：

A_{\mathrm{trap}}\approx -1

于是初值或历史项中残留的不协调误差会近似满足：

e_k\approx -e_{k-1}

这就是一拍正、一拍负的数值震荡。换句话说，开关突变并不是凭空制造了一个新的物理振荡，而是让离散初值没有完全落在新系统的正确解轨道上；梯形法又缺少对这个寄生交替模态的阻尼，于是它被保留下来。

从这个角度看，数值震荡可以理解为：切换瞬间的状态和历史项没有与新拓扑重新协调，激发了梯形法的寄生模态。 这和前面的频率解释并不矛盾。频率解释说它表现为接近 Nyquist 频率的交替分量；初值解释说这个交替分量是怎样在切换瞬间被激发出来的。

后向欧拉为什么更稳

同样对：

\frac{\mathrm{d}y}{\mathrm{d}t}=-ay

使用后向欧拉：

y_k=y_{k-1}+T(-ay_k)

整理得：

y_k=\frac{1}{1+aT}y_{k-1}

也就是：

A_{\mathrm{BE}}=\frac{1}{1+q}

无论 $q$ 多大， $A_{\mathrm{BE}}$ 都是正数，而且 $q$ 越大， $A_{\mathrm{BE}}$ 越接近 0。对于突变后的刚性衰减过程，后向欧拉会表现出很强的数值阻尼。

这也是为什么一些 EMTP 类程序会在检测到开关不连续或数值震荡时，临时使用后向欧拉、CDA（Critical Damping Adjustment，临界阻尼调整）或类似策略，然后再切回梯形法。

当然，后向欧拉不是白捡的好处。它只有一阶精度，数值阻尼也可能把真实的高频暂态一起抹掉。所以它更像是一块“刹车片”，而不是所有场景下的最佳积分器。

积分因子：0.5 还是 1

现在问题来了：既然梯形法精度好但阻尼弱，后向欧拉阻尼强但精度低，那有没有中间方案？

有。这就是常说的 theta method，或者 $\theta$ 积分。

把积分公式写成统一形式：

y_k=y_{k-1}+T\left[(1-\theta)f_{k-1}+\theta f_k\right]

当：

$\theta=0$ ：显式欧拉；
$\theta=0.5$ ：梯形法；
$\theta=1$ ：后向欧拉；
$0.5<\theta<1$ ：带数值阻尼的隐式积分。

对测试方程 $y'=-ay$ ，它的放大因子为：

A_{\theta}=\frac{1-(1-\theta)q}{1+\theta q}

当 $q$ 很大时：

A_{\theta}\rightarrow -\frac{1-\theta}{\theta}

这个式子很有意思：

$\theta=0.5$ 时，高频极限是 $-1$ ，几乎没有高频阻尼；
$\theta=1$ 时，高频极限是 $0$ ，高频分量被强烈压掉；
$0.5<\theta<1$ 时，高频分量会被不同程度地阻尼。

所以，“积分因子到底取 0.5 还是 1”并不是一个纯数学洁癖问题，而是一个工程权衡：

想保留二阶精度和较小数值阻尼，就靠近 $0.5$ ；
想抑制开关突变后的高频 chatter，就靠近 $1$ ；
想折中，就选 $0.5$ 和 $1$ 之间的某个值。

阻尼积分、可调积分、alpha/Beta 积分

在不同资料和软件里，这类思想会有不同名字：

theta method；
generalized trapezoidal method；
阻尼积分；
可调积分；
alpha 积分或 alpha-beta 积分。

名字不完全一样，但背后的思路很接近：不要把积分面积永远固定成左右各一半，而是允许当前端点占更多权重，从而给高频数值分量加一点阻尼。

如果沿用上面的 $\theta$ 写法，梯形法就是 $\theta=0.5$ ，后向欧拉就是 $\theta=1$ 。阻尼积分可以理解为：

0.5<\theta<1

有些工程书或程序会使用 $\alpha$ 、 $\beta$ 来描述类似的权重。读到这些名字时，不要先被术语吓住，先去看它到底把 $f_{k-1}$ 和 $f_k$ 的权重怎么分配。只要看清权重，就能判断它更接近梯形法，还是更接近后向欧拉。

工程上怎么理解

这篇文章不是要否定梯形法。恰恰相反，梯形法在电磁暂态仿真中非常重要。它的问题不是“普通情况下不稳定”，而是：

开关动作会制造不连续；
不连续会激发很高频的数值分量；
梯形法对最高频附近的分量阻尼很弱；
于是可能出现一步正、一步负的数值震荡；
后向欧拉或阻尼积分可以压制这类分量，但会牺牲一部分精度或真实高频响应。

所以，积分方法不是越高阶越好，也不是阻尼越强越好。真正的工程问题是：当前仿真更怕什么？

如果更怕相位误差和过度阻尼，梯形法很合适。
如果更怕开关突变后的 chatter，后向欧拉、CDA、theta 积分或可调积分就值得考虑。
如果系统里有大量电力电子开关，光靠固定步长梯形法往往不够，还需要开关时刻插值、事件定位、状态重初始化等配套策略。