上一篇文章我们从常微分方程的角度，粗略讲了显式欧拉、隐式欧拉和梯形积分法。那一篇更多是在数学意义上回答一个问题：连续时间里的微分方程，为什么可以被一步一步算出来？
但是真正进入仿真软件以后，读者看到的通常不是一行微分方程，而是一个个模块。比如在 Simulink 里，最常见、也最容易被忽略的模块，就是积分算子。
Simulink 积分算子示意
如果输入是x(t)，输出是y(t)，那么这个模块表达的连续时间关系就是：
y(t)=∫0tx(τ)dτ+y(0)写成微分形式就是：
dtdy=x(t)写成拉普拉斯域，就是大家熟悉的：
Y(s)=s1X(s)所以，一个看起来很简单的1/s，其实已经把“状态”“初值”“积分”“微分方程求解”这些概念全都装进去了。理解这个模块如何在数字仿真环境里被计算，是理解后面 Dommel 算法、电感电容等效电路、状态更新公式的第一步。
仿真软件到底在算什么
先不要急着想 Simulink 内部有多少复杂求解器。只看最简单的情况：输入x，经过积分器1/s，输出y。
连续时间里，我们可以说任意时刻都有：
y˙(t)=x(t)但是数字计算机不可能真的在“任意时刻”计算。它只能在一系列离散时间点上计算：
t0, t1, t2, ⋯, tk, tk+1如果固定步长为T，那么：
tk=kT于是连续信号x(t)和y(t)，在计算机里就变成了两个序列：
x[0],x[1],x[2],⋯y[0],y[1],y[2],⋯这时候，积分器的任务就从“求连续积分”变成了“根据过去和当前的采样值，更新下一步状态”。这也是仿真软件做数值积分的核心。
积分方法的简单回顾
我们上一篇已经介绍了如何从根源推导数值积分方法比如梯形法则，这里简单回顾一下，连续积分可以理解成面积：
y(tk)−y(tk−1)=∫tk−1tkx(t)dt现在问题就变成了：右边这个面积怎么算？
最朴素的办法，是假设这一小段时间里x(t)近似等于左端点x[k−1]，那么：
∫tk−1tkx(t)dt≈Tx[k−1]于是得到：
y[k]=y[k−1]+Tx[k−1]这就是显式欧拉积分。它很直观：拿上一步的斜率，往前走一步。
如果改用右端点x[k]近似这一段面积：
∫tk−1tkx(t)dt≈Tx[k]就得到：
y[k]=y[k−1]+Tx[k]这就是隐式欧拉积分。
如果用梯形面积近似：
∫tk−1tkx(t)dt≈2T(x[k]+x[k−1])就得到：
y[k]=y[k−1]+2T(x[k]+x[k−1])这就是梯形积分。它比单纯取左端点或右端点更平衡，也正是后面电磁暂态仿真里非常重要的一条路。
从s域到z域
现在再回到最开始那个积分算子。连续系统里，输入x(t)和输出y(t)满足：
Y(s)=s1X(s)也就是说，连续域里的积分器，本质上就是一个s1算子。
如果我们要把这个连续积分器搬到数字仿真环境里，一个很直接的想法是：能不能把s换成某个关于z的表达式？
比如，(先不解释它从哪来) 我们直接采用梯形法常用的替换：
s≈T21+z−11−z−1那么连续积分器里的s1就会变成：
s1≈2T1−z−11+z−1把它代回最开始的输入输出关系：
Y(s)=s1X(s)在z域里就得到：
Y(z)=2T1−z−11+z−1X(z)两边同乘1−z−1：
(1−z−1)Y(z)=2T(1+z−1)X(z)展开：
Y(z)−z−1Y(z)=2TX(z)+2Tz−1X(z)再回到离散时间序列，就是：
y[k]−y[k−1]=2T(x[k]+x[k−1])整理一下：
y[k]=y[k−1]+2T(x[k]+x[k−1])这不就是前面刚刚看到的梯形积分吗？
换句话说，如果我们事先知道：
s≈T21+z−11−z−1那么从连续域的s1积分器，到离散域的梯形积分递推，就可以很快推出来。
那问题来了：为什么s可以这样近似？
现在反过来看。梯形积分的出发点是把一小段面积近似成梯形：
∫tk−1tkx(t)dt≈2T(x[k]+x[k−1])因为y是x的积分，所以有：
y[k]−y[k−1]=2T(x[k]+x[k−1])对这个递推式做z变换：
(1−z−1)Y(z)=2T(1+z−1)X(z)于是：
X(z)Y(z)=2T1−z−11+z−1而连续域里：
X(s)Y(s)=s1所以：
s1≈2T1−z−11+z−1反过来就得到：
s≈T21+z−11−z−1这就是梯形法，也就是 Tustin 变换，对应的s→z关系。
所以，所谓s到z的转换，并不是凭空规定一个代换公式。它背后真正发生的事情是：某一种积分方法，先给出了s1在z域里的离散实现；然后我们再把这个离散实现整理成s的替换关系。
如果事先把常见积分方法对应的s和z关系总结出来，那么以后看到连续域里的传递函数、积分算子或微分算子，就可以更快地把它搬到离散域里。这就是下面这张表的意义。
常见 S 到 Z 转换表
下面这张表先采用固定步长T，并采用常见的递推约定。不同软件、不同模块、不同输入输出对齐方式，可能会在是否带一个采样延迟上有细微差别；后续做具体电路等效时，要以实际离散方程为准。
| 方法 | 积分递推形式 | s1 的离散近似 | s 的替换关系 | 常见特点 |
|---|
| 显式欧拉 Forward Euler | y[k]=y[k−1]+Tx[k−1] | 1−z−1Tz−1 | Tz−1 | 形式简单，使用上一拍输入，稳定性较弱 |
| 隐式欧拉 Backward Euler | y[k]=y[k−1]+Tx[k] | 1−z−1T | T1−z−1 | 数值阻尼较强，稳定性好，但可能需要联立求解 |
| 梯形法 / Tustin | y[k]=y[k−1]+2T(x[k]+x[k−1]) | 2T1−z−11+z−1 | T21+z−11−z−1 | 二阶精度，频域映射性质好，电磁暂态中非常常见 |
| 二阶 Gear / BDF2 | 2T3y[k]−4y[k−1]+y[k−2]=x[k] | 3−4z−1+z−22T | 2T3−4z−1+z−2 | 多步隐式方法，数值阻尼更强，常用于抑制振荡 |
这里要特别注意：表里的s→z并不是凭空规定出来的，而是从不同的积分近似方法推出来的。换句话说，如果你在连续系统里看到一个1/s，你想把它放进数字仿真环境，就必须决定这个1/s到底用哪一种离散积分公式来实现。
为什么这和电磁暂态仿真有关
电磁暂态仿真里最基本的动态元件，是电感和电容：
uL=LdtdiLiC=CdtduC看到这里，可以先停一下。这个形式是不是和前面的积分器非常像？
前面我们讲的是：
dtdy=x(t)写到s域就是：
Y(s)=s1X(s)而电感、电容只是把变量换成了电压和电流。比如电容关系：
iC=CdtduC可以改写成：
dtduC=C1iC再写成积分形式，比如电容：
uC(t)=C1∫iC(t)dt看，这和y=x/s其实是同一件事：输入是电流，输出是电压，中间差了一个C1比例系数和一个积分算子s1。
电感也是类似的：
dtdiL=L1uL所以从数学本质上说，电磁暂态仿真并不是突然冒出一套完全不同的理论。它首先是在求解一组由电感、电容、线路和控制环节组成的常微分方程，或者微分代数方程。只要进入数字仿真环境，就必须把这些连续关系改写成离散递推。
选择显式欧拉、隐式欧拉、梯形法，或者 Gear 方法，本质上就是选择不同的 ODE 数值积分方式，进而得到不同的离散等效模型。
Dommel 算法之所以经典，一个重要原因就是它把这种离散化后的动态元件，整理成了电路工程师熟悉的形式：等效导纳 + 历史电流源。这样，电感、电容、线路等动态元件就可以被放进节点导纳矩阵里统一求解（下一章会讲到）。
这一篇先不展开 RLC Dommel 等效电路。你只需要记住一件事：从1/s积分器到离散递推，再到s→z替换，是连续仿真模型进入数字计算机的入口。后面讲 Dommel 算法时，很多看起来突然出现的“历史项”，其实都来自这里。
小结
这一篇只做一件事：把 Simulink 里看似简单的积分算子拆开。
- 连续时间里，积分算子表示Y(s)=s1X(s)。
- 数字仿真环境里，计算机只能在采样点上更新状态。
- 不同的面积近似方式，会得到不同的离散递推公式。
- 每一种递推公式，都对应一种s→z替换关系。
- 电磁暂态仿真中的电感、电容离散化，本质上也是这个问题。
读到这里，读者可能会发现，所谓s域、z域、数值积分、仿真软件求解器，其实不是几套互不相干的知识。它们都在回答同一个工程问题：连续世界里的微分方程，如何被数字计算机一步一步算出来。
参考资料
- 华冬英、李祥贵，《微分方程的数值解法与程序实现》。
- Neville Watson and Jos Arrillaga 著，陈贺、白宏、项祖涛译，《电力系统电磁暂态仿真 Power Systems Electromagnetic Transients Simulation》。
支付宝
微信