第一讲微分方程

第一节 微分方程的数值计算

1-1-1 引言

电力系统的动态过程通常用微分方程来描述。在微分方程中，自变量是时间，因变量是系统中各物理变量。对于n 阶常微分方程，可以采用引入新变量的方法把它转化为n 个一阶常微分方程来求解。例如，对于二阶常微分方程

22

2d y K y dt

=- （1-1）

可以表示成：

2

dz K y dt

dy z dt

⎧=-⎪⎪⎨

⎪=⎪⎩ （1-2）

一般来说，依据电工学和动力学建立起来的电力系统数学模型是可以用来数字仿真的。但是，在复杂电力系统中，由于数学模型很复杂，仿真的计算量很大，影响了仿真效率。有些环节对于所研究的问题影响很小，在建立数学模型时可以把它们忽略。这种简化严格来说应该用实际实验，或者用精确的物理模拟实验来检验，以确定简化的合理性。数学模型中的参数是数字仿真的基础，一般可以通过各种测量和实验得到，但这是一件十分烦琐的工作。目前，有些参数只能凭经验给出，其合理性有待实验的严格检验。

电力系统仿真的另一个重要问题是数学模型的求解问题。对于微分方程，除少数可以得到解析解以外，大多数只能采用数值解法。早在18世纪末，就有很多人提出过求解微分方程精确而有效的数值积分方法。

微分方程的积分是一簇曲线。通常，在初值确定之后，数值解才能够确定。

得到微分方程的数值近似解有两种基本的方法。一种方法是把近似解表示成有限个独立函数之和，另一种方法是差分法。

差分方法是寻求在一系列离散点上的近似值，这些离散点称为结点。在多数情况下，这些结点是等距的，即

n t nh =

（1-3）

h 称为步长。差分方法是一种递推算法，它使求解过程能顺着结点的顺序一步一步向前

推进，即可用前一个结点上的值(单步法)或者前面几个结点上的值(多步法)来计算当前步上的近似值。 1-1-2 单步法

所谓单步法就是用前一结点上的值来计算当前结点上的近似值的方法。 1．欧拉法

欧拉法是一种最简单的单步法。

对一阶微分方程

(,)dy

f t y dt

= （1-3）

假定0y 已给定，可计算出00(,)f t y 。把一阶微分方程写成数值积分法的计算形式

1

10(,)t t y y f t y dt =+⎰

（1-4）

如果1t 十分靠近0t ，则有

100010(,)()y y f t y t t =+-

（1-5）

当然，也可以把一阶微分方程写成

10

0010

(,)y y f t y t t -=-

（1-6）

即用0t 时的导数表示()y t 在区间[0t ，1t ]内的平均增量。 结点1n t +和结点n t 之间的关系可表示为

11(,)()n n n n n n y y f t y t t ++=+-

（1-7）

尽管这种方法精度不高，但却提出了一种设想，可以通过递推计算方法得出微分方程的近似解，这种方法称为欧拉法。因为是用折线来近似表示曲线(,)f t y ，故又可称为折线法。

2．隐式梯形法

从欧拉法可以看出，为提高精度，可以用00(,)f t y 和11(,)f t y 的平均值作为()y t 在区间[0t ，1t ]内的平均增量值。

10110010(,)(,)

2

y y f t y f t y t t -+=-

（1-8）

11001010(,)(,)

()2

f t y f t y y y t t +=+

-

（1-9）

在等式的两边均含有未知量1y ，故称为隐式梯形法。隐式梯形公式不能用递推的方式直接做数值计算。如果用欧拉法求得的值作为预测值，则有：

10001011001010(,)()(,)(,)

()2

y y f t y t t f t y f t y y y t t =+-⎧⎪

⎨+=+-⎪⎩ （1-10）

这种方法称为改进欧拉法。

利用欧拉法和改进欧拉法求微分方程的数值解，如果选取步长相同，改进欧拉法的计算量大，但改进欧拉法的精度高。反过来说，如果两种方法要求精度相同，则改进欧拉法可以选取较大步长，总的计算量可以节省，舍入误差也可较小。

3．龙格一库塔法

隐式梯形法是把00(,)f t y 和11(,)f t y 的平均值作为()y t 在区间[0t ，1t ]内的平均增量值。同理可以设想在[0t ，1t ]区间内，取不同的(,)f t y 后再求加权平均值，用它作为该区间的平均增量，即

11

N

n n i i i y y R K +==+∑

（1-11）

其中

1001

1

01

1

(,)

(,)

2,3,,i i i ij n ij j j j K hf t y K hf t h B y B K i N

--====++=∑∑

1

1N

i

i R

==∑

这种计算式称为自启动的N 阶龙格一库塔法。当N =l 时为欧拉法，当N =2时为改进欧拉法。常采用的各阶龙格一库塔法计算公式的系数如表l 所示。

表1 各阶龙格一库塔法的系数表

N

ij B

i R

通常采用的名称

1 0

11

R =

欧拉法

2 1212

211

B B =

=

120,1R R == 1122R R == 改进欧拉法 3

1212

31321,2

B B B ==-=

1126

233

R R R ===

三阶龙格一库塔法

4 1212

131322

4142430,0,0,1

B B B B B B ===

===

11461233

R R R R ====

四阶龙格一库塔法

四阶龙格一库塔法

112341

(22)6

n n y y K K K K +=++++

（1-12）

其中

1(,)n n K hf t y =

2111

(,)22n n K hf t h y K =++

3211

(,)22n n K hf t h y K =++

43(,)n n K hf t h y K =++