数值仿真入门

Dommel EMTP 算法是求解电磁暂态问题最被广泛采用的办法，但是讲解 Dommel 算法之前，其实需要对微分方程的数值解法有一定的基础。不然上来就学习 RLC 的离散化，可能对初学者有点劝退，难以理解。而且电网也不只有 RLC 元件，还有非线性元件、需要用状态空间描述的外部控制环路等等。所以，虽然是做电气领域的仿真算法，但是万变不离其宗，本质上还是数学问题。

但是微分方程的数值解法，又是一门深奥的学科，笔者在入门的时候也走过不少弯路，因此萌生了写这篇文章的想法。我们最终的目的是深刻理解并灵活应用 Dommel 算法去求解电磁暂态问题，所以本文不会深入“常微分方程的数值解”这片深海，只是从电气工程师，仿真工程师角度，带大家入门，希望能帮助初学者脱离苦海。

这里推荐两边书作为本文的延伸阅读，首先是华冬英，李祥贵的《微分方程的数值解法》，初学者建议读第一、二章即可，浅尝辄止。然后去读 Neville Watson and Jos Arrillagaz 著，陈贺、白宏、项祖涛译的《电力系统电磁暂态仿真 Power Systems Electromagnetic Transients Simulation》。这样会对 Dommel 算法有一个比较深刻的理解。

本文不会深入 RLC 电路基于 Dommel 算法的等效电路的说明，仅讨论一般常微分方程 (ODE) 的数值解的问题。

数学上的 ODE 数值解法

其实理想情况下的完美解是通过初等方法求出的解析解，比如$y^{\prime}=\sin(x)$, 我们很容易可以得到它的解析解就是$y=-\cos{(x)}+C,C\in\mathbb{R}$. 但是，如果是$y^{\prime}=\sin(y^2)$或者$y^{\prime}=\sin(xy)$就完全没有办法，或者很困难求得解析解。

而经常做仿真的读者也知道，一个完整的物理系统往往是很复杂的，不用说推导解析解，就是写出状态空间形式都很麻烦。而且在有控制反馈的闭环系统中，随手添加一个单位延时、引入滤波、超前、滞后环节都是都很正常的，这都会增加微分变量进而增加系统的复杂程度。所以，严格推导解析解是很困难的，由此也衍生出模型的扫频近似（比如基于 Vector Fitting、零极点匹配的电路有理近似）、特定物理问题的解析解等等方向。

那么在工程上，我们可以通过一些近似方法解出特定条件下的近似值，假设$x$是时间，假设第 0 秒是起始点（工程上时间是负的没有意义），我们人为设定$y$的初始值是$\alpha$，有没有办法可以解出第$n$秒时 y 的近似值呢？这个其实就是常微分方程的初值问题，归纳总结为如下形式：

当$x$不是时间，也就是说不是时域的时候，$x$的范围可以引申为$\quad a<x\leqslant b$, 这个在数学上其实是没有区别的，只不过电气工程上更关心时域问题。$y$的初值可能是电感的电流，或者电容的电压，不过不正确的初值可能会引发其他的问题，这里先不讨论。

好，那么我们下一步怎么处理 (1) 呢？这里介绍容易让读者理解的两个路线，一个是针对$\frac{\mathrm{d}y}{\mathrm{d}x}$，从差商近似的角度；另一个就是，对等式两端进行积分，针对$y(x_{i+1})-y(x_i)=\int_{x_i}^{x_{i+1}}f(x,y(x))\mathrm{d}x$角度入手。

有经验的读者其实可能已经知道欧拉、梯形法则，了解下一步离散化的方法了。但是笔者这里建议先忘记欧拉、梯形法则，我们从源头看看 (1) 到底是怎么解的。

差商近似 —— 由求导的思想引出欧拉方法

要求解$\frac{\mathrm{d}y}{\mathrm{d}x}$，或者$y'$，很容易联想到$y^{\prime}(x_i)\approx\frac{y(x_{i+1})-y(x_i)}{h}$，这个基于图解法很好理解。抛开图解法，它的来源其实是泰勒公式：

从计算机数值求解的角度，我们观察一下泰勒公式，下一个步长的值，$y(x_0+h)$是我们要求解的，$y(x_0)$和$x_0$的初值是人为设定的，步长$h$也是可以设定的，而根据公式 (1)，$y^{\prime}(x_0)$又是可以直接求解的，我们不知道的是高阶导形式$y^{(n)}(x_0)$。

理论上来说，只要有任意$y^{(n)}(x_0)$，我们的泰勒近似就可以达到任意精度。但是实际操作中，高阶导的求解往往非常困难复杂，所以我们牺牲一些精度，在本例中，$y^{\prime\prime}(x_0)$往后我们全都放弃，重新整理泰勒公式可得：

很容易通过上述公式求解出我们要的$y(x_0+h)$。它的一般形式就是：

如果一个数值方法的局部截断误差 (LTE)=$O(h^{p+1})$，则称该数值方法是$p$阶的。所以，这个方法的精度就是一阶的。至此，我们通过差商近似的方法，推导出了$f(x,y)$的数值解，这个方法就是鼎鼎大名的显示欧拉方法，也可以叫做一阶向前差商近似。

我们来做一个小练习，方程出处为《微分方程的数值解法》例 2.1.1，方程的精确解是$y=\sqrt{1+2x}$

这里给出由 AI 产生的 matlab 代码

 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

clc; clear; close all;

% 参数设置
h = 0.1;
x = 0:h:1;
N = length(x);

% 初始化
y_euler = zeros(1, N);
y_euler(1) = 1;          % 初值 y(0)=1

% 显式欧拉法
for n = 1:N-1
    y_euler(n+1) = y_euler(n) + h * ( ...
        y_euler(n) - 2*x(n)/y_euler(n) );
end

% 解析解
y_exact = sqrt(1 + 2*x);

% 误差
error = y_euler - y_exact;

% ===== 作图 =====
figure;

subplot(2,1,1)
plot(x, y_exact, 'k-', 'LineWidth', 2); hold on;
plot(x, y_euler, 'ro--', 'LineWidth', 1.5);
grid on;
xlabel('x');
ylabel('y');
legend('Exact solution', 'Explicit Euler');
title('Explicit Euler vs Exact Solution');

subplot(2,1,2)
plot(x, error, 'b*-', 'LineWidth', 1.5);
grid on;
xlabel('x');
ylabel('Error');
title('Numerical Error (Euler - Exact)');

可以观察到，显示欧拉求解的效果并不好，而且误差明显随着步长的前进在逐步增加，读者可以试试，最后结果会发散，这明显是不能接受的。另外还有一个问题，初值不知道怎么办？显然，如果初值给错了，在发散的数值解法下，所有解都没有意义。我们需要对算法进行改进。

复盘一下，我们推导的源头相当于从$y^{\prime}(x_i)$入手，自然而然的联想到$y^{\prime}(x_i)\approx\frac{y(x_{i+1})-y(x_i)}{h}$，但是实际上，基于泰勒公式，$y^{\prime}(x_i)$还可以有其他的形式，比如：

我们把一阶向后差商改写一下，变成$y^{\prime}(x_{i+1})\approx\frac{y(x_{i+1})-y(x_{i})}{h}$，整理一下就是所谓的隐式欧拉公式：

称它为隐式，是因为$y(x_{n+1})$在等式右边，需要迭代求解。初学者学到这里可能就很困惑，我们先不用管迭代求解是什么，总之先记住，解包含$y(x_{n+1})$的方法统称为隐式方法即可。

而二阶中心差商，出现了$y(x_{i-1})$，通俗的理解就是，这个方法不仅需要我们提供初值，还需要提供一个历史项，这里不深入讨论，这类方法一般叫做多步法。

目前为止，还没有出现大家耳熟能详的“梯形法则”，那么我们接下来讨论从积分角度对 (1) 求解。

积分近似 —— 由求积分的思想引出梯形法则

梯形积分法则很容易理解，这里借用百度百科的图示，一目了然，使用梯形面积去近似积分：

但是梯形积分法其实远远没有你想的那么简单，在非线性 ODE 解法中，完整形态的梯形法则其实是一种隐式的需要迭代的算法。而在数字软件仿真环境下，往往不会使用隐式迭代的算法，而是采用双线性变换去离散近似积分。这里读者有一个印象就好，暂时不需要深入理解，我们后边会给出具体例子。

既然是积分近似，我们就要消除求导项，所以对 (1) 两端进行积分：

显然，使用梯形法则，可以将右侧积分项近似为：

终于，可以得到新的表达式：

但是很明显，等式右侧有$y_{i+1}$，这个是我们要求的，而我们的解，却又包含了我们本身要求的未知数，这似乎不合理？数学上定义这种情况，或者说这种方法，为隐式方法，通常需要迭代求解，我们来解释一下什么意思。

熟悉 Dommel 算法的读者，可能就会产生疑问，似乎套用梯形法则的时候，从来不存在迭代，隐式的情况（这个其实是笔者初学的问题）？笔者私以为 Dommel 算法其实是基于双线性变换，它虽然也是基于梯形法则，但是是线性 ODE 的特例，所以不存在“隐式”，“迭代”的情况（这个后边给例子）。因此简单的把它归类为“梯形法则”其实是不严谨的，因为梯形法则其实是需要基于迭代的隐式算法，理论上它可以求解任何非线性 ODE，只是 RLGC 组成的电路应用了它而已。

我们观察一下这张图，可以发现$f(b)$，或者$f(x_{i+1},y_{i+1})$，是我们未知的，但是$f(a)$，或者$f(x_{i},y_{i})$，是我们已知的。我们虽然不能求解出精确的$y_{i+1}$，但是不是还有欧拉法么，我们可以求解出一个不是那么精确的“$y_{i+1}$”，记为$y^{(0)}_{i+1}$，这样一来，套用梯形法则的条件就完备了！套用一次梯形法则，我们可以得到精确了一点的$y^{(1)}_{i+1}$，再以此类推，得到$y^{(2)}_{i+1}$、$y^{(3)}_{i+1}$。。。，一直到我们满意为止。整个过程如上图所示，这个就是所谓的迭代。我们直接上 matlab 代码来描述这个过程：

 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

clc; clear; close all;

% ============================
% 参数设置
% ============================
a = 0; b = 1;     % 求解区间
N = 10;          % 步数，保证稳定
h = (b-a)/N;

x = zeros(1, N+1);
y = zeros(1, N+1);

for i = 1:N+1
    x(i) = a + (i-1)*h;
end

y(1) = 1.0;          % 初值
maxerr = 0.0;
epsilon = 1e-4;      % Picard 收敛精度

% ============================
% 主循环：梯形法 + 迭代
% ============================
for i = 1:N

    % ---- 欧拉预测 ----
    ytemp1 = y(i) + h * f(x(i), y(i));

    k = 0;
    % ---- 迭代 ----
    while true
        k = k + 1;
        if k ~= 1
            ytemp1 = ytemp2;
        end
        ytemp2 = y(i) + h/2 * ( f(x(i), y(i)) + f(x(i+1), ytemp1) );
        fprintf('i=%d, k=%d, |ytemp1 - ytemp2| = %.6e\n', i, k, abs(ytemp1 - ytemp2));
        if abs(ytemp1 - ytemp2) < epsilon
            break;
        end
    end

    y(i+1) = ytemp2;

    % ---- 误差 ----
    err = abs(y(i+1) - yexact(x(i+1)));

    if err > maxerr
        maxerr = err;
    end
end

fprintf('The max error is %.6e\n', maxerr);

% ============================
% 解析解 & 误差
% ============================
y_exact = yexact(x);
error = y - y_exact;

% ============================
% 作图（上下两幅）
% ============================
figure;

subplot(2,1,1)
plot(x, y_exact, 'k-', 'LineWidth', 2); hold on;
plot(x, y, 'bo--', 'LineWidth', 1.5);
grid on;
xlabel('x'); ylabel('y');
legend('Exact solution', 'Implicit Trapezoidal (Picard)');
title('Implicit Trapezoidal Method vs Exact Solution');

subplot(2,1,2)
plot(x, error, 'r*-', 'LineWidth', 1.5);
grid on;
xlabel('x'); ylabel('Error');
title('Numerical Error (Trapezoidal - Exact)');

% ============================
% 函数定义
% ============================
function val = f(x, y)
    val = y - 2*x/y;
end

function val = yexact(x)
    val = sqrt(1 + 2*x);
end

可以观察到，精度已经大大提升了，为了节约计算资源，迭代跳出循环的条件可以改为一定次数（比如 50 次）或者更改最大误差容忍度。工程上往往只迭代一步，即用欧拉估算一步以后，在梯形法则一步，称为“改进的欧拉算法”。

到目前为止，我们已经对微分方程的数值算法有了一个初步的认识。深入了解的话，基于泰勒公式的高阶展开，会演变出各种高阶的算法，比如大名鼎鼎的 Runge-Kutta 算法，用高 - 低阶精度算法混合的变步长算法等等，感兴趣的读者可以延申阅读《微分方程的数值解法》。

但是作为电磁暂态仿真、电气仿真工程师的我们，目前入门的知识已经暂时够用了，数学的难题有数学专家去解决，我们接下来要讲解如何工程化，下一篇文章将讲解数字仿真环境下的离散算法，进而引出大名鼎鼎的 Dommel 算法。

想说的话

• 读完本文，你可能对具体实操求解微分方程还是懵懵懂懂，这个没关系，最重要的是理解可以从“差商”和“积分”两种近似角度，去求解带$\frac{\mathrm{d}y}{\mathrm{d}x}$的方程。而目前大部分的数值解法的思想源头都是泰勒展开，目前只要了解了这个，就足够了。

• 笔者认为理解最原始 ODE 求解方法是很重要的，因为电磁暂态仿真本质上是对电感 L，电容 C引入的微分方程的求解。同时也要了解“梯形法则”最原始的形态是什么样的，它是可以求解非线性ODE 的隐式算法，而 Dommel 算法其实是线性的特例。

• 感兴趣的读者，可以试试用差商近似的方法，看看能不能推导出基于积分思想得到的梯形法则算法，即公式 (2)。这个非常有用，它涉及到后边会提到的一个很“神奇”的积分方法的推导思想，数学界称为$\theta$-method，《电磁暂态仿真》称为“可调积分”，业界也有人称为$\alpha$积分方法，它可以解决 Dommel 引入的数值震荡问题，这个将来会单开一篇讲。

• 《微分方程的数值解法》例 2.1.1，也就是本文用的例子，笔者认为不是特别好，因为这个方程是一个不稳定方程，它并不适合初学者，任何复现原著代码的人都会发现问题，感兴趣的读者可以试试，增加仿真时间，你会发现连 ode45 都解不了，这个方程的数值误差是累积的，最后都会发散。如果换用隐式欧拉，甚至还会出现步长越小，越“发散”的现象，这个困惑了笔者很久。但是笔者最终选用了这个例子（一方面是因为当初自学的时候，搜遍全网，有人和我有一样的困惑，但是没有好的解答），是因为这个例子能很好的理解“收敛性”和“稳定性”，并且可以打破我们常规认识：“步长越小越好”，进而引出绝对稳定域的概念，在特殊情况下，如果稳定域是一个空心圆，步长小反而会发散。这里读者可以先不用管，有个印象就好，以后单开一篇讲。

从工程仿真资源利用率的角度看，选用合适的算法，合适的步长，才能达到“精度”和“速度”的最佳平衡点，而这都需要坚实的数学理论做基础。