http://www.ma-xy.com
目录
第一章 常微分方程 1
1.1 常微分方程建模 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 数学摆的微分方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.2 传染病模型 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.3 人口增长的微分方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 常微分方程基本理论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3
常微分方程数值方法
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.1 常微分方程初值问题 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3.2 常微分方程边值问题 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3.3 MATLAB 解常微分方程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4 常微分方程稳定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.1 李雅普诺夫意义下的稳定 . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.4.2 李雅普第二方法 (直接法) . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
1.5 混沌理论 Chaos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.5.1 混沌系统的引入 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
1.5.2 混沌的基本理论 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.5.3 混沌的基本性质 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
1.5.4 混沌的数值方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.5.5 混沌系统的应用 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1
http://www.ma-xy.com
第一章 常微分方程
1.1 常微分方程建模
1.1.1 数学摆的微分方程
数学摆是质量为 m 的物体 M 系于一长为 l 的线上，在重力的作用下沿圆周运动，如图 (1.1)
所示，求物体的运动方程。
图 1.1: 单摆示意图
解：像一般分析的那样，先假设，再确定变量，然后是建立坐标系、受力分析，最终我们会
有它的运动方程。
假设 1: 质量为 m 的物体 M 仅受重力和拉力作用，忽略其它力；假设 2：物体 M 在初始点
处的速度 v 是 0；假设 3：物体 M 视为质点；假设 4：物体初始位置在固定点水平右方 (保持细
线拉直)；就此题来看，我们应建立二维坐标系，并且用极坐标系可能更方便一些，因为我们有了
极径 l。
建系：以线端点 o 为坐标原点，水平向右方向为 x 轴正向，竖直向下 (重力方向) 为 y 轴方
向，建立平面直角坐标系 xoy.
设置变量：设线长为 l ，线与 y 轴夹角为 φ，物体 M 的质量为 m，初始时刻为 t
0
(不妨设
物体 M 的坐标为 (x, y), 其实不必要)。
受力分析：物体 M 的受力分析示意图如图 (1.2) 所示，我们要在受力分析的基础上，确定
各时刻 t 的 φ 的大小，即确定曲线 φ(t)。
1
http://www.ma-xy.com
1.1 常微分方程建模 第一章 常微分方程
图 1.2: 单摆受力分析
其在 M 点处受重力 G 和拉力 T 的作用，速度为 v，有
G = mg
G cos φ = T
所以是 G sin φ 提供加速度，故
˙v = G · sin φ/m = g · sin φ
而由运动方程 v = ω · l，我们有
v = l · ˙φ (φ(t
0
) = φ
0
)
所以，可以确定最终的运动方程为
l · ¨φ = g · sin φ
¨φ = g/l · sin φ
或者写为
d
2
φ
dt
2
= g/l · sin φ (1.1)
给定初始条件：φ(t
0
) = φ
0
，其中：t
0
是初始时间，φ 为初始 (刚开始下落时) 的夹角大小。
此时，我们的目标就变为求一个 φ(t) 函数，φ(t) 应满足方程 (1.1) 以及 φ(t
0
) = φ
0
的初始条件。
1.1.2 传染病模型
第一次参加建模是 2014 年小美赛，当时要处理的问题正是传染病“埃博拉”。经过查找文献，
找到了 SIR 模型，懵懵懂懂的将模型改进为 SIRAB(当时实在不知道该用个啥样的名字)，但建
立完 SIRAB 之后，问题也随之而来，这是个什么模型，怎么解？
对于传染病而言，我们更关心的是随着时间 t 的增长，感染的人数是如何变化的，如何对传
染病进行控制，所以至少我们会有下面的图 (1.3)
http://www.ma-xy.com 2 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.1 常微分方程建模
图 1.3: 感染人数增长示意图
我们不妨将其离散化，问题是：如果我们知道了 t 时刻的感染人数，那 t + 1 时刻的人数会
是多少。我们考虑一个地区 A，A 在 t 时刻共有 N(t) 个人，如果地区 A 不是封闭的话，t + 1
时刻 N (t) 可能增加也可能减少。为了简单，我们假设地区 A 的人口总数恒为 N 。如果我们知道
了 A 地区 t 时刻的染病人数，如何求 t + 1 时刻的染病人数？设染病人数为 I(t)，那么未感染人
数就变为了 N − I(t)，假设未感染的人群都是易感染者，记为 S(t)(当然，实际情况可能是未感
染人群中部分人已经有免疫能力了，那么他就不易变为感染者)。我们建立下面的传染病 SI 模型
I(t + ∆t) = I(t) + αS(t)∆t
S(t + ∆t) = S(t) − αS(t)∆t
I(t) + S(t) = I(t + 1) + S(t + 1) = N
其中：α 为单位时间内未感染者变为感染者的可能，记为发病率。将上式转化为
˙
I(t) = αS(t) ≈ ∆I(t)
˙
S(t) = −αS(t) ≈ ∆S(t)
I(t) + S(t) = N
为方便理解，
˙
I(t) ≈ ∆I(t) 记为 t 时刻传染者新增的人数，αS(t) 为 t 时刻未感染者 S(t) 转化为
感染者的人数。
对于上面的 SI 模型，我们有许多可以改进方法，比如，我们可以做如下改进：
(1) 前面，我们假设地区 A 在个时刻 t 内仅有两种人：感染者 I 和未感染者 S。实际上，还
可能有感染之后恢复的人，并且这类人不会再染病，我们记这类人为 R，设染病者的恢复率为 β，
则 A 地区不同人群的人口转移如图 (1.4) 所示
图 1.4: SIR 感染者转移示意图
由此可建立微分方程模型
˙
I(t) = αS(t) − βI(t)
˙
S(t) = −αS(t)
˙
R(t) = βI(t)
I(t) + S(t) + R(t) = N
http://www.ma-xy.com 3 http://www.ma-xy.com
http://www.ma-xy.com
1.1 常微分方程建模 第一章 常微分方程
(2) 还有一种可能是：在 [t
0
− t
1
] 时间段内不存在 R 类人，在 t
1
时刻 (疫苗到达时) 出现了 R 类
人。
(3) 我们还可以假设这个传染病是有潜伏期的，未感染者 S 先变为潜伏者 E(E 不表现为发
病)，然后潜伏者才变为感染者 S，其转移情况如图 (1.5) 所示
图 1.5: SIRE 感染者转移示意图
(4) 我们在来假设一部分感染者 S 死亡了，记死亡者为 D，可以建立 SEIRD 传染病模型，
其转移情况如图 (1.6) 所示
图 1.6: SEIRD 感染者转移示意图
当然，我们可以将人群的分类分的越来越细。在建完传染病模型之后，我们要考虑模型中的
α, β, γ 等参数如何求解，暂时将模型中的参数统一记为 θ ，下面，我们分两种情况来讨论 θ 的求
解：
1. 如果我们已经有各种类型人的数据了，比如：已经获得了 t = [1, T ] 时刻的 S, I, R 三类人
的人数，那么，我们可以将数据带入建立的 SIR 模型当中来求解参数 θ，这是一个曲线拟
合的工作。
2. 如果我们没有获得数据，可以分析上面系统 (SIR) 的稳定性，上面建立的模型明显是一个
线性系统，我们可以在系统分析方面做些文章。
1.1.3 人口增长的微分方程
有前面的两个模型做引导，人口增长模型本不应该再讨论，但是，在后面介绍的其它方程当
中 (偏微分方程、积分方程和随机微分方程等)，我们都需要拿人口增长模型来做引导，所以，这
里也简单介绍一下。
离散时间来看，已知 t 时刻的人口数，需要求 t + 1 时刻的人口数。设 t 时刻人口数为 x(t)，
我们可以想到 [t, t + 1] 时刻内，新增人口 ∆x 应该与 x(t) 有关，x(t) 越大，∆x 也越大，即
∆x ∝ x(t)。不妨设其关系为线性关系，线性关系的系数设为 α，即 ∆x = αx(t)，于是有
x(t + 1) = x(t) + αx(t)
上面是在离散时间上讨论的，下面，我们在连续时间 [t, t + ∆t] 内进行考虑，当 ∆t 很小时，有
x(t + 1) ≈ x(t) + αx(t)∆t
http://www.ma-xy.com 4 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.2 常微分方程基本理论
故
˙x = αx(t)
x(t
0
) = x
0
或者写为
dx
dt
= αx(t)
x(t
0
) = x
0
其中：t
0
为初始时刻，x
0
为人口基数，或者初始人口。
1.2 常微分方程基本理论
用函数与方程的思想将上面的模型泛化，有
F (x, y, y
′
, y
′′
, . . . ) = 0
其中：F 为函数。
对上面的这类常微分方程 (ODE) 问题，我们有解析方法和数值方法两种方法进行求解。下
面，我们将介绍一些解析法中的基本概念即解的存在唯一性。我们将微分方程分类，以便于对不
同类型的方程进行求解，仅就方程而言：
1、常微的阶数。常微中 y 的最高阶导数的阶为常微的阶。例如：y
′′′
+ y
′′
+ y
′
+ y = 0 的阶
数为 3。
2、常微的显隐性。F (x, y, y
′
, . . . , y
(n)
) = 0 为 n 阶隐式方程；y
(n)
= f(x, y, . . . , y
(n−1)
) 为 n
阶显式方程。
3、常微的线性与非线性。y
′′
+ ay
′
+ by = f(x) 为线性常微；y
′′
+ ae
y
′
+ b cos y = f(x) 为非
线性常微。
4、常微的常系数与变系数。y
′′
+ ay
′
+ by = f(x) 为常系数；y
′′
+ a(x)y
′
+ b(x)y = f(x) 为
变系数。
5、常微的齐次性与非齐次性。y
′′
+ ay
′
+ by = 0 是齐次 (自治)，方程变化仅依赖 y 自身；
y
′′
+ ay
′
+ by = f(x) 是非齐次，y
′
依赖于 x。
6、初边值条件。初值条件用于指定某常微系统的初始条件，如 y(0) = y
0
；边值条件用于指
定某常微系统在时刻 a, b 的系统信息 (可以是 y(a) 的值，也可以是其导数值 y
′
(a) 等等)，例如：
y(a) = α; y(b) = β。
7、解的显隐性：设方程解为 φ(x, y)，φ 为显函数则为显式解。
上面，我们介绍了常微分方程以及常微分方程的一些基本定义，后面，我们应该考虑的问题
是：1、常微分方程解是否存在；2、解存在的话是否唯一；3、是否有解析解；4、数值求解方法；
5、数值解的存在唯一性；6、解的稳定性等等。
http://www.ma-xy.com 5 http://www.ma-xy.com
http://www.ma-xy.com
1.2 常微分方程基本理论 第一章 常微分方程
8、解的存在唯一性。考虑下面的常微问题解的存在唯一性
˙
y
=
f
(
x, y
)
y(x
0
) = y
0
如果右端 f (x, y) 在闭矩形域 R
R = {(x, y)|x
0
− a ⩽ x ⩽ x
0
+ a, y
0
− b ⩽ y ⩽ y
0
+ b}
上满足如下条件：1、f 在 R 上连续；2、f 在 R 上关于变量 y 满足利普西茨 (Lipschitz) 条件，
即存在利普西茨常数 L，使得对 R 上任何两个点 (x, y) 和 (x, y)，有下面的不等式成立
|f(x, y) − f(x, ¯y)| ⩽ L|y − ¯y|
则常微在区间 x
0
−h
0
≤ x ≤ x
0
+h
0
上存在唯一解 y = φ(x)，φ(x
0
) = y
0
。其中：h
0
= min
a,
b
M
，
M = max
x,y ∈R
|f(x, y)|。
定理 (解对初值的连续性定理) 若 f 在区域 G 内连续，而且 f 关于 y 满足局部利普西茨条
件，则 y = φ(x) 在它存在的范围内连续。
定理 (解对初值的可微性定理) 若 f 及
∂f
∂y
在区域 G 内连续，则解 y = φ(x) 是连续可微的。
9、常微分方程组
有了方程，多个方程组合起来就可以形成方程组。既然如此，我们也可以将方程组的概念引
入到常微分方程当中，形成微分方程组。前面建立的 SIR 模型就是微分方程组的形式，我们再给
一个微分方程组的示例：
dx
dt
= ˙x = f
1
(t, x, y, z)
dy
dt
= ˙y = f
2
(t, x, y, z)
dz
dt
= ˙z = f
3
(t, x, y, z)
上面的例子是一个 3 元方程组。我们将其时间量 t 离散化来看，如果有上面方程组的测量数
据，其形式应该如下表 (1.1) 所示
表 1.1: 数据形式示例
t t
1
t
2
. . . t
n
x · · · ·
y · · · ·
z
· · · ·
http://www.ma-xy.com 6 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.3 常微分方程数值方法
将上面的例子泛化，写出常微分方程组的一般形式，有
˙x
1
= f
1
(t, x
1
, x
2
, . . . , x
n
)
˙x
2
= f
2
(t, x
1
, x
2
, . . . , x
n
)
. . .
˙x
n
= f
n
(t, x
1
, x
2
, . . . , x
n
)
(1.2)
其中：x
i
(i = 1, 2, . . . , n) 为待求的关于时间 t 的函数；f
i
(i = 1, 2, . . . , n) 为已知的函数形式。
将上面的微分方程组 (1.2) 写为矩阵形式，有
˙
x = F(t, x)
其中：F 为标量函数。
1.3 常微分方程数值方法
1.3.1 常微分方程初值问题
我们先从简单的一阶常微分方程来看。常微分方程的初值问题描述为：求一个函数 y，y 在
x
0
= a 处的值为 y
0
，而且 y 满足方程 y
′
(x) = f (x, y(x))，即
y
′
= f(x, y)
y(a) = y
0
a < x < b
其解的示意图如图 (1.7) 所示
图 1.7: 常微初值问题解的示意图
注：y 是 x 的函数，y
′′
亦为 x 的函数。
对于此问题，我们应该讨论问题解的存在性与唯一性 (这个之前介绍过)。在解存在唯一的基
础上，我们来考虑其数值解法。下面，我们介绍初值问题的欧拉法。
我们先将 x 在 [a, b] 区间内离散化为 n 个点，这里为了方便，我们实现 n 等分。令 h =
b−a
n
，
有
x
i
= a + i
b − a
n
i = 0, 1, . . . , n
假设已知 x
i
, y
i
，现在的问题是如何求 y
i+1
。
y
′
(x
i
) = f (x
i
, y(x
i
))
http://www.ma-xy.com 7 http://www.ma-xy.com
http://www.ma-xy.com
1.3 常微分方程数值方法 第一章 常微分方程
而
y
′
(x
i
) =
y(x
i+1
) − y(x
i
)
h
−
h
2
y
′′
(ξ
i
)
故
y
i+1
≜ y(x
i+1
) = y(x
i
) + hf (x
i
, y
i
) +
h
2
2
y
′′
(ξ
i
)
由此，可以进一步思考推广：既然 y
′
(x) = f (x, y)，那么，我们对等式两边进行相同的操作
后，等式仍然成立。接下来，我们应该考虑的是：对于同一个问题，我们会有各种不同的解法，我
们应该如何评价它们呢？对于此问题，这里不做讨论。
常用的常微初值问题的数值解法是 (4 阶)Runge-Kutta 方法，R-K 方法是单步方法，即求
y
i+1
时只需要有 y
i
即可，其计算精度高，但计算量大，我们在后面的“嫦娥 3 号”最优控制的
直接参数解法中，构建非线性方程时引入了 R-K，所以这里不做介绍。
1.3.2 常微分方程边值问题
为简单，我们来看一个二阶常微分方程边值问题：求一个函数 y，y 在 x = a 处有已知信息，
在 x = b 处亦有已知信息，信息可能是 y 在 a 处的一阶导、二阶导或者三阶导，并且要求 y 满
足方程 y
′′
= f(x, y, y
′
)，a < x < b。例如：
y
′′
= f(x, y, y
′
)
y(a) = α
y(b) = β
a < x < b
下面，来介绍常微边值问题的数值解法 - 差分法。该方法用数值微分公式代替微分方程和边
值条件中的导数，略去误差项，将常微离散为一个差分方程组，之后将方程组的解作为方程的近
似解。
在点 (x
i
, y
i
) 处有
y
′′
(x
i
) = f (x
i
, y
i
, y
i
(x
i
))
而数值微分公式为
y
′′
(x
i
) =
y(x
i−1
) − 2y(x
i
) + y(x
i+1
)
h
2
−
h
2
12
y
(4)
(ξ
i
)
将二者结合，略去高阶项 y
(4)
(ξ
i
)，有
y(x
i−1
) − 2y(x
i
) + y(x
i+1
)
h
= f(x
i
, y
i
, y
′
i
) i = 1, 2, . . . , n − 1
上式即为差分方程组。
当 y 的边值信息已知时 (y(a) = α, y(b) = β)，则差分方程组为
y(a) − 2y(x
1
) + y(2)
h
= f(x
1
, y
1
, y
′
1
) i = 1
y(x
i−1
) − 2y(x
i
) + y(x
i+1
)
h
= f(x
i
, y
i
, y
′
i
) i = 2, 3, . . . , n − 1
y(x
n−2
) − 2y(x
n−1
) + y(x
b
)
h
= f(x
i
, y
i
, y
′
i
) i = n
http://www.ma-xy.com 8 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.3 常微分方程数值方法
解上述方程组，得到 {y
i
}
n
i=1
即为 y 的离散值。
1.3.3 MATLAB 解常微分方程
常用函数说明
1、MATLAB 中提供 ode23、ode45、ode113、ode15s、ode23s、ode23t、ode23tb 函数用于
求解常微分方程 (组) 初值问题。关于函数的说明，可以参考表 (1.2)
表 1.2: MATLAB 常微初值问题函数表
函数 问题类型 精确度 说明
ode45 非刚性 中等 采用算法为 4-5 阶 Runge-Kutta 法，是大多数情况下首选的
函数
ode23 非刚性 低
基于 Bogacki-Shampine2-3 阶 Runge-Kutta 公式，在精确度
要求不高的场合，以及对于轻度刚性方程，ode23 的效率可能
好于 ode45
ode113 非刚性 低到高
基于变阶次 Adams-Bashforth-MoutlonPECE 算法。在对误差
要求严格的场合或者输入参数 odefun 代表的函数本身计算量
很大情况下比 ode45 效率高。ode113 可以看成一个多步解算
器，因为它会利用前几次时间点上的解计算当前时间节点的
解。因此它不适应于非连续系统
ode15s 刚性 低到中
基于数值差分公式 (反向差分公式，BDFs 也叫 Gear 方法)，
因此效率不是很高。同 ode113 一样，ode15s 也是一个多步解
算器，当 ode45 求解失败，或者非常慢，并且怀疑问题是刚性
的，或者求解微分代数问题时可以考虑用 ode15s
ode23s 刚性 低
基于修正的二阶 Rosenbrock 公式。由于是单步解算器，当精
度要求不高时，它的效率可能会高于 ode15s。ode23s 可以解
决一些 ode15s 求解起来效率不太高的刚性问题
ode23t 适度刚性 低 ode23t 可以用来求解微分代数方程
ode23tb 刚性 低 当方程是刚性的，并且求解要求精度不高时可以使用
说明：常微的刚性和非刚性，直观上讲，在短时间内常微解 y(x) 有巨大波动变化的，称为
刚性问题 (sti)。对刚性问题，如果使用 ode45 的话，一般会花费很大代价。
上述表 (1.2) 中函数命令的一般格式为
[t,y] = ode∼(odefun,tspan,y0,options)
其中：
odefun
为方程函数句柄；
tspan
为方程的时间长度；
y0
为初始条件；
options
为参数设置。
2、MARLAB 中提供 dde23 和 ddesd 函数用于处理延迟常微分方程 (DDE)。前者用来求解
状态变量延迟数为常数的微分方程 (组)，后者用来求解延迟非常数的微分方程 (组)，而 ddensd
http://www.ma-xy.com 9 http://www.ma-xy.com
http://www.ma-xy.com
1.3 常微分方程数值方法 第一章 常微分方程
用于求解中性类型延迟微分方程。延迟微分方程在常微分方程的基础上添加时延，一般形式为
y
′
(t) = f (t, y(t), y(t − τ
1
), y(t − τ
2
), . . . , y(t − τ
n
))
其中：τ
1
, τ
2
, . . . , τ
n
> 0 是时延，y(t − τ
1
), y(t − τ
2
), . . . , y(t − τ
n
) 是时延项。
延迟微分方程中，y 在 t 时刻的状态增量 y
′
(t) 不仅与此时的状态 y(t) 有关，而且还与之前
的状态 y(t − τ
1
), y(t − τ
2
), . . . , y(t − τ
n
) 有关。如果你不是很熟悉 DDE，可以继续阅读后面章节，
在后面章节中你会慢慢体会到延迟的思想。
dde23 函数命令格式为
sol = dde23(ddefun,lags,history,tspan,options)
其中：ddefun 为时滞方程函数句柄；lags 为时间延迟向量，其中，τ
k
保存在 lags(k) 中；history
是描述 t ≤ t
0
时状态变量值的函数，可以是句柄也可以是常数；tspan 是时间区间；options 是参
数设置结构体；sol 为输出结构体，包括：sol.x，sol.y，sol.yp，solver 等。
ddesd 函数命令格式为
sol = ddesd(ddefun,delays,history,tspan,options)
其中：上述这个函数可以处理的延迟微分方程的一般形式为：
y
′
(t) = f (t, y(t), y(d(1)), y(d(2)), . . . , y(d(k)))
ddefun 为时滞方程函数句柄，例如：dydt=ddefun(t,y,z)，t 对应当前时刻 t，y 是一个 y(t) 的近
似列向量，z 的第 j 个列向量对应 y(d(j)) 的估计；delays 为函数句柄，返回 d(j) 列向量，d(j)
可以依赖于 t 和 y(t)。
ddensd 函数命令格式为
sol = ddensd(ddefun,dely,delyp,history,tspan,options)
其中：上述这个函数可以处理的延迟微分方程的一般形式为:
y
′
(t) = f (t, y(t), y(dy
1
), ..., y(dy
p
), y
′
(dyp
1
), ..., y
′
(dyp
q
))
t is the independent variable representing time；dy
i
is any of p solution delays；dyp
j
is any of q
derivative delays.
3、MATLAB 中提供了 bvp4c，bvp5c 函数用于求解常微分方程边值问题。bvp4c 和 bvp5c
函数具有相同的命令格式
sol = bvp4/5c(odefun,bcfun,solinit,options)
solinit = bvpinit(x, yinit, params)
其中：bcfun 为边值函数句柄；solinit 为结构体：solinit.x 为取值范围，solinit.y 为初始猜测值，
solinit.params 为参数的初始猜测；options 为参数结构体；x 为时刻点，即边值时刻 x=[a,b]，如
果是多点边值，x=[a,c,b]；yint 为解的初始猜测值 y(a)。
http://www.ma-xy.com 10 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.3 常微分方程数值方法
函数示例
这里，我们先给出 Euler 方法的程序。用 Euler 方法求解如下常微分问题
x
′
= −3x + 6t + 5
x(0) = 3
上述问题的解析解为 x = 2e
−3t
+ 2t + 1。求解程序如下
1 L=1; h=0.05; t=0:h : L ;
2 x_Euler = z e r os ( l ength ( t ) ,1) ; x_Euler ( 1 ) =3;
3 x_pc = x_Euler ; x_ode45=x_Euler ;
4 f o r n=1: len gth ( t )−1
5 %欧 拉 法
6 x_Euler (n+1) = x_Euler (n)+h*(−3*x_Euler (n)+6*t (n)+5) ;
7 %预 测−校 正 法
8 k1 = h*(−3*x_pc(n) + 6* t (n) +5) ;
9 k2 = h*(−3*(x_pc( n)+k1 ) + 6*( t (n)+h)+5) ;
10 x_pc(n+1) = x_pc( n)+(k1+k2 ) / 2 ;
11 end
12 %ode45
13 [ t , x_ode45]=ode45 (@( t , x) [−3*x+6*t +5] , t , x_ode45 (1) ) ;
14 %解 析 解
15 x_exact = 2*exp(−3*t )+2*t+1;
16 %画 图
17 f i g u r e
18 plot ( t , x_Euler , ’ xk ’ , t , x_pc , ’ ok ’ , t , x_ode45 , ’+k ’ , t , x_exact , ’ k ’ , ’ MarkerSize ’ ,8 , ’ LineWidth ’ ,1)
19 axis ( [ 0 1 2.3 3 . 1 5 ] ) , s e t ( gca , ’ Fon tsize ’ , 8 ) , xla b el t , yl a be l x
20 legend ( ’ 欧 拉 法 ’ , ’ 预 测−校 正 法 ’ , ’ ode45 ’ , ’ 解 析 解 ’ , ’ l oc a ti o n ’ , ’North ’ )
21
求解结果如图 (1.8) 所示
图 1.8: Euler 方法求解常微分方程的结果
下面，我们给出一些 MATLAB 求解常微分的示例：(1) 用 ode45 求解含参常微分方程初值
http://www.ma-xy.com 11 http://www.ma-xy.com
http://www.ma-xy.com
1.3 常微分方程数值方法 第一章 常微分方程
问题
y
′′
1
− µ
1 − y
2
1
y
′
1
+ y
1
= 0
y
1
(0) = 1, y
1
(0)
′
= 0
0 < t < 30
其中：µ > 0 是一个标量参数。设 y
′
1
= y
2
，重写上面的方程，有
y
′
1
= y
2
y
′
2
= µ(1 − y
2
1
)y
2
− y
1
y
1
(0) = 1, y
1
(0)
′
= 0
0 < t < 30
求解程序为
1 tspan = [ 0 , 3 0 ] ;
2 y0 = [ 1 , 0 ] ;
3 odefun = @(mu)@( t , y) [ y(2) ; mu*(1−y (1) ^2) *y (2)−y (1) ] ;
4 mu = 2 ;
5 [ t , y ] = ode45 ( odefun (mu) , tspan , y0 ) ;
6 p l o t ( t , y ( : , 1 ) , ’−o ’ , t , y ( : , 2 ) , ’−o ’ )
7 t i t l e ( ’ S o lution of van der Pol Equation (\mu = 1) with ODE45 ’ ) ;
8 legend ( ’y_1 ’ , ’y_2 ’ )
9
求解结果如图 (1.9) 所示
图 1.9: ode45 求解带参常微分方程
(2) 用 ode15s 求解刚性常微分初值问题。承继上面的问题，当 µ = 1000 时，在短时间内，
http://www.ma-xy.com 12 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.3 常微分方程数值方法
函数值变化非常大，问题变为刚性问题，我们用 ode15s 求解如下刚性问题
y
′
1
= y
2
y
′
2
= 1000(1 − y
2
1
)y
2
− y
1
y
1
(0) = 1, y
2
(0) = 0
0 < t < 3000
求解程序如下
1 mu = 1000;
2 tspan = [ 0 , 3 0 0 0 ] ;
3 y0 = [ 1 , 0 ] ;
4 [ t , y ] = ode15s ( odefun (mu) , tspan , y0) ;
5 p l o t ( t , y ( : , 1 ) , ’−o ’ )
6
求解结果如图 (1.10) 所示
图 1.10: ode15s 刚性常微分初值问题
(3) 用 ode15i 求解隐式常微分方程。上面展示的是显式常微分问题：y
′
= ( t, y) ，并且 ode45
和 ode15s 的微分方程函数句柄 odefun 要求方程为显式。下面，我们将展示隐式微分方程的求
解，隐式常微分问题可以表示为
f(t, y, y
′
) = 0
y(0) = y0
对于隐式微分方程，我们一般用 ode15i 进行求解。作为示例，我们用 ode15i 求解 Weissinger 隐
式问题：
ty
2
(y
′
)
3
− y
3
(y
′
)
2
+ t(t
2
+ 1)y
′
− t
2
y = 0
求解程序如下
http://www.ma-xy.com 13 http://www.ma-xy.com
http://www.ma-xy.com
1.3 常微分方程数值方法 第一章 常微分方程
1 w eis s ing e r = @( t , y , yp) [ t *y^2 * yp^3 − y^3 * yp^2 + t *( t ^2 + 1)*yp − t ^2 * y ] ;
2 t0 = 1 ;
3 y0 = s q r t (3/2) ;
4 yp0 = 0 ;
5 [ y0 , yp0 ] = de c i c ( we i s s i n g e r , t0 , y0 , 1 , yp0 , 0) ;%为 ode15i 提 供 恰 当 的 初 始 值
6 tspan = [ 1 , 1 0 ] ;
7 [ t , y ] = ode15i ( weissinger , tspan , y0 , yp0 ) ;
8 ytrue = s q rt ( t .^2 + 0 .5 ) ;
9 p l o t ( t , y , t , ytrue , ’ o ’ )
10
求解结果如图 (1.11) 所示
图 1.11: ode15i 求解隐式常微分方程
(4) 上面介绍了一些求解 ODE 初值问题的方法，下面我们来求解常微分方程边值问题。求
解含参常微分方程边值问题，作为示例，我们求解 Mathieu 方程
y
′′
+ (λ − 2q cos 2x)y = 0
y
′
(0) = 0
y
′
(π) = 0
y(0) = 1
其中：λ 为参数。求解第四特征值 (q = 5) 的 Mathieu 方程，求解程序如下
1 lambda = 15;
2 q = 5 ;
3 mat4init = @(x) [ cos (4*x) ; −4*si n (4*x) ] ;%
4 s o l i n i t = b v p i n i t ( l i n s p a c e ( 0 , pi , 1 0 ) , mat4init , lambda) ;%求 解 迭 代 初 始 值 。 挑 选 一 个 满 足 y ’( 0 )
=0,y ’ ( p i ) =0,y(0 )=1的 函 数 （ cos4x ） ,−4 s i n (4 x ) 为 其 导 数
5 mat4ode = @(q)@(x , y , lambda) [ y(2) ; −(lambda − 2*q*cos (2*x) ) *y(1) ] ;%微 分 方 程 匿 名 函 数 句 柄
6 mat4bc = @(ya , yb , lambda) [ ya (2) ; yb(2) ; ya ( 1 )−1 ] ;%边 值 条 件
7 s o l = bvp4c( mat4ode ( q ) , mat4bc , s o l i n i t ) ;
8 f p r i n t f ( ’The f o urth eigenvalue i s approximately %7.3 f . \ n ’ , . . .
9 s o l . parameters )
10 xint = l i n s p a ce ( 0 , pi ) ;
http://www.ma-xy.com 14 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.3 常微分方程数值方法
11 Sxint = deval ( sol , xint ) ;
12 f ig u r e , plot ( xint , Sxint ( 1 , : ) )
13 a xis ( [ 0 pi −1 1 . 1 ] )
14 t i t l e ( ’ Eigen funct ion of Mathieu ’ ’ s equation . ’ )
15
求解结果如图 (1.12) 所示
图 1.12: Mathieu 方程
(5) 求解延迟微分方程问题 DDE。固定时滞的时滞微分方程描述为
y
′
(t) = f (t, y(t), y(t − τ
1
), ..., y(t − τ
k
))
其中：τ
i
为时滞长度，y(t − τ
i
) 为时滞项。因为时滞长度为常数，所以被称为固定时滞微分方程。
如果 τ
i
不固定，而是与时间 t 或者 y 有关 (表示为 d(t))，则称为非固定延迟微分方程。求解下
面的固定延迟微分方程：
y
′
1
(t) = y
1
(t − 1)
y
′
2
(t) = y
1
(t − 1) + y
2
(t − 0.2)
y
′
3
(t) = y
2
(t)
0 < t < 5
小于初始时间的历史信息是：当 t ⩽ 0 时，y
1
(t) = 1, y
2
(t) = 1, y
3
(t) = 1。将下面的程序保存为
函数，以便后面的调用
1 function dydt = ddex1de ( t , y ,Z)
2 % 延 迟 方 程 的 函 数
3 ylag1 = Z ( : , 1 ) ;%对 所 有 延 迟 为 delags (1) 的 状 态 变 量 的 求 解
4 ylag2 = Z ( : , 2 ) ;
5 %y1( t −1)： 状 态 变 量 为 y1 ， 延 迟 了 d e l a g s ( 1 ) ， 故 用 ylag1 (1) 表 示 。
6 dydt = [ ylag1 ( 1 )
7 ylag1 ( 1 ) + ylag2 (2)
8 y ( 2 ) ] ;
9 end
10
http://www.ma-xy.com 15 http://www.ma-xy.com
http://www.ma-xy.com
1.3 常微分方程数值方法 第一章 常微分方程
求解程序为
1 tspan = [ 0 , 5 ] ;
2 delags = [1 , 0 . 2 ] ;
3 h is t o ry = ones (3 , 1 ) ;
4 s o l = dde23 (@ddex1de , delags , history , tspan ) ;
5 f i g u r e ;
6 p l o t ( s o l . x , s o l . y)
7 t i t l e ( ’An example o f Wille ’ ’ and Baker . ’ ) ;
8 x la b el ( ’ time t ’ ) ;
9 y la b el ( ’ s o l u ti o n y ’ ) ;
10
求解结果如图 (1.13) 所示
图 1.13: 固定时滞的时滞微分方程
(6) 求解下面非固定延迟微分方程，
y
1
(t)
′
= y
2
(t)
y
′
2
(t) = −y
2
(exp(1 − y
2
(t))) ∗ y
2
(t)
2
∗ exp(1 − y
2
(t))
0.1 < t < 5
其中：延迟函数 d(t) = exp(1 − y
2
(t)) 。该延迟微分方程有解析解 y
1
(t) = log(t), y
2
(t) = 1/t，可
以作为时间小于初始时间 t < 0.1 的历史信息。将下面的代码保存为函数，以便后面调用。
1 function v = ddex3hist ( t )
2 % 历 史 函 数
3 v = [ l o g ( t ) ; 1./ t ] ;
4 end
5 function d = ddex3delay ( t , y )
6 % 延 迟 函 数
7 d = exp (1 − y (2) ) ;
8 end
9 function dydt = ddex3de ( t , y ,Z)
10 % 延 迟 方 程
11 %由 于 只 有 一 个 延 迟 项 ， 因 此Z只 有 一 列 。
http://www.ma-xy.com 16 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.3 常微分方程数值方法
12 %y2( exp(1−y2 ( t ) ) ) 延 迟 了 exp(1−y2 ( t ) ) , y2 又 是 第 二 个 状 态 变 量 ， 因 此 用Z（ 2 ） 表 示 。
13 dydt = [ y (2) ; −Z(2) *y ( 2 ) ^2*exp (1 − y (2) ) ] ;
14 end
15
求解程序如下
1 t0 = 0 . 1 ;
2 t f i n a l = 5 ;
3 tspan = [ t0 , t f i n a l ] ;
4 s o l = ddesd (@ddex3de , @ddex3delay , @ddex3hist , tspan ) ;
5 % Exact s ol u ti o n
6 texact = l i ns p ac e ( t0 , t f i n a l ) ;
7 yexact = ddex3hist ( t e x act ) ;
8 f i g u r e
9 p l o t ( texact , yexact , s o l . x , s o l . y , ’ o ’ )
10 legend ( ’y_1, exact ’ , ’y_2 , exact ’ , ’y_1, ddesd ’ , ’y_2, ddesd ’ )
11 t i t l e ( ’D1 problem of Enright and Hayashi ’ )
12
求解结果如图 (1.14) 所示
图 1.14: 非固定时滞的时滞微分方程
(7) 求解下面的中性时滞微分方程：
y
′
(t) = 1 + y(t) − 2y(t/2)
2
− y
′
(t − π)
0 ≤ t ≤ π
其中：y(t) 的历史信息是：t ⩽ 0 时，y(t) = cos(t)。将下面的程序保存为函数，以便后面调用。
1 f u n c t i o n yp = ddefun ( t , y , ydel , ypdel )
2 yp = 1 + y − 2* ydel ^2 − ypdel ;
3 end
4 function dy = dely ( t , y )
5 dy = t /2;
6 end
7 function dyp = delyp ( t , y )
http://www.ma-xy.com 17 http://www.ma-xy.com
http://www.ma-xy.com
1.4 常微分方程稳定性 第一章 常微分方程
8 dyp = t−p i ;
9 end
10 function y = hi s to r y ( t )
11 y = cos ( t ) ;
12 end
13
求解程序如下
1 tspan = [ 0 pi ] ;
2 s o l = ddensd ( @ddefun , @dely , @delyp , @history , tspan ) ;
3 tn = l i n s p ac e (0 , p i ) ;
4 yn = deval ( s o l , tn ) ;
5 p l o t (tn , yn) ;
6 xlim ( [ 0 p i ] ) ;
7 ylim ([ −1.2 1 . 2 ] ) ;
8
求解结果如图 (1.15) 所示
图 1.15: 中性时滞微分方程
1.4 常微分方程稳定性
微分方程 (组) 的解析法 (初等方法) 经历了很漫长的研究。19 世纪中下，刘维尔的工作表
明：绝大多数微分方程 (组) 不能用初等方法进行求解。前面，我们谈到微分方程 (组) 的解存在
唯一性及其解析方法。但刘表明，绝大多数微分方程求不出解析解。为此，我们有必要研究一下
解 (x(t)) 的其它性质。法国数学家庞加莱 (1852-1912) 发展了解的定性理论，俄罗斯数学家李雅
普诺夫 (1857-1918) 创立了解的稳定性理论。二者的共同特点是在不求解方程 (组) 的情况下，直
接依据方程本身的结构特点来研究定解的性质。
下面，我们来介绍微分方程组的稳定性。考虑下面的微分方程组
˙
x = F (t, x)
其中：F 对 x ∈ D ⊂ R
n
和 t ∈ R 连续，并且对 x 满足利普西茨条件，即解存在唯一。
http://www.ma-xy.com 18 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.4 常微分方程稳定性
系统 (方程组) 中可以有一个变量或者多个变量。稳定性问题的实质是考察系统 x(t)(为简便，
将 x(t) 记为 x(t)) 由初始状态 x
0
扰动 x
0
+ ∆x 引起的受扰运动 x(t; t
0
, x
0
) 能否趋近于或者返回
平衡状态。
平衡状态是指：不受外界干扰的情况下，系统状态 x(t) 不随时间 t 变化而变化；稳定状态是指：
系统状态 x(t) 不随时间 t 变化而变化。由于我们并不考虑外界干扰，故平衡态为 ˙x(t) =
dx
dt
= 0 。
若已知 F (x, t) 的函数形式，我们可以用 F (x, t) = 0 来求平衡态，不妨设平衡态为 x
e
。
再用数学语言描述一下稳定性：设方程对初值 (t
0
, x
0
) 存在唯一解 x
0
(t; t
0
, x
0
)，而对其它初值
(t
0
, x
1
) 存在唯一解 x
1
(t; t
0
, x
1
)。稳定性是指：当 ||x
1
−x
0
|| 很小的时候，||x
1
(t; t
0
, x
1
)−x
0
(t; t
0
, x
0
)||
是否也很小。如果从系统运动图像来看，可以表现为 x(t) 对初值的敏感性问题：如系统稳定性
示意图如图 (1.16) 所示
图 1.16: 系统稳定性示意图
1.4.1 李雅普诺夫意义下的稳定
1. 李雅普稳定是指：初始条件 x
0
在平衡状态 x
e
附近的任意一个初始状态的运动轨迹
x(t; t
0
, x
0
) 均维持在平衡状态 x
e
的附近。
设系统的初始状态为 (t
0
, x
0
)，平衡态 x
e
，初始状态的运动轨迹为 x(t; t
0
, x
0
)。如果对于每个
实数 ε > 0，都存在另一个实数 δ(ε, t
0
) > 0 ，x
0
位于以 x
e
为球心，δ 为半径的闭球域 S(δ) 内，
即 x
0
∈ S(δ)，或者写为
||x
0
− x
e
|| ≤ δ(ε, t
0
) t = t
0
由球域内的任意初始状态 x
0
出发的运动轨迹 x(t; t
0
, x
0
)，在 t → ∞ 时，都位于以 x
e
为球心，以
∀ε > 0 为半径的球域 S(ε) 内，即
||x(t; x
0
, t
0
) − x
e
|| ≤ ε t ≥ t
0
则 x
e
为李雅普稳定，如图 (1.17) 所示，轨迹 x(t) 围绕 x
e
波动且不会逃出 S(ε)。
图
1.17:
李雅普稳定示意图
http://www.ma-xy.com 19 http://www.ma-xy.com
http://www.ma-xy.com
1.4 常微分方程稳定性 第一章 常微分方程
2. 李雅普渐进稳定是指：从域 S(δ) 出发的轨迹 x(t) 不仅不会超出 S(ε)，而且当 t → ∞ 时，
x(t) 收敛于 x
e
。即李雅普渐进稳定首先要求是李雅普稳定，然后，t → ∞ 时，有
lim
t→∞
||x(t; x
0
, t
0
) − x
e
|| → 0
若 δ 与 t
0
无关，且上式的极限过程与 t
0
无关，则称平衡状态 x
e
是一致渐进稳定，如图
(1.18) 所示
图 1.18: 李雅普渐进稳定示意图
3. 李雅普大规模渐进稳定是指：将李雅普渐进稳定的初始条件扩展到整个状态空间时的渐进
稳定，换句话说，就是无论什么样的初始条件，都会稳定。对 ∀x
0
∈ S(δ), δ → ∞, S(δ ) → ∞，都
有
lim
t→∞
||x(t; x
0
, t
0
) − x
e
|| → 0
4. 李雅普不稳定是指：如果对某个实数 ε > 0，和任意实数 δ > 0，不管 ε, δ 有多小，在 S(δ )
内，都存在一个初始状态 x
0
，x
0
出发的轨迹都超出 S(ε)，则平衡状态 x
e
称为李雅普不稳定。
1.4.2 李雅普第二方法 (直接法)
上面，我们介绍了稳定性的概念，但是根据定义来判别 x
e
的稳定性是非常局限的，下面，我
们介绍用李雅普第二方法来判断系统的稳定性。
李雅普函数 V 虚构一个能量函数 V ，V 一般与 x, t 有关，记为 V (x, t)，注意：x 是 t 的函数
x(t)。为了简单，不妨假设 V 不是时变的，即 V 不显含 t，记为 V (x) := V (x(t))。
V (x) 为 n 维矢量 x 定义的标量函数，x ∈ Ω，且在 x = 0 时，V (x) = 0。我们先来定义能
量函数 V 的正定性：
1、如果 V (x) > 0，则 V 为正定。例如：V (x) = x
2
1
+ 2x
2
2
，这里，系统共有两个变量 x
1
, x
2
；
2、如果 V (x) ≥ 0，则 V 为半正定；
3、如果 V (x) < 0，则 V 为负定；
4、如果 V (x) ≤ 0，则 V 为半负定；
5、如果 V (x) > 0 或者 V (x) < 0，则 V 为不定；
上面，我们简单介绍了李雅普函数，下面，就用李雅普函数 V (x) 来判定系统的稳定性。设
系统状态方程为
˙x = F (x, t)
http://www.ma-xy.com 20 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
系统的平衡状态为 F (x, t) = 0，即 x
e
。我们有如下李雅普稳定判定方法：
1. 若
˙
V =
dV
dt
≥ 0，则 x
e
为李雅普稳定。
2. 若
˙
V =
dV
dt
< 0，或者写为
˙
V ≤ 0 和 ∀x ∈ Ω，
˙
V (x(t; x
0
, t
0
)) 不恒等于 0，则 x
e
为李雅普
渐进稳定。
3. 若
˙
V < 0，并且，当 ||x|| → ∞ 时，V (x) → ∞，则 x
e
为李雅普大范围渐进稳定。
4. 若
˙
V > 0，则 x
e
为李雅普不稳定。
上述判据是判定系统是否稳定的充分条件。接下来我们应该思考的是：如何选取 (构造) 函
数 V (x)？其实，V (x) 的选取并不唯一，并且很主观，很方便，没有统一的选取方法。在非线性
系统下，可以依据雅可比矩阵法和变量梯度法进行选取。
1.5 混沌理论 Chaos
1.5.1 混沌系统的引入
在介绍混沌理论之前，我们先来看两个例子：logistic 映射和 Lorenz 方程。(1)logistic 映射
为
x
t+1
= λx
t
(1 − x
t
)
毫无疑问，这是一个含参 λ 的一维动力系统 (离散常微分方程)，在给定初值条件 t
0
= 1 ，x(1) = x
0
时，会有它的运动轨迹 x
t
。下面，我们来研究一下运动轨迹的几种情况，如图 (1.19) 所示
(a) logistic 映射对初值的敏感性 (b) logistic 映射对参数的敏感性
图 1.19: logistic 映射对参数和初值的敏感性
上面的前 3 张图 (1.19a) 是为了说明 logistic 映射对参数 λ 的敏感性，后两张图 (1.19b) 是
为了说明 logistic 映射对初始条件 x
0
的敏感性。由上面这个例子我们可以看出，虽然已经确定
了运动方程 ˙x = f(x, t)，但是其运动轨迹 (方程的解)x
t
却对初始条件及参数极为敏感，其轨迹
曲线表现出不可预测，类似于随机性运动。虽然依据 ˙x = f (x, t) 以及 x(t
0
) = x
0
，我们能够推算
出来未来任意时刻 t 的运动状态 x
t
，但初始数据测量的不准确导致了整体测量的误差极大，甚至
不可预测。20 世纪 70 年代后的研究表明：大量非线性系统尽管状态方程是确定的，却普遍存在
初值敏感问题，运动状态类似于随机运动。下面的工作就是简单研究一下这个初值敏感的确定系
统 (混沌系统)。
http://www.ma-xy.com 21 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
(2) 在介绍混沌 (chaos) 之前，我们再来看一个著名的初值敏感的例子 - Lorenz 方程。美国
气象学家 Lorenz 在研究大气对流时，建立了 Lorenz 方程，下面用一个简单的模拟模型来导出
Larenz 方程。
图 1.20: 大气对流模拟示意图
模拟模型的示意图如图 (1.20)。一个环形管内装满水，从下方点火加热，管内的水发生对流，
设水是不可压缩的，对流时各处的瞬时速度是一样的，速度为 x(θ, t)。由于不可以压缩，我们有
x 与 θ 无关，即写为 x(t)，令管的半径为 1，管的截面积为 1 单位，ν 为粘度，−νx 为流动阻力，
再设管子上的温度为 T ，温度按高度线分布，即
T
1
= T
0
+ ∆T cos θ
其中：T
0
与 ∆T 为常数。
管内流体的温度 T 为
T = T
1
+ y(t) sin θ − z(t) cos θ
其中：z(t) 为管与液体间温度差的对称部分，y(t) 为温差左右不对称部分，图 (1.20) 中左半部分
液体较凉，右半部份液体较热，于是有
dx
dt
= −νx + k y
考虑一段管子 dθ 的热平衡，流入与流出的液体的温差为
∂T
∂θ
dθ，带走的热量为 κx
∂T
∂θ
dθ，管
壁传给流体的热量为 (T
1
− T )dθ，流体随时间吸收的热量为
∂T
∂t
dθ。于是，由热能守恒有
∂T
∂t
= −x
∂T
∂θ
+ κ(T
1
− T )
即
˙y(t) sin θ − ˙z(t) cos θ = −x[−∆T sin θ + y cos θ + z sin θ] + κ(−y sin θ + z cos θ)
令 sin θ cos θ 前的系数左右相等，有
˙y = ∆T x − xz − κy
˙
z
=
xy
−
κz
不妨令
x = κx, y =
νκ
k
y, z =
νκ
k
z
http://www.ma-xy.com 22 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
再令
t =
1
κ
, ρ =
∆T k
κν
, σ =
ν
κ
由此，可以得到 Lorenz 方程
˙x = σ(y − x)
˙y = x(ρ − z) − y
˙z = xy − βz
我们设置参数值为 σ = 10，ρ = 28，β =
8
3
，x
0
= 12，y
0
= 2，z
0
= 9 计算得到它们在三维
空间上的曲线，以及 xy, yz, xz 面的投影和 x, y, z 变量各自的运动轨迹，如下面图 (1.21) 所示
(a) Lorenz 一维轨迹线 (b) Lorenz 二维相面
(c) Lorenz 三维相空间
图 1.21: Lorenz 仿真图
改变系统的初始值，会得到另一个完全不一样的结果：洛伦兹方程的解对初始值十分敏感，
现对 y 的初始值稍加修改，将 2 改为 2.01 和 1.99，让后求解 z 的数值解，如图 (1.22) 所示，随
着时间的推移，三条曲线的吻合程度越来越差，差距越来越大，变化也越来越不明显，成为混沌
状态。
http://www.ma-xy.com 23 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
图 1.22: Lorenz 初值敏感性实验
注：上面轨迹中的不完全自我重复，轨迹不相交缺永不停止转动的双螺旋线即为“Lorenz 吸
引子”，在 Lorenz 方程之后，法国天文学家 Henon 提出了 Henon 映射方程。
1.5.2 混沌的基本理论
混沌现象 混沌是指确定的动力系统因为初值敏感而表现出的不可测的，类似随机的运动。哈密
顿原理将动力学系统划分为可积与不可积两部分，可积系统的典型行为是周期运动和准周期运
动，而混沌是不可积系统的典型行为。俄国数学家 Lyapunov 发明了运动系统的稳定性，同时，还
提出了描述混沌运动的重要指标量 - Lyapunov 指数，该指数是系统稳定性的关键。庞加莱在研
究三体问题时发现了混沌，并用拓扑学的方法来解决问题。现代动力学系统理论的稳定性理论、
分岔理论、奇点理论和吸引子等皆源于庞的研究。19 世纪末，庞在研究三体运动时，发现了三体
系统的不可积和轨道的复杂性，直到 20 世纪 60 年代初，三位数学家 A. 柯尔莫格洛夫、V. 阿诺
尔德和 J. 莫赛证明了 KAM 定理后，才从一定意义上回答了这个问题。
20 世纪 50-60 年代为混沌理论发展的初期 (Lorenz 系统)，20 世纪 70 年代是快速发展时期，
混沌和吸引子等概念被提出，20 世纪 80 年代混沌系统的定量研究成为热点。1980 年，数学家
Mandelbrot 给出了第一张混沌图像，Mandelbrot 集被公认为是混沌的一种标志。混沌系统的性
质量：分数维、lyapunov 指数、Kolmorov 熵等的出现使混沌理论进入应用阶段。
混沌的定义 混沌至今还没有一个统一的定义。下面，我们给出两个影响广泛的定义：Li-Yorke
和 Devaney 的混沌定义。
1975 年，马里兰大学数学家 Yorke 和他的学生 Li(李天岩) 在论文《Period three means
chaos》中给出了混沌的第一个数学定义：
定义 (Li 混沌) 设 f 是 [a, b] 上的连续自映射，其点映射形式为 x
t+1
= f(x
t
, λ)，x
n
∈ [a, b]，
f 为混沌映射 (混沌方程)，如果
(1)f 的周期点的周期无上界，即 f 存在一切周期的周期点。
(2) 存在不可数的子集 s ⊂ [a, b]，s 中无周期点，且满足
¬∀x, y ∈ s, x = y，有
lim
n→∞
inf |f
n
(x) − f
n
(y)| = 0
http://www.ma-xy.com 24 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
∀x, y ∈ s, x = y，有
lim
n→∞
sup |f
n
(x) − f
n
(y)| > 0
®∀x ∈ s，f 的任意周期点 p，有
lim
n→∞
sup |f
n
(x) − f
n
(p)| > 0
其中：f
n
(x) = f (f
n
− 1(x)) = f (f (. . . f(x)))。
Li-Yorke 的混沌定义中的¬和表明子集 s 中的点 x, y 相当分散又相当几种，®表明子集 s
不会趋于任意周期点。
定理 (Li-Yorke 定理) 设 f 在 [a, b] 上连续自映射，若 f 有 3 这个周期点，则对于任意正
整数 n，f 有 n 周期点。
由 Li-Yorke 定理和 Li-Yorke 混沌的定义，我们可以看到：如果 f 周期点 3，那么 f 有任意
整数周期点，故 f 混沌。
另一个著名的混沌定义是 Devaney 的混沌定义：
定义 (Devaney 混沌) 设 D 是紧度量空间 ([a, b])，f 为 D → D 的连续映射，并且满足：
¬f 对初值敏感依赖，即 ∃δ > 0, ∀ε > 0, ∀x ∈ D，在 x 的 ε 邻域内，∃y, n 使得
d[f
n
(x), f
n
(y)] > δ
f 有拓扑传递性，对 D 上任意开集 X, Y ，∃k > 0，f
k
(X) ∩ Y = ϕ。
®的周期点在 D 内稠密。
则 f 为 Devaney 意义下的混沌。
条件¬意味着无论初值依赖多近，在 f 的多次作用之后，二者之间的距离 d 都会超出 δ。条
件意味着任意点的邻域在 f 的多次作用下将遍及整个度量空间 D，这说明 f 不可能分解为两
个在 f 下互不影响的子系统。条件®意味着混沌系统存在着规律性成分。
上面，我们给出了混沌的定义，在给出定义之后，接下来我们应该考虑的内容是：
1. 如何判定一个系统是混沌的？
2. 混沌系统有哪些特征 (性质)？
3. 我们可以用混沌系统来做什么？
对于 1 和 2，其实二者之间是具有联系的，因为我们可以根据混沌系统的特性来判断一个系
统是否是混沌系统。为此，我们有必要先来了解一下混沌系统的性质。
http://www.ma-xy.com 25 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
1.5.3 混沌的基本性质
吸引子 观察 Lorenz 系统的状态轨线图 (1.21)，会发现状态的轨线好像被什么吸住了，但又在
两个“蝴蝶翅膀”上蹦来蹦去。没错，能“吸住”轨线的那个东西 (set) 就是所谓的吸引子了，它
使得一切运动趋向吸引子，而一旦到达吸引子内，运动就开始互相排斥。一般来讲，把 t → ∞ 时
状态的“归宿”称为吸引子 (吸引子是空间中的一个点，或者一个点集)。通常，吸引子被分为两
类：1、平凡吸引子；2、奇异吸引子。而平凡吸引子又分为：1. 不动点吸引子 (平衡运动)；2. 极
限环 (周期运动)；3. 整数维环面 (准周期运动)。吸引子是相空间中的一个集合 (set)。下面给出
吸引子的数学定义：
设 A 是映射 f 的不变集，即 lim
t→∞
f
t
(A)，f (A) = A，其中：f
t
表示 f 作用 t 次。若对 A
的每个 ε 邻域 S(ε) = {x|d(x, A) < ε}，存在 A 的一个 δ 邻域 S(δ)，使得 ∀p ∈ S
δ
，有 f
t
(p) ∈ S
ε
，
t ∈ N，则 A 为稳定集。
若 A 为稳定集，又 ∃a > 0，∀Q ∈ S
δ
，在 t → ∞ 时，d(f
t
(Q), A) → 0，则 A 为渐进稳定
集。若 A 为闭的渐进稳定不变集，则 A 为吸引集；吸引集中含稠密轨线时，即 ∃s ∈ A，使得
{f
t
(s), t = 1, 2, . . . } 的闭包等于 A，则 A 称为吸引子。
注：1. 吸引子是一个集合。
2. 2 维时，相空间即为相平面。
3. 吸引子是 (相) 空间中的一个集合。
4. 相空间一般来说就是状态空间 Ω = {x(t)}。
相空间 Phase Space 因为后面要讨论相空间重构的问题，所以这里简单介绍一下相空间。相：系
统的具体某一状态 x(t)。相空间即为状态空间，运动系统所有可能的状态取值的集合 Ω = {x(t)}。
我们为什么要讨论吸引子？吸引子其实就代表了系统，如果系统具有奇异吸引子，那么这个
系统就是混沌的。研究发现，奇异吸引子具有分数维 (与整数维相对)、正的 Lyapunov 指数和正
的 Kolmogrov 熵等特性，这些性质是平凡吸引子不具备的。由此，我们可以通过计算吸引子的
特征量 (分数维、Lyapunov 指数和 Kolmogrov 熵等) 来判断一个系统是否有奇异吸引子，从而
判断一个系统是否是混沌系统。下面，我们来介绍吸引子的 3 种特征量，以及它们在相空间重构
基础上的数值解法。
吸引子的 2 个特征量 1. 描述邻近轨道发散率的 Lyapunov 指数或者最大 Lyapunov 指数谱。我
们考虑一个一维映射 x(t + 1) = f (x(t))，假设初始条件 x(t
0
) 附近有一点 x(t
0
) + δx(t
0
)，则经过
n 次迭代后，有
δx(t
n
+ 1) = δx(t
n
)f
′
(x(t
n
))
其中：t
0
, t
n
分别为初始时刻和现在时刻
|δx(t
n
)| = |δx(t
0
)|
n−1
i=0
|f
′
(x(t
i
))| =: |δx(t
0
)|e
λt
n
http://www.ma-xy.com 26 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
将上式后面的 λ 分离出来，有
λ
(
x
0
) =
lim
t
n
→∞
1
t
n
n−1
i=1
ln
|
f
′
(
x
(
t
i
))
|
这里的 λ(x
0
) 即为 Lyapunov 指数。当 λ(x
0
) < 0 时，系统有稳定不动点，当 λ(x
0
) = 0 时，对
应分岔点或周期解，当 λ(x
0
) > 0 时，系统为混沌系统。
相似的方法，我们可以得到 n 维映射的 λ
i
(i = 1, 2, . . . , n)，当 max(λ
i
) > 0 时，系统为混沌
系统。
2. 反应信息产生频率的 Kolmogorov 熵。谈到熵，我们自然想到信息论 (概率论) 中关于熵
的定义。考察 n 维动力系统 ˙x = F (x, t)，F = (f
1
, f
2
, . . . , f
n
)
T
，系统在吸引子的轨道 x(t) 把相
空间分割成边长为 r 的 n 维超立方体，每隔时间 ∆t 去测量系统的变化。设 p(i
1
, i
2
, . . . , i
n
) 是
x(r = ∆t) 在超立方体 i
1
、x(t = 2∆t) 在超立方体 i
2
，. . . ，x(r = n∆t) 在超立方体 i
n
中的联合
概率分布，下面简写为 p(·)，则 Kolmogorov 熵为
K
1
= − lim
∆t→0
lim
r→0
lim
n→∞
i
n
p(·) log
2
p(·)
q 阶 Renyi 熵定义为
K
q
= − lim
∆t→0
lim
r→0
lim
n→∞
1
n∆t
1
q − 1
log
2
i
n
p
q
(·)
其中：p(·) 为联合概率分布。显然，我们有
lim
q→ 1
K
q
= K
1
1.5.4 混沌的数值方法
前面，我们给出 Lyapunov 指数 (谱) 和 kolmogorov 熵的定义，但这仅仅是解析方向上的。
对于实际获得的系统状态数据 data，如何在数值上求解吸引子的特征量是我们要重点解决的。下
面，我们将介绍这些特征量的数值求解方法，并基于此，来判定系统是否混沌。因为关联维度、L
指数和 K 熵的数值解法都是基于相空间重构的，所以，我们先来介绍观测数据 data 的相空间重
构。
相空间重构 考虑如下系统的实际测量数据，
˙x = F (x, t)
我们设系统实际测量的状态数据为 D = {x
t
}
N
t=1
(姑且先考虑一维系统 x
t
∈ R)。Packard 等人认
为从一个系统的状态数据 D 可以重构出系统的相空间 (状态空间)，因为 D 本身蕴含了参与动力
系统的全部变量和有关信息，通过将 D 在某些固定时间的延迟点上的观测量看成新的坐标，由
它们共同确定多维空间中的一个点。F.Takens 和 R.Mane 给出了相空间重构的延迟嵌入定理：在
无噪声的情况下，引入两个参数 m, τ ，将观测到的 D 从 t 时刻开始，每隔时间 τ 取个数，连续
取 m 次，形成一个子向量
x
t
= {x
t
, x
t−τ
, . . . , x
t−(m−1)τ
} t = N, N − 1, . . . , (m − 1)τ + 1
http://www.ma-xy.com 27 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
我们令 N
0
= (m − 1)τ + 1，t = N
0
, N
0
+ 1, . . . , N 可以知道，子向量 x
t
是 m 维向量，并且子
向量集 {x
t
}
N
N
0
形成一个 m 维空间，只要 m ≥ 2d + 1，动力系统的几何结构就完全打开，其中：
d 为吸引子的维数，τ 为延迟时间，m 为嵌入维数。状态空间 Φ ⊂ R
m
中吸引子的几何特征与原
动力系统等价，重构的相空间和原系统是微分同胚的。T.Sauer 将延迟嵌入定理扩展到有噪声的
情况。
嵌入维数
嵌入维数 m 的数值确定方法一般有：1. 饱和关联维数 G-P 方法 2. 伪邻近点法 (虚假邻近
点)3. Cao 方法 4. C-C 方法。
虚假邻近点法 F.Takens 和 R.Mane 等证明，当 m ≥ 2d + 1 时，系统的几何结构可以完全打开。
我们不妨逐步增大 m，观察什么时候某些几何结构 (关联积分、关联维数 d 和 L 指数等) 不再随
m 的增加而改变，就停止。虚假邻近点法假设 m 与 τ 无关，基本思想是：当嵌入维数从 m 变
为 m + 1 时，考虑相点 x
n
中哪些是真实的邻近点，哪些是虚假的邻近点，当不存在虚假邻近点
时，几何结构被完全打开。那么，我们要考虑的问题是：¬如何定义虚假临近点?对于重构的相
空间 Φ = {x
n
}
N
n=N
0
中的点，我们不可能要求任意一个相点 x
n
都不存在任何虚假临近点，那么，
我们如何设置虚假邻近点的比率呢？并且注意，在该算法中要求延迟时间 τ 已知。
Step1. 初始化。
初始化系统观测数据 D = {x
t
}
N
t=1
，嵌入维数 m，延迟时间 τ ，最大嵌入维数 m
max
。
Step2. 重构相空间。
依据现在的嵌入维数 m，延迟时间 τ，重构相空间 Φ = {x
n
}
N
n=N
0
。我们记此时的相空间为
{x
n
}
m
。
Step3. 对相空间中所有相点 x
n
计算最邻近点。
对每个 x
n
∈ R
m
，找到与之相邻的最邻近点 x
ηn
，并记录二者之间
d
m
n
= ||x
ηn
− x
n
||
m
2
Step4. 再计算临近点距离。
令 m = m + 1，重构相空间，再次计算所有相点的最邻近点距离
d
m+1
n
= ||x
ηn
− x
n
||
m+1
2
一定要注意，这里不是重新求最邻近点，而是求 x
ηn
在 m + 1 情况下的距离。
Step5. 判断 x
ηn
是否为 x
n
的虚假邻近点。
一般的，当 d
m+1
n
比
m
n
大很多时，判定 x
ηn
是 x
n
的虚假邻近点。为此，我们需要来定义
“大很多”是多少。一个简单的想法是使用阈值，我们给出下面两种判别规则。
判别规则 1：
(d
m+1
n
)
2
− (d
m
n
)
2
(d
m
n
)
2
≥ R
tol
http://www.ma-xy.com 28 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
其中：R
tol
为阈值，一般取值 10-15 之间。
判别规则 2：
d
m+1
n
1
N
N
k=1
(x
k
−
1
N
N
k=1
x
k
)
2
≥ A
tol
其中：A
tol
为阈值，一般取值 2。
上面的两个规则可以一起使用，也可以分开使用。当满足规则时，判定 x
ηn
是虚假邻近点。
Step6. 计数并计算比例 ratio。
记录所有相点 x
n
的含虚假临近点 x
ηn
的数量 C
m
，计算如下比例
ratio =
C
m
N − N
0
Step7. 终止条件。
当 ratio < 5% 或者 C
m
= C
m+1
或者达到最大嵌入维度 m
max
时，终止并输出 m；否则，
置 m = m + 1，返回 Step2.
改进的虚假邻近点法-Cao 方法 从几何观点看，虚假临近点法是一种较好的方法，但在判断虚
假邻点时，阈值的不同会导致结果的不同，这使得这种方法有较大的主观性，为此，我们给出下
面的改进虚假邻近点法-Cao 方法。
Step1. 初始化。
初始化系统观测数据 D = {x
t
}
N
t=1
，嵌入维数 m，延迟时间 τ ，最大嵌入维数 m
max
。
Step2. 重构相空间。
依据现在的嵌入维数 m，延迟时间 τ，重构相空间 Φ = {x
n
}
N
n=N
0
。我们记此时的相空间为
{x
n
}
m
。
Step3. 对相空间中所有相点 x
n
计算最邻近点。
对每个
x
n
∈
R
m
，找到与之相邻的最邻近点 x
ηn
，并记录二者之间
d
m
n
= ||x
ηn
− x
n
||
m
∞
Step4. 再计算临近点距离。
令 m = m + 1，重构相空间，再次计算所有相点的最邻近点
d
m+1
n
= ||x
ηn
− x
n
||
m+1
∞
一定要注意，这里不是重新求最邻近点，而是求 x
ηn
在 m + 1 情况下的距离。
Step5. 计算 a(n, m)。
对每个相点 x
n
∈ R
m
计算 a(n, m)。
a(n, m) =
d
m+1
n
d
m
n
Step6. 计算 E(m), E
∗
(m)。
http://www.ma-xy.com 29 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
计算所有相点 x
n
的 a(n, m) 的均值
E(m) =
1
N − N
0
− 1
N
n=N
0
a(n, m)
E
∗
(m) =
1
N − N
0
− 1
N
n=N
0
|x
ηn −mτ
− x
n−mτ
|
Step7. 计算 E
1
(m), E
2
(m)。
E
1
(m) =
E(m + 1)
E(m)
E
2
(m) =
E
∗
(m + 1)
E
∗
(m)
Step8. 终止条件。
如果 E
1
(m) = E
2
(m) 或者 m > m
max
，终止并输出 m；否则，置 m = m + 1，返回 Step2.
我们可以发现：对于随机时间序列，∀m, E
2
(m) ≡ 1；对于确定时间序列，E
2
(m) 与 m 有关，
E
2
(m) 将不能为常数，所以，改进虚假邻近点法的同时也区分了确定性系统和随机系统。
延迟时间
延迟时间 τ 的数值确定方法一般有：1. 线性相关法 2. 自相关法 3. 平均互信息法 4. C-C 方
法 5. 去偏复自相关法。
自相关函数法 这里，仍然假设嵌入维数 m 和延迟时间 τ 无关。我们的任务是：在给定 m 后如
何求系统的延迟时间 τ 。考虑动力系统的原始观测数据 {x
n
}
N
n=1
，定义系统的平均值 ¯x 为
E(x) = ¯x =
1
N
N
n=1
x
n
给定 τ 的情况下，定义系统的自相关函数 C
l
(τ) 为
C
l
(τ) =
1
N
N
n=1
(x
n+τ
− ¯x)(x
n
− ¯x)
1
N
N
n=1
(x
n
− ¯x)
2
=
E[(x
n+τ
− ¯x)(x
n
− ¯x)]
E(x
n
− ¯x)
2
我们取 C
l
(τ) 第一次为 0 时的 τ 为系统的延迟时间。其实，这里的均值和自相关函数是来自
统计中的定义，如果你学过时间序列。自相关虽然易于计算 (matlab 函数命令为 autocorr)，但
是它只考虑了 {x
n
}
N
n=1
中的线性关系，这是一种局限。当然，我们除了可以选择自相关函数来求
解 τ 之外，还可以选择其它的相关函数来求解。
http://www.ma-xy.com 30 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
平均互信息 考虑动力系统的原始观测数据 {x
n
}
N
n=1
，定义此“时间序列”的平均互信息函数为
I(τ) =
N
i=1
p(x
n
, x
n+τ
) log
2
p(x
n
, x
n+τ
)
p(x
n
)p(x
n+τ
)
其中：p(x
n
, x
n+τ
)，p(x
n
) 和 p(x
n+τ
)w 为概率，可以通过密度估计进行计算。
我们选取 I(τ ) 的第一个最小点时的 τ 为延迟时间，此时产生的冗余最小，独立性最大。平
均互信息是信息论中的概念，我们当然还可以选取其它的概念来求 τ 。注意，由于 τ 的求解不涉
及 m，所以我们先求 τ，再用虚假临近点等方法求 m，这里的 τ 和 m 无关，所以可以分开求解。
后面我们将介绍一种 C-C 方法，这种方法假设 τ 和 m 有关。
C-C 方法 在研究吸引子的动力学特性的过程中，由 Grassbeerger 和 Procaccia 提出的相关维
的概念经常被使用。下面介绍在嵌入时间序列时用到的关联积分 C(r, τ, m, N ) 的概念
C(r, τ, m, N ) =
1
(N − N
0
)(N − N
0
+ 1)
N
i=N
0
N
j=i+1
H(r − ||x
i
− x
j
||) (1.3)
其中：H 是 Hearside 函数，即
H(x) =
0 x ≤ 0
1 x > 0
x
i
, x
j
∈ R
m
为相点。
N
i=N
0
N
j=i
H(r − ||x
i
− x
j
||) 表示除 x
i
本身外，相空间其它相点 x
j
到 x
i
的距离小于 r 的点数。由于距离对称性，所以将 j = i 变为 j = i + 1，然后总数乘以 2。
C-C 算法的研究与函数 (检验统计量)S(m, N, r, t) = C(m, N, r, t) − C
m
(1, N, r, t) 的值有关。
为了研究时间序列的动力学特性, 以及找到合适的延迟时间，我们需要将整个时间序列分为 t 个
子序列。具体方法如下：
给定一个时间序列 (动力系统观测数据){x
i
}
N
i=1
，当 t = 1 时，将序列 {x
i
}
N
i=1
分为 1 个子序
列，也就是序列本身，则
S(m, N, r, t = 1) = C(m, N, r, t = 1) − C
m
(1, N, r, t = 1)
当 t = 2 时，将时间序列分为 2 个子序列
{x
1
, x
3
, . . . , x
N−1
}
{x
2
, x
3
, . . . , x
N
}
则 (两个子序列的检验统计量求平均)
S(m, N, r, 2) =
1
2
{[C
1
(m, N /2, r, 2) − C
m
1
(1, N /2, r, 2)]
+ [C
2
(m, N /2, r, 2) − C
m
2
(1, N /2, r, 2)]}
http://www.ma-xy.com 31 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
更一般的说，将整个时间序列分为 t 个子序列，分别为：
{x
1
, x
1+t
, . . . , x
1+(m−1)t
}
{x
2
, x
2+t
, . . . , x
2+(m−1)t
}
.
.
.
{x
t
, x
2t
, . . . , x
mt
}
则 (t 个子序列的检验统计量求平均)
S(m, N, r, t) =
1
t
t
s=1
{[C
s
(m, N /t, r, t) − C
m
s
(1, N /t, r, t)]}
最后，当 N → ∞ 时，我们有
S(m, r, t) =
1
t
t
s=1
{[C
s
(m, r, t) − C
m
s
(1, r, t)]}
如果时间序列是独立同分布的, 那么对固定的 m, t，当 N → ∞ 时，对 ∀r，均有 S(m, r, t) ≡ 0。
但实际的序列是有限的，并且序列元素间可能相关，实际得到的 S(m, r, t) 一般不为 0。这样，局
部最大时间间隔可以取 S(m, r, t) 穿越零点或者对所有的半径 r 相互差别最小的时间点，因为这
意味着这些点几乎是均匀分布的。我们选择对应值最大和最小半径 r，定义差量为
∆S(m, t) = max S(m, r
j
, t, N) − min S(m, r
j
, t, N)
该量度量了关于半径 r 的最大偏差。所以局部最大时间 t 应该是 S(m, r, t) 的零点和 ∆S(m, t)
的最小值。但是 S(m, r, t) 的零点对所有的 m, r 几乎相等，∆S(m, t) 的最小值对所有的 m 几乎
相等。时间延迟 τ 对应着这些局部最大时间 t 中的一个。
根据 BDS 统计结论，取 m = 2, 3, 4, 5，r
j
= 0 .5 × iσ，i = 1, 2, 3, 4，其中：σ = std(x
i
)，计算
¯
S(t) =
1
16
5
m=2
4
j=1
S(t, m, r
j
)
∆
¯
S(t) =
1
4
5
m=2
∆S(m, t)
S
cor
(t) = ∆
¯
S(t) + |
¯
S(t)|
经过上面的计算，我们可以输出
¯
S(t), ∆
¯
S(t), S
cor
(t) 的图像，找到
¯
S(t) 的第一个零点或者
∆
¯
S(t) 的第一个极小值对应的 t
−
即为延迟时间 τ 。找到 S
cor
(t) 的全局最小值对应的时间 t
+
是
嵌入窗 τ
w
(平均轨道周期) 的最优估计。而对于嵌入窗 τ
w
，有
τ
w
= ( m − 1)τ
于是，由
t
+
= ( m − 1)t
−
http://www.ma-xy.com 32 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
可以求出嵌入维数 m。
在 τ 和 m 确定之后，我们的相空间也就重构好了。下面，在重构相空间的基础上计算系统
的几何不变量：关联积分 C、关联维度 d、L 指数和 K 熵等。
关联维度
下面，我们用 G-P 方法来计算动力系统的关联维度 d。在重构相空间的基础上，我们定义关
联积分 C(r|τ, m, N ) (1.3)，关联积分描述了距离小于 r(r 为外来参数) 的对点数的分布情况，如
果在 r 的某一区间段内，有 C∞r
d
，则 d 为关联维度。
Step1. 初始化。
初始化系统观测数据 D = {x
t
}
N
t=1
嵌入维数 m，延迟时间 τ ，r
min
，r
max
。
Step2. 重构相空间。
依据现在的嵌入维数 m，延迟时间 τ ，重构相空间 Φ = {x
n
}
N
n=N
0
。其中：n = N
0
, N
0
+1, . . . N ，
N
0
= ( m − 1)τ + 1，x
n
∈ R
m
。
Step3. 计算关联积分。
对每一个 r ∈ [r
min
, r
max
]，计算关联积分 (1.3)
C(r|τ, m, N) =
1
(N − N
0
)(N − N
0
+ 1)
N
i=N
0
N
j=i+1
H(r − ||x
i
− x
j
||)
Step4. 绘制 ln C(r) 与 ln r 的图像。
Step5. 确定关联维度 d。
d 为 ln C(r) 与 ln r 的图像中直线部分的斜率。
对于上面介绍的
G-P
方法，如果要改进，很自然想到改进其关联积分
C
(
r
|
τ, m, N
)
的定义。
将其关联积分重新定义为
C(r|τ, m, w, N) =
1
(N − N
0
− w + 1)(N − N
0
− w + 2)
N−w
i=N
0
N
j=i+w
H(r − ||x
i
− x
j
||)
其中：要求 w > τ(2/N)
2/m
，特别地，取 w = τ 即可，w 并没有严格的要求。
上面新定义的关联积分 C(r|τ, m, w, N ) 忽略了相空间中非常靠近的点对关联积分的贡献，只
考虑了 |i − j| ≥ w 的点。
Lyapunov 指数
Lyapunov 指数的数值确定方法一般有：1. Wolf 法、2. Jacobian 法、3. 小数据量法、4. L
指数谱矩阵计算法 BBA(J.P.Eckmann 和 R.Brown)。下面，我们来介绍 Wolf 法和小数据量法。
Wolf 法 Wolf 法是 Wolf 于 1985 年提出的，用于求解最大 Lyapunov 指数的数值方法，也称
轨线法。
Step1. 初始化。τ, m, T
1
, T
2
.
Step2. 重构相空间。
http://www.ma-xy.com 33 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
重构相空间 Φ = {x
n
}
N
n=N
0
。其中：n = N
0
, N
0
+ 1, . . . N ，N
0
= ( m − 1)τ + 1，x
n
∈ R
m
。
Step3. 计算 L
1
。
以初点 x
1
为基点，选取距离 x
1
最邻近的点 x
n1
为端点，记
L
1
= ||x
1
− x
n1
||
Step4. 计算 L
′
1
。
取 T
1
，基点变为 x
1+T
1
，端点变为 x
n1+T
1
，记
L
′
1
= ||x
1
− x
n1+T
1
||
Step5. 计算指数增长率 λ
1
。
λ
1
=
1
T
1
log
2
L
′
1
L
1
Step6. 计算 L
2
。
以点 x
1+T
1
为基点，选取距离 x
1+T
1
最邻近的点 x
n2
为端点，并且使 x
n1+T
1
x
1+T
1
与 x
n2
x
1+T
1
之间的夹角 θ 尽可能小，记
L
2
= ||x
1+T
1
− x
n2
||
Step7. 计算 L
′
2
。
取 T
2
，基点变为 x
1+T
1
+T
2
，端点变为 x
n2+T
2
，记
L
′
2
= ||x
1+T
1
+T
2
− x
n2+T
2
||
Step8. 得到新的指数增长率 λ
2
。
λ
2
=
1
T
2
log
2
L
′
2
L
2
Step9. 重复上述步骤，知道时间步长达到点集 {x, n = 1, 2, . . . , N − (m − 1)τ 的终点，记演化总
步数为 M ，我们用指数增长率 λ
i
, i = 1, . . . , M 的平均值作为最大 L 指数的估计
ˆ
λ(m) =
1
M
M
i=1
λ
i
=
1
M
M
i=1
1
T
i
log
2
L
′
i
L
i
Step10. 终止条件。
如果
ˆ
λ(m + 1) =
ˆ
λ(m)，则终止并输出
ˆ
λ(m)；否则，置 m = m + 1，返回 Step2.
小数据量法 在介绍小数据量法之前，我们先来定义时间序列的平均周期 P 。设系统观测数据为
{x
t
}
N
t=1
，经过 FFT 变换后，其频率为 f
1
, f
2
, . . . , f
N
，其幅值为 A
1
, A
2
, . . . , A
N
，则 {x
t
}
N
t=1
的
平均频率 (功率谱的平均频率) 为
F =
N
i=1
f
i
A
i
N
i=1
A
i
http://www.ma-xy.com 34 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
{x
t
}
N
t=1
的平均周期 p 可以用功率谱的平均频率的倒数估计，即
p =
1
F
=
N
i=1
A
i
N
i=1
f
i
A
i
下面，给出小数据量法的算法流程：
Step1. 初始化。
初始化系统观测数据 D = {x
t
}
N
t=1
，τ, m, p，这里的 p 为 {x
t
}
N
t=1
的平均周期。
Step2. 重构相空间。
重构相空间 Φ = {x
n
}
N
n=N
0
。其中：n = N
0
, N
0
+ 1, . . . N ，N
0
= ( m − 1)τ + 1，x
n
∈ R
m
。
Step3. 找到每个相点 x
n
的最邻近点 x
ηn
，记二者之间的距离为
d
n
= ||x
ηn
− x
n
||
2
为了避免 x
n
与 x
ηn
在同一轨线上，我们要求
|ηn − n| > p
Step4. 设置离散步长 t。
t = 1, 2, . . . , min(M − ηn, M − n)，其中：M 为相空间中点的个数，M = N − N
0
+ 1.
Step5. 计算 d
n
(t)。
对每个相点 x
n
，计算
d
n
(t) = ||x
n+t
− x
ηn +t
||
Step6. 对每个 t，求所有 x
n
的 ln d
n
(t) 的平均值 y(t)
y(t) =
N
n=N
0
(ln d
n
(t))
Step7. 求最大 L 指数的估计。
用 OLS 等方法求 y(t) 的回归直线，其斜率 λ 即为最大 L 指数的估计。
Kolmogorov 熵
Kolmogorov 熵的数值确定方法一般有：G-P 法、STB 法，这里不做介绍。
1.5.5 混沌系统的应用
时间序列建模步骤
混沌现象一般发生在非线性系统下颚确定性系统当中，对于具体的某一个测量数据 {x
n
}
N
n=1
，
x
n
∈ R，我们给出该时间序列的一般建模步骤，如图 (1.23) 所示
http://www.ma-xy.com 35 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
图 1.23: 时间序列建模流程图
从上面的建模流程图中，可以看到，我们需要的工具有：序列降噪方法、非线性检验、平稳
时间序列建模、确定性检验、混沌性检验、相空间重构以及预测方法。前面主要介绍的就是混沌
系统的相空间重构以及混沌性检验方法，下面，我们来介绍非线性检验和确定性检验两种技术，
关于序列降噪方法和平稳时间序列建模 (ARMA 等方法)，在后面的数学建模实战部分会有一定
的介绍。
非线性检验
常见的非线性检验方法有：AAFT、IAAFT、CAATT 和 GCR。下面，我们就来介绍基于
替代数据思想的 GCR 方法。
我们的目标是检验数据来自线性系统还是非线性系统，一种非常好的方法就是替代数据法。
替代数据的基本思想是：原假设 H0：{x
n
}
N
n=1
来自某线性系统族。如果在这个 H0 的基础上，
数据所表现出来的性质与线性系统相违背，那么，就可以拒绝原假设 H0。我们从线性系统族中
生成一些替代数据，如生成 M 个时间序列 {˜x}
i
，i = 1, . . . , M，计算 M 个时间序列的统计量
Q
i
, (i = 1, . . . , M )，统计量是样本的函数，可以是 L 指数、关联维度和 Takens 统计量，这之后
在计算测量数据 {x
n
}
N
n=1
的统计量 Q
0
。最后，判断 Q
0
是否在 Q
i
, (i = 1, . . . , M ) 的分布当中。
上面的替代数据的思想正是假设检验的思想：在原假设 H0 成立时，Q 统计量应该是什么样
的分布 (上面采用生成数据 {˜x}
i
，i = 1, . . . , M 来计算)，而样本的实际统计量为 Q
0
，判断 Q
0
是否落在 Q 分布的拒绝域内。我们要讨论的是：如何从线性系统族中生成一些替代数据？
下面，我们来讨论如何生成 M 个时间序列 {˜x}
i
，i = 1, . . . , M 。如果我们要求 {˜x
n
}
i
与 {x
n
}
来自同一分布，并且有近似的自相关性，我们可以考虑就在原序列 {x
n
} 的基础上，值不变，打乱
顺序即可生成 M 个时间序列 {˜x
n
}
i
，这样 {˜x}
i
来自同一分布。然而，如何才能使 {˜x
n
} 与 {x
n
}
的自相关性相同呢？我们可以将这一要求转化为目标，利用前面介绍的 GA 等智能优化算法的思
http://www.ma-xy.com 36 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
想，生成多个 {˜x
n
}，计算它们的自相关性，进化那些自相关性与 {x
n
} 弱的样本 {˜x
n
}。
min
α
max{|
˜
C(τ) − C(τ )|}
其中：α 为序列 {x
n
} 的一种排列方式，例如 1, 2, 3... 和 3, 2, 4, ...。
˜
C(τ) 为样本 {˜x
n
} 在 τ 下的
相关函数值，C(τ) 为样本 {x
n
} 在 τ 下的相关函数值。
当然，在其它文献当中，可能定义 max{|
˜
C(τ) − C(τ )|},τ = 1, 2, . . . , τ
max
为成本函数 (损
失函数)，这也未尝不可。关于上面的 min max 问题的求解，利用 GA 等算法是简单的，但是在
{x
n
} 的 n 较大时，要考虑计算量的问题。
对于前面介绍的替代数据法，还存留一个问题，那就是检验统计量
Q
的选取及检验。
Q
一
般选取 L 指数、关联维度 d 和 K 熵等系统的几何特征。也可以选取平均意义下的一步预测误差
等。
关于检验统计量 Q 的分布，我们有
|Q
0
−
¯
Q
s
|
σ
s
˙∼ N (0, 1)
其中：Q
0
为原序列 {x
n
} 的检验统计量值，
¯
Q
s
为 M 个替代时间序列 {˜x
n
} 的统计量平均值，σ
s
为统计量方差。上面的分布渐进服从标准正态分布。
确定性检验
如果系统的将来值可以由过去值完全确定，则称系统为确定性系统，例如：x
n+1
= 4 x
n
。我
们可以写出确定性系统的一般形式
x
n+1
= F (x
n
)
这里的传递函数 F 完全已知。但是当确定性系统具有混沌特性时，一个小的测量误差会导致预测
误差的极大变化，导致系统的不可测性。虽然不可测，但是混沌系统与随机系统仍是不同的，因为
随机系统的未来值是完全不可预测的。对此，我们可以这样来做：设系统的观测序列为 {x
n
}
N+l
n=1
，
用前 N 个观测值来重构相空间，并预测后 l 个值，记为 {ˆx
n
}
N+l
N
，将其与 {x
n
}
N+l
N
比较，然后
设定判别规则。
当然，关于预测问题，在 1 步预测时，如果认为 1 次预测有太大的随机性，也可以进行多次
一步预测，然后取多次预测的平均值作为最终预测值。在多步预测时，有两种方案：¬利用 BP
等方法一次进行多步预测；一步一步的预测：将前一步预测值作为已知值，重构相空间，然后
进行 1 步预测。
王海燕的《非线性时间序列分析》一书中还提到过用递归图来判定系统的确定性。下面，我
们直接介绍递归图的画法以及如何判断系统的确定性。
设观测序列为 {x
n
}
N+l
n=1
，在 m, τ 给定之后，重构相空间为 Φ = {x
n
}
N
n=N
0
，相点 x
n
∈ R
m
为系统在 n 时刻的状态。计算 i, j 时刻相点 x
i
, x
j
的距离 d
ij
= ||x
i
− x
j
||，设定距离的阈值为
r，以 i 为横坐标，以 j 为纵坐标，当 d
ij
< r 时，在 (i, j) 处画一个点。观察递归图的对角线两
侧是否有与对角线近似平行的小短线，如果有，则为混沌系统 (确定性系统)；如果没有，则为随
机系统。注意：上面方法中的量和关联积分中的是相似的。
http://www.ma-xy.com 37 http://www.ma-xy.com
http://www.ma-xy.com
1.5 混沌理论 CHAOS 第一章 常微分方程
预测方法
零阶局部预测 对于混沌时间序列 (混沌系统)，我们在相空间重构的基础上来尝试的进行预测。
对于原序列 {x
n
}
N
n=1
我们要预测其下一个状态值 x
N+1
，那么，在相空间重构上，预测 x
N+1
就相
当于预测相点 x
N+1
(相点 x
N+1
中只有最后一个分量 x
N+1
是未知的)。现在，共有 N − N
0
个相
点 x
n
∈ R
m
，我们要根据这些相点预测下一个 x
N+1
。基本思路是：相点 x(N ) 靠近 x(N + 1)，我
们在相空间 Φ 中找到 k 个与 x(N ) 靠近的点 x(N
1
), x(N
2
), . . . , x(N
k
)，那么，x(N
1
+ 1), x(N
2
+
1), . . . , x(N
k
+1) 这 k 个相点也靠近 x(N +1)，于是，我们可以让 k 个相点的平均值作为 x(N +1)
的估计值
ˆ
x(N + 1) =
1
k
k
i=1
x(N
i
+ 1)
上述方法称为零阶局部预测。
加权零阶局部预测 但是，不得不考虑 k 个相点与 x(N ) 接近的程度，这样，就可以设置加权
零阶局部预测，让接近程度大的相点的权重大。用 d
i
表示 x(N
i
) 与 x(N ) 之间的距离，并令
d
min
= min{d
i
}，则加权零阶局部预测为
ˆ
x(N + 1) =
k
i=1
x(N
i
+ 1)e
−l(d
i
−d
min
)
k
i=1
e
−l(d
i
−d
min
)
=
k
i
=1
e
−l(d
i
−d
min
)
k
i=1
e
−l(d
i
−d
min
)
x(N
i
+ 1)
其中：l 为外来参数，具有主观性；权重为
w
i
N
=
e
−l(d
i
−d
min
)
k
i=1
e
−l(d
i
−d
min
)
可以尝试着用局部回归核函数来替代权重，也可以用小波核来代替。由于参数 l 具有主观性，
向晓东等人提出使用
w
i
N
=
d
M
−d
i
d
M
−d
min
k
i=1
d
M
−d
i
d
M
−d
min
来作为权重。
加权 1 阶局部预测 像局部常数核权重回归和局部多项式核权重回归那样，对于零阶局部预测，
我们有 1 阶局部预测和多阶局部预测，下面，我们来看一下加权 1 阶局部预测。
前面，我们直接用 x(N ) 来替代 x(N + 1)，现在，假设 x(N + 1) 与 x(N ) 之间是多项式关
系，姑且先假设为 1 阶多项式 (即线性)，则有
x(N + 1) = a + bx(N )
那么，接下来的工作就是求参数 a, b。在相空间 Φ 中找到 k, (k > m + 1) 个与 x(N) 邻近的相点
x(N
1
), x(N
2
), . . . , x(N
k
)，则有
a + b
x(N
i
) =
x(N
i
+ 1)
http://www.ma-xy.com 38 http://www.ma-xy.com
http://www.ma-xy.com
第一章 常微分方程 1.5 混沌理论 CHAOS
其中：
x(N
i
) = (x(N
1
), x(N
2
), . . . , x(N
k
))
T
。
于是，我们利用 (加权)OLS 等方法可以求解 a, b
a
∗
, b
∗
= arg min
a,b
k
i=1
w
i
[x(N
i
+ 1) − a − bx(N
i
)]
2
最终，我们有
ˆ
x(N + 1) = a
∗
+ b
∗
x(N)
基于最大 L 指数的预测 假设相空间 Φ 的最大 L 指数为 λ
1
，找到相点 x(N) 的最邻近点 x(Nη)
以及 x(N η + 1)，由最大 L 指数的物理意义
e
λ
1
=
||x(N + 1) − x(Nη + 1)||
||x(N) − x(N η)||
有
||x(N + 1) − x(Nη + 1)|| = de
λ
1
由于 x(N + 1) 中只有最后一个分量 x
N+1
是未知的，因此上式可解。
http://www.ma-xy.com 39 http://www.ma-xy.com
