系统仿真技术
第 2章 经典的连续系统仿真建模方法学
剡昌锋 刘军
兰州理工大学机电工程学院
Tel:13369316141(剡） 13893121598（刘）
E_mail,changf_yan@163.com
Lj_gabe@sohu.com
2.1 离散化原理及要求
?问题：数字计算机在数值及时间上的 离散性 -
---被仿真系统数值及时间上的 连续性?
?连续系统的仿真，从本质上：对原连续系统
从时间、数值两个方面 对原系统进行离散化
并选择合适的数值计算方法来 近似 积分运算
? 离散模型 ≈原连续模型？
相似原理
?设系统模型为：,其中 u(t)为输
入变量,y(t)为系统变量；令仿真时间间隔
为 h，离散化后的输入变量为,系统变
量为,其中 表示 t=nh。
?如果,且
?即,
?（对所有 n=0,1,2,… ）
?则可认为两模型等价。
(,,)y f y u t=&
)(? ntu
)(? nty nt
)()(? nn tutu ? )()(? nn tyty ?
?( ) ( ) ( ) 0u n n ne t u t u t= - ?( ) ( ) ( ) 0y n n ne t y t y t= -
u(t)
h
y(t)
-
+
图 2.1 相似原理
原连续模型 (,,)y f y u t=&
仿真模型 ? ? ?(,,)ny f y u t=&
)(? nty
0)( ?ny te
)(? ntu
对仿真建模方法三个基本要求
? （ 1） 稳定性,若原连续系统是稳定的，则离散化
后得到的仿真模型也应是稳定的。
? （ 2） 准确性,有不同的准确性评价准则，最基本
的准则是,
? 绝对误差准则,
? 相对误差准则,
? 其中 ? 规定精度的误差量。
???? )()(?)( nnny tytyte
???? )(? )()(?)(
n
nn
ny ty
tytyte
对仿真建模方法三个基本要求 (续 )
? 3） 快速性,若第 n步计算对应的系统时间间
隔为 计算机由计算需要的时间为 Tn，
若 Tn=hn 称为 实时仿真, Tn?hn称为 超实
时仿真 Tn?hn 称为 亚实时仿真,对应于 离
线仿真
,1 nnn tth ?? ?
数值积分算法
? 对,已知系统变量 y的初始条
件,要求 y随时间变化的过程 ―― 初值问
题
? 计算过程：由初始点 的
? 欧拉法
? 对任意时刻 tn+1
? 截断误差正比于
(,,)y f y u t=&
00()y t y=
00()y t y= )( 00 ytf,
0
0( ) (,)
t
ty t y f t y d t=+ ò
1 1 0 0 0( ) ( )y y t y t f t y= @ + D,
1 1 1( ) ( ) ( )n n n n n n ny y t y t t f t y+ + += @ + -,
h2
f ( t,y )
f ( t
0
,y
o
)
tt
t
0 t 1
数值积分算法 (续）
?梯形法,
?是隐函数形式。预报－ — 欧拉法估计初值，
校正－ — 用梯形法校正,
?校正公式
?预报公式
?反复迭代，直到满足
?经典的数值积分法可分为两类,
?单步法与多步法
)],(),([21)( 1111 ???? ???? nnnnnnn ytfytfhytyy
)],(),([21 )( 11)1( 1 innnnnin ytfytfhyy ???? ???
),()( 1 nnnin ytfhyy ????
??? ??? inin yy 111
2.2 龙格库塔法
? 2.2.1龙格 -库塔法基本原理
?对
?若令,
?则有
? 的数值求解：称作“右端函数”计算问题。
在 附近展开台劳级数，只保留 项，则有,
（ 1）
? ???? 1 ),()()( 1 nnttnn dtytftyty
)( nn tyy ? ? ?? 1 ),(Q n
n
t
tn dtytf
nnnn yyty Q)( 11 ??? ??
nQ
t0 h2
2
0001 0)(2
1),( h
t
f
dt
dy
y
fhytfyy
t?
?
?
? ????
龙格 -库塔法基本原理（续）
?假设这个解可以写成如下形式,
?其中
?对 式右端的函数展成台劳级数，保留 h项，
可得,
?代入，则有,
1 0 1 1 2 2()y y a k a k h= + +
),( 001 ytfk ? k f t b h y b k h2 0 1 0 2 1? ? ?( )，
k2
02 0 0 1 2 1(,) ( ) t
ffk f t y b b k hty抖@ + +抖
01 0 1 0 0 2 0 0 1 2 1(,) [ (,) ( ) ]t
ffy y a h f t y a h f t y b b k hty抖= + + + +抖
龙格 -库塔法基本原理（续）
?将 (2)式与 (1)式进行比较，可得,
?四个未知数 但只有三个方
程，因此有无穷多个解。
?若限定,则
?计算公式,
?其中
a a a b a b1 2 2 1 2 21 1 2 1 2+ = = =,/,/
，，，,2121 bbaa
a a1 2? a a b b1 2 1 212 1? ? ? ?，
)(2 2101 kkhyy ???
)(),( 1002001 hkyhtfkytfk ????,，
龙格 -库塔法基本原理（续）
?若写成一般递推形式，即为,
?其中
?截断误差正比于 h3，称为二阶龙格 -库塔法
（简称 RK-2）。
?截断误差正比于 h5的四阶龙格 --库塔法（简称
RK-4）公式,
?其中,
)(2)( 2111 kkhyyty nnn ???? ??
),(),( 121 hkyhtfkytfk nnnn ????,
)22(6)( 432111 kkkkhyyty nnn ?????? ??
),(1 nn ytfk ? )
22( 12 k
hyhtfk
nn ???,
)22( 23 khyhtfk nn ???, )( 34 hkyhtfk nn ???,
2.2.2龙格 --库塔法的特点
? 1.形式多样性
? 例,非唯一解，可以得到许多种
龙格 --库塔公式,(中点公式 )
? 其中
? 各种龙格 ---库塔法可以写成如下一般形式,
? 其中,
2121 bbaa,，，
hkyy nn 21 ???
),(1 nn ytfk ? )
22( 12 k
hyhtfk
nn ???,
?
?
? ??
s
i
iinn kChyy
1
1
)(
1
1
??
?
???
i
j
jijnini kbhyhatfk, si,，,?21?
龙格 --库塔法的特点 (续）
? 式中各系数满足以下关系
? s称为级数，表示每步计算右端函数 f的最少次数。
可以证明,1阶公式至少要计算一次,2阶公
式 ； …,； 4阶公式 ；依此类推。有时
为了某种特殊需要，可以选择 的计算公式。
a
a b i s
C
i ij
j
i
i
i
s
1
1
1
1
0
2 3
1
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
，,,?
smin ? 2 smin ? 4
minss ?
龙格 --库塔法的特点 (续）
? 2.单步法
? 在计算 时只用到,而不直接用
等项。优点：存储量减小，可以自启动
? 3.可变步长
? 步长 h在整个计算中并不要求固定，可以根
据精度要求改变
? 但是在一步中，为计算若干个系数,则
必须用同一个步长 h。
1?ny ny 21 ?? nn yy,
ik
龙格 --库塔法的特点 (续）
? 4.速度与精度
?四阶方法的 h可以比二阶方法的 h大 10倍，每
步计算量仅比二阶方法大一倍
? 高于四阶的方法由于每步计算量将增加较
多，而精度提高不快。
2.2.3实时龙格－库塔法
?实时仿真,要求仿真模型的运行速度往往与
实际系统运行的速度保持一致。
?一般的数值积分法 难以满足 实时仿真的要求，
这不仅仅是因为由这些方法所得到的模型的
执行速度较慢，而且这些方法的机理不符合
实时仿真的特点。
?考虑系统 ))(( tuyf
dt
dy,?
实时龙格－库塔法 (续）
? RK-2公式如下,
? 一个计算步内分两子步,
? tn时刻：利用当前的 un,yn计算 k1----计算一次右端
函数 f需
? tn+h/2时刻：应计算 k2，尽管此时 yn +1/2已经得到，
但 un +1则无法得到。 (若对 un +1也进行预报 ―― 加大
仿真误差 )。仿真执行延迟 h/2―― 输出要迟后半个
计算步距
),,(
),,(
)(
2
1112
1
211
hkyutfk
yutfk
kk
h
yy
nnn
nnn
nn
??
?
???
??
?
h / 2
实时龙格－库塔法 (续）
RK-2的计算流程
并输出计算 1?ny1?nu采入
1k计算
nt
2
ht
n ?
1??? nn tht
21
ht
n ??
21 ?? ?? nn tht
2k计算 1k计算下一个
实时龙格－库塔法 (续）
?实时 2阶龙格－库塔法,
? tn时刻：应计算 k1，利用当前的 un,yn，需
要 ； tn+h/2时刻，应计算 k2，此时 yn +1/2已
经得到,un +1/2也可得到,k2的计算就不会引
入新的误差。计算一次右端函数 需要,
可实时输出 yn +1。
)
2
,,(
),,(
12/12/12
1
21
k
h
yutfk
yutfk
hkyy
nnn
nnn
nn
??
?
??
??
?
h/2
f h/2
实时龙格－库塔法 (续）
实时 RK-2公式计算流程
)2/( htu n ?采入 并输出计算 1?ny
nt 2htn ? 1??? nn tht
21
ht
n ?? 21 ?? ?? nn tht
1k计算 2k计算 1k计算下一个
2.3 线性多步法
? 2.3.1线性多步法基本原理
?基本原理,利用一个多项式去匹配变量若干
已知值和它们的导数值。
?设,时刻的 和
已知 ;
?预报：由 和 来计算
?校正：若 也已知，由它们来计算
t t tn n n k,,,? ? ?1 1? y y yn n n k,,,? ? ?1 1? ?,?,,?y y yn n n k? ? ?1 1?
y y yn n n k,,,? ? ?1 1? ?,?,,?y y yn n n k? ? ?1 1? y yn k n k? ?,?
kny? ?yn k?
线性多步法 (续 )
? 采用的多项式具有以下形式 (m阶 )
? 其中,是待定系数,
? 在 时刻，,可得到,
? （ 1）
y t d
t t
h
d
y t id
t t
h
id
m i
n k
i
i
i
i
m
i
m
m h i
n k
i
m i
h i
i
i
m
( )
? ( )
?
??
?
?
?
?
? ?
? ?
??
?
?
?
?
? ? ?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
??
? ?
?
?
00
1
1
1
1 1
1
di h tt kn ?? ??
tn k? 0?? y dn k? ? 0
?y dn k h? ? ? 1 1
线性多步法 (续 )
?由 和 确定,
需要 m+1个独立方程。该 m+1个方程可由以
下等式导出,
? （ 2）
y y yn n n k,,,? ? ?1 1? ?,?,,?y y yn n n k? ? ?1 1? ),2,1,0( mid
i ??
y t y d
t t
h
d j
y t y id
t t
h
id j
m n k j i
n k n k j
i
i
i
i
m
i
m
m n k j h i
n k n k j
i
m
i
h i
i
i
m
( )
? ( ) ?
? ?
??
?
?
?
?
? ?
? ? ?
??
?
?
?
?
? ? ?
?
?
?
?
?
?
?
? ?
? ? ?
??
? ?
? ? ?
?
?
?
?
??
? ?
00
1
1
1
1 1
1
线性多步法 (续 )
? 1、预报公式
? 令 m=2k-1，从（ 2）式得到如下方程组,
d d d d yp p p mp n k0 1 2 1? ? ? ? ? ? ??
d d d d yp p p m mp n k0 1 2 22 4 2? ? ? ? ? ? ??
d d d d yp p p m mp n k0 1 2 2 33 3 3? ? ? ? ? ? ??
?
d d md h yp p mp n k1 2 12? ? ? ? ? ? ?? ?
d d m d h yp p m m p n k1 2 1 22 2 2? ? ? ? ? ? ?? ? ?? ?
d d m d h yp p m m p n k1 2 1 32 3 3? ? ? ? ? ? ?? ? ?? ?
?
? 将其写成矩阵形式,（ 4）
? 其中上标 p表示预报。
? 其解为, （ 5）
V d Zp p p?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
??
??
??
??
??
??
?
?
?
?
?
?
n
kn
kn
kn
n
kn
kn
kn
m
k
k
k
k
m
m
m
m
m
m
yh
yh
yh
yh
y
y
y
y
d
d
d
d
d
d
d
d
kmkk
m
m
m
kkkk
?
?
?
?
?
?
?
?
?
??
?
?
?
?
??
?
?
?
3
2
1
3
2
1
2
1
1
2
1
0
12
12
12
32
32
32
3210
3233210
2232210
3210
1
33331
22221
11111
（ 3）
? ?d V Zp p 1 p? ?
其解为,? ?
d V Zp p 1 p? ?
由于 为常数阵，其逆存在,Z向量中的各元素为
已知值，因而 d向量的各元素值可计算得到，从而
由,得到下一时刻的预报值 。
Vp
y d y dn k n k h? ?? ? ?0 1 1,?
缺点：只有 是所需要的，其它元素的计算成
为多余,得不到 与 和
显式表达式
d d0 1,
y yn k n k? ?,? y y yn n n k,,,? ? ?1 1?
?,?,,?y y yn n n k? ? ?1 1?
?定义,（ m+1） ?1的列向量
? (6)
?定义辅助变量 (7)
?此式可改写为 (8)
?向量 的元素可划分为两个组
? (9)
T)0,,0,1( ??0e
? ?y n k? ?? ?e d e V Z0T p 0T p 1 p
? ? ? ?? p T ? ?e V0T p 1
? ?V ep T 0? p ?
?p
? ?? p p p k p p p k p Ta a a b b b? 1 2 1 2,,,,,,? ? ?
?例,k=3，则 (8)式为,
?可计算得到,
1 1 1 0 0 0
1 2 3 1 1 1
1 2 3 2 2 2 2 3
1 2 3 3 3 2 3 3
1 2 3 4 4 2 4 3
1 2 3 5 5 2 5 3
1
0
0
0
0
0
2 2
3 3 2 2
4 4 3 3
5 5 4 4
1
2
3
1
2
3
? ?
? ?
? ?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
a
a
a
b
b
b
p
p
p
p
p
p
? ? ? ?? p p p p p p p Ta a a b b b? ? ? ? ? ?1 2 3 1 2 3 18 9 10 9 18 3,,,,,,,,,,? ?
? 只依赖于 k，即先前 和 的个
数，而与它们的数值无关。这样，可以预
先求解（ 8）式得到。
?从而得到 的显式表达式,
?p yn k j? ? ?yn k j? ?
? ?y n k p T? ? ?e d Z0T p p?
yn k?
y a y h b yn k jp n k j jp n k j
j
k
j
k
? ? ? ? ?
??
? ? ?? ?
11
?例：试推导用 预报 公式
?条件：已知
?由
y y yn k n k n k? ? ? ? ? ?1 1 2,?,? yn k?
? ?V ep T 0? p ?
V p ?
?
?
?
?
?
?
?
?
?
?
1 1 1
0 1 2
0 1 4
y y yn k n k n k? ? ? ? ? ?1 1 2,?,?? ? ? ?? p T p p pa b b?
1 1 2,,
? ?V p T ? p
p
p
p
a
b
b
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1 0 0
1 1 1
1 2 4
1
0
0
1
1
2
a b bp p p1 1 21 3 2 1 2? ? ? ?,/,/
? ?y y h y yn k n k n k n k? ? ? ? ? ? ?? ? ?1 12 1 23 ? ?
2、校正公式
?预报公式 ―― 显式公式，未包括 。
?校正,对该预报值应进行校正，即先预报得
到,然后再用此值推出
?由 和 以及 来预
报,
?可令 m=2k-1，从（ 2）式得到如下方程组,
?yn k?
~yn k? ?yn k?
y y yn n n k,,,? ? ?1 1?,~ kny ? ?,?,,?y y yn n n k? ? ?1 1?
?yn k?
kn
c yd
?? ~0
?将其写成矩阵形式,
d d d d yc c c mc n k0 1 2 1? ? ? ? ? ? ??
d d d d yc c c m mc n k0 1 2 22 4 2? ? ? ? ? ? ??
d d d d yc c c m mc n k0 1 2 2 33 3 3? ? ? ? ? ? ??
?
d d md h yp c mc n k1 2 12? ? ? ? ? ? ?? ?
d d m d h yc c m mc n k1 2 1 22 2 2? ? ? ? ? ? ?? ? ?? ?
d d m d h yc p m mc n k1 2 1 32 3 3? ? ? ? ? ? ?? ? ?? ?
?
V d Zc c c?
校正公式 (续 )
?其中上标 c表示校正，可得
? (14)
? ?d V Zc c 1 c? ?
1 0 0 0 0
1 1 1 1 1
1 2 2 2 2
1 3 3 3 3
1
0 1 2 3
0 1 2 2 3 2 2
0 1 2 3 3 3 3
2 3
2 3
2 3
2 1
2 1
0
1
2
3
1
2
?
?
?
?
?
?
?
?
?
?
?
?
m
m
m
m
m
k
k
k
m
k k k k
m
m
m
d
d
d
d
d
d
d
d
? ? ?
? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ?
? ?
? ?
? ?
? ?
? ?
y
y
y
y
y
h y
h y
h y
n k
n k
n k
n k
n
n k
n k
n k
1
2
3
1
2
3
?
?
?
?
?
?定义,为（ m+1） ?1的列向量，
上标 T表示转置。将 左乘（ 13）式可得,
? （ 15）
?定义 （ 16）
?可改写为 （ 17）
? （ 18）
e 1T ? (,,,,)0 1 0 0?
e1T
? ?? ? ?? ?h y n k? e d e V Z0T c 1T c 1 c
? ? ? ?? c T ? ?e V1T c 1
? ?V ec T 1? c ?
? ?? c c c kc c c kc Ta a a b b b? 0 1 1 2,,,,,,? ? ?
?例,k=3
?同样,只依赖于 k，即先前 和 的个
数，而与它们的数值无关。这样
? （ 19）
?从而 （ 20）
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
??
0
0
0
0
1
0
35255210
34244210
33233210
32222210
111210
000111
3
2
1
3
2
1
445
334
223
2
p
p
p
p
p
p
b
b
b
a
a
a
?c yn k j? ? ?yn k j? ?
? ?? ? ??h y n k c T? e d Z1T c c?
? ( ? )y h a y h b yn k jc n k j jc n k j
j
k
j
k
? ? ? ? ?
??
? ? ? ??1
10
?例：已知,预估,然后用
校正 。
预估
预报公式为
校正
校正公式为
y y yn k n k x k? ? ? ? ? ?1 2 3,,kny ?~
21,,~ ????? knknkn yyy ?yn k?
1 1 1
1 2 4
1 3 9
0
1
2
1
2
3
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ?
? ?
? ?
d
d
d
y
y
y
p
p
p
n k
n k
n k
? ? ? ? ? ?? p T p p pa a a? ? ?1 2 3 3 3 1,,,,
321 33~ ??????? ??? knknknkn yyyy
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
2
1
2
1
0
~
421
111
001
kn
kn
kn
c
c
c
y
y
y
d
d
d
? ? ? ? ? ?2123210,2,,,???? cccTc aaa?
)2~( 2211231 ?????? ????? knknknhkn yyyy?
2.3.2 线性多步法误差分析
? 为了便于分析，对预报公式和校正公式，定义统一的表达式,
? （ 21）
? -------显式预报 显式预报
? 时称为后向差分公式（ BDF）
? 同时均不等于 0时为隐式校正公式
? k称为公式的阶次 。假设变量各时间的 精确值 已经得到，将
其代入（ 21）式，可得,
? （ 22）
? ?
? ?
???? ??
k
i
k
i
ikniikni yhy
0 0
0???
00 ?? ?yn k? 00 ?? yn k?
021 ???? k??? ?
00,??
? ?
? ? ??
???k
i
k
i knikni
ihtyhihty
0 0
)()( ???
线性多步法误差分析（续）
?在 附近，将每个函数展开成泰勒级数,
? （ 23）
?对所有 i(i=0,1,2,…,k), 将（ 23）式代入（ 22）
式，合并同类项，可得
?
? （ 24）
tn k?
???? ????? ???? )()()()( !2 )(!1 2 knihknihknkn tytytyihty
???????? ????? ???? )()()()( !2 )(!1 2 knihknihknkn tytytyihty
????? ????? ???? )()()()( )(2210 knmmmknknkn tyChtyChtyhCtyC
线性多步法误差分析（续）
?其中
? （ 25）
?如 均为 0，则称为 p阶公式
? （ 26）
kC ???? ????? ?2100
? ?)()2()1( 10!0121!1111 kkkC ?????? ????????? ??
? ?)2()2()1( 21!112221!2122 kk kkC ?????? ????????? ??
? ?)2()2()1( 2221!213231!3133 kk kkC ?????? ????????? ??
?
? ?)2()2()1( 1211)!1( 121!1 kmmmkmmmmm kkC ?????? ??? ????????? ??
C C C C p0 1 2,,,,?
?? ????? ????
? ? ??? ?
)()()( )1(11
0 0 kn
p
p
pk
i
k
i knikni
tyChihtyhihty ??
线性多步法误差分析（续）
? 下面，如果我们能证明上一节推导出的公式若能满
足求 均为 0的条件，则就得出了这些公
式的截断误差满足（ 26）式。
? 以三阶公式为例，将（ 25）式与
相关表达式表示成右端为零，可得,
C C C C p0 1 2,,,,?
32103210,,,,,,,????????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
???
???
0
0
0
0
0
0
0
0
37271703210
36261603210
35251503210
34241403210
33231303210
32221203210
11113210
00001111
3
2
1
0
3
2
1
0
666
555
44455
33344
22233
22
?
?
?
?
?
?
?
?
?先讨论预报公式，由于 和,这意
味着要将上述矩阵的第 1列移到等式的右边，
并去掉第 5列。为使矩阵成为方阵，将其最后
两行也去掉。
10 ??? 0
0??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
???
???
0
0
0
0
0
0
0
0
37271703210
36261603210
35251503210
34241403210
33231303210
32221203210
11113210
00001111
3
2
1
0
3
2
1
0
666
555
44455
33344
22233
22
?
?
?
?
?
?
?
?
0
1
22
2
3 3 2 2 2
3
4 4 3 3 3
0
5 5 4 4 4
1
5 5 5
2
6 6 6
3
1 1 1 1 0 0 0 0
0 1 2 3 1 1 1 1
0 1 2 3 0 2 1 2 2 2 3
0 1 2 3 0 3 1 3 2 3 3
0 1 2 3 0 4 1 4 2 4 3
0 1 2 3 0 5 1 5 2 5 3
0 1 2 3 0 6 1 6 2 6 3
0 1 2 3 0 7 1 7 2 7 3
a
a
a
a
b
b
b
b
轾轾
犏犏
犏犏
犏犏
犏犏
创
犏犏
犏犏
创
犏犏
ê犏
创
ê犏
ê犏
创
ê犏
ê犏
创 ê犏
ê犏
创 ê犏
臌 臌
0
0
0
0
0
0
0
0
轾
犏
犏
犏
犏
犏
犏
犏
=
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
臌其结果与用于推导预报公式的矩阵方程 完全一样。
? 对校正公式，可采用类似的办法，只是,这样要将第 5
列移到等式的右边；若假定,则可去掉第 4列。同样
为得到方阵，去掉最后两行，结果就是 3阶校正公式。
10 ???
?3 0?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
???
???
0
0
0
0
0
0
0
0
37271703210
36261603210
35251503210
34241403210
33231303210
32221203210
11113210
00001111
3
2
1
0
3
2
1
0
666
555
44455
33344
22233
22
?
?
?
?
?
?
?
?
这就表明，上一节导出的预报与校正公式的截断误差系数可以用下式来计算
? ?)2()2()1( 21!11211)!1( 111 kpppkppppp kkC ?????? ????????? ????? ??（ 27）
2.4 稳定性分析
?仿真方法选择的基本要求,仿真计算不改变
原系统的绝对稳定性。
?原系统是稳定的。观察欧拉法仿真递推公式
?故有 （ 28）
? yn (n=0,1,2,?)为它的一个仿真解,
???? iydtdy ???,/ Re ? ?? ? 0
),(1 nnnn ythfyy ???
nnn yhyy ???? 1
稳定性分析（续）
?设 为其准确解，即
? （ 29）
?用 (29)式减去 (28)式，可得,
?即 （ 30）
?特征方程为 （ 31）
?显然，为了使扰动序列 ?n不随 n增加而增长，
必须要求,
?我们称它所对应的域就是该算法的稳定域,
h? 2?1/??，即 h小于等于系统时间常数的两
倍。
nny ??
)()()( 11 nnnnnn yhyy ???? ????? ??
nnn h ???? ??? 1
0)1(1 ???? nn h ???
z h? ? ?( )1 0?
11 ?? h?
确定数值积分法稳定域的一般方法
? 测试方程,
? 数值积分公式 (33)
? 其中 是一个关于 高阶多项式函数,
? 则只有当 时，算法才稳定。
0 -1 -2
Im ? h
Re ? h
图 2.8 稳定域
???? iydtdy ???,/
nn yhpy ??? )(1 ?
p( )h? h?
1)( ??hp
第 2章 经典的连续系统仿真建模方法学
剡昌锋 刘军
兰州理工大学机电工程学院
Tel:13369316141(剡） 13893121598（刘）
E_mail,changf_yan@163.com
Lj_gabe@sohu.com
2.1 离散化原理及要求
?问题：数字计算机在数值及时间上的 离散性 -
---被仿真系统数值及时间上的 连续性?
?连续系统的仿真，从本质上：对原连续系统
从时间、数值两个方面 对原系统进行离散化
并选择合适的数值计算方法来 近似 积分运算
? 离散模型 ≈原连续模型？
相似原理
?设系统模型为：,其中 u(t)为输
入变量,y(t)为系统变量；令仿真时间间隔
为 h，离散化后的输入变量为,系统变
量为,其中 表示 t=nh。
?如果,且
?即,
?（对所有 n=0,1,2,… ）
?则可认为两模型等价。
(,,)y f y u t=&
)(? ntu
)(? nty nt
)()(? nn tutu ? )()(? nn tyty ?
?( ) ( ) ( ) 0u n n ne t u t u t= - ?( ) ( ) ( ) 0y n n ne t y t y t= -
u(t)
h
y(t)
-
+
图 2.1 相似原理
原连续模型 (,,)y f y u t=&
仿真模型 ? ? ?(,,)ny f y u t=&
)(? nty
0)( ?ny te
)(? ntu
对仿真建模方法三个基本要求
? （ 1） 稳定性,若原连续系统是稳定的，则离散化
后得到的仿真模型也应是稳定的。
? （ 2） 准确性,有不同的准确性评价准则，最基本
的准则是,
? 绝对误差准则,
? 相对误差准则,
? 其中 ? 规定精度的误差量。
???? )()(?)( nnny tytyte
???? )(? )()(?)(
n
nn
ny ty
tytyte
对仿真建模方法三个基本要求 (续 )
? 3） 快速性,若第 n步计算对应的系统时间间
隔为 计算机由计算需要的时间为 Tn，
若 Tn=hn 称为 实时仿真, Tn?hn称为 超实
时仿真 Tn?hn 称为 亚实时仿真,对应于 离
线仿真
,1 nnn tth ?? ?
数值积分算法
? 对,已知系统变量 y的初始条
件,要求 y随时间变化的过程 ―― 初值问
题
? 计算过程：由初始点 的
? 欧拉法
? 对任意时刻 tn+1
? 截断误差正比于
(,,)y f y u t=&
00()y t y=
00()y t y= )( 00 ytf,
0
0( ) (,)
t
ty t y f t y d t=+ ò
1 1 0 0 0( ) ( )y y t y t f t y= @ + D,
1 1 1( ) ( ) ( )n n n n n n ny y t y t t f t y+ + += @ + -,
h2
f ( t,y )
f ( t
0
,y
o
)
tt
t
0 t 1
数值积分算法 (续）
?梯形法,
?是隐函数形式。预报－ — 欧拉法估计初值，
校正－ — 用梯形法校正,
?校正公式
?预报公式
?反复迭代，直到满足
?经典的数值积分法可分为两类,
?单步法与多步法
)],(),([21)( 1111 ???? ???? nnnnnnn ytfytfhytyy
)],(),([21 )( 11)1( 1 innnnnin ytfytfhyy ???? ???
),()( 1 nnnin ytfhyy ????
??? ??? inin yy 111
2.2 龙格库塔法
? 2.2.1龙格 -库塔法基本原理
?对
?若令,
?则有
? 的数值求解：称作“右端函数”计算问题。
在 附近展开台劳级数，只保留 项，则有,
（ 1）
? ???? 1 ),()()( 1 nnttnn dtytftyty
)( nn tyy ? ? ?? 1 ),(Q n
n
t
tn dtytf
nnnn yyty Q)( 11 ??? ??
nQ
t0 h2
2
0001 0)(2
1),( h
t
f
dt
dy
y
fhytfyy
t?
?
?
? ????
龙格 -库塔法基本原理（续）
?假设这个解可以写成如下形式,
?其中
?对 式右端的函数展成台劳级数，保留 h项，
可得,
?代入，则有,
1 0 1 1 2 2()y y a k a k h= + +
),( 001 ytfk ? k f t b h y b k h2 0 1 0 2 1? ? ?( )，
k2
02 0 0 1 2 1(,) ( ) t
ffk f t y b b k hty抖@ + +抖
01 0 1 0 0 2 0 0 1 2 1(,) [ (,) ( ) ]t
ffy y a h f t y a h f t y b b k hty抖= + + + +抖
龙格 -库塔法基本原理（续）
?将 (2)式与 (1)式进行比较，可得,
?四个未知数 但只有三个方
程，因此有无穷多个解。
?若限定,则
?计算公式,
?其中
a a a b a b1 2 2 1 2 21 1 2 1 2+ = = =,/,/
，，，,2121 bbaa
a a1 2? a a b b1 2 1 212 1? ? ? ?，
)(2 2101 kkhyy ???
)(),( 1002001 hkyhtfkytfk ????,，
龙格 -库塔法基本原理（续）
?若写成一般递推形式，即为,
?其中
?截断误差正比于 h3，称为二阶龙格 -库塔法
（简称 RK-2）。
?截断误差正比于 h5的四阶龙格 --库塔法（简称
RK-4）公式,
?其中,
)(2)( 2111 kkhyyty nnn ???? ??
),(),( 121 hkyhtfkytfk nnnn ????,
)22(6)( 432111 kkkkhyyty nnn ?????? ??
),(1 nn ytfk ? )
22( 12 k
hyhtfk
nn ???,
)22( 23 khyhtfk nn ???, )( 34 hkyhtfk nn ???,
2.2.2龙格 --库塔法的特点
? 1.形式多样性
? 例,非唯一解，可以得到许多种
龙格 --库塔公式,(中点公式 )
? 其中
? 各种龙格 ---库塔法可以写成如下一般形式,
? 其中,
2121 bbaa,，，
hkyy nn 21 ???
),(1 nn ytfk ? )
22( 12 k
hyhtfk
nn ???,
?
?
? ??
s
i
iinn kChyy
1
1
)(
1
1
??
?
???
i
j
jijnini kbhyhatfk, si,，,?21?
龙格 --库塔法的特点 (续）
? 式中各系数满足以下关系
? s称为级数，表示每步计算右端函数 f的最少次数。
可以证明,1阶公式至少要计算一次,2阶公
式 ； …,； 4阶公式 ；依此类推。有时
为了某种特殊需要，可以选择 的计算公式。
a
a b i s
C
i ij
j
i
i
i
s
1
1
1
1
0
2 3
1
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
，,,?
smin ? 2 smin ? 4
minss ?
龙格 --库塔法的特点 (续）
? 2.单步法
? 在计算 时只用到,而不直接用
等项。优点：存储量减小，可以自启动
? 3.可变步长
? 步长 h在整个计算中并不要求固定，可以根
据精度要求改变
? 但是在一步中，为计算若干个系数,则
必须用同一个步长 h。
1?ny ny 21 ?? nn yy,
ik
龙格 --库塔法的特点 (续）
? 4.速度与精度
?四阶方法的 h可以比二阶方法的 h大 10倍，每
步计算量仅比二阶方法大一倍
? 高于四阶的方法由于每步计算量将增加较
多，而精度提高不快。
2.2.3实时龙格－库塔法
?实时仿真,要求仿真模型的运行速度往往与
实际系统运行的速度保持一致。
?一般的数值积分法 难以满足 实时仿真的要求，
这不仅仅是因为由这些方法所得到的模型的
执行速度较慢，而且这些方法的机理不符合
实时仿真的特点。
?考虑系统 ))(( tuyf
dt
dy,?
实时龙格－库塔法 (续）
? RK-2公式如下,
? 一个计算步内分两子步,
? tn时刻：利用当前的 un,yn计算 k1----计算一次右端
函数 f需
? tn+h/2时刻：应计算 k2，尽管此时 yn +1/2已经得到，
但 un +1则无法得到。 (若对 un +1也进行预报 ―― 加大
仿真误差 )。仿真执行延迟 h/2―― 输出要迟后半个
计算步距
),,(
),,(
)(
2
1112
1
211
hkyutfk
yutfk
kk
h
yy
nnn
nnn
nn
??
?
???
??
?
h / 2
实时龙格－库塔法 (续）
RK-2的计算流程
并输出计算 1?ny1?nu采入
1k计算
nt
2
ht
n ?
1??? nn tht
21
ht
n ??
21 ?? ?? nn tht
2k计算 1k计算下一个
实时龙格－库塔法 (续）
?实时 2阶龙格－库塔法,
? tn时刻：应计算 k1，利用当前的 un,yn，需
要 ； tn+h/2时刻，应计算 k2，此时 yn +1/2已
经得到,un +1/2也可得到,k2的计算就不会引
入新的误差。计算一次右端函数 需要,
可实时输出 yn +1。
)
2
,,(
),,(
12/12/12
1
21
k
h
yutfk
yutfk
hkyy
nnn
nnn
nn
??
?
??
??
?
h/2
f h/2
实时龙格－库塔法 (续）
实时 RK-2公式计算流程
)2/( htu n ?采入 并输出计算 1?ny
nt 2htn ? 1??? nn tht
21
ht
n ?? 21 ?? ?? nn tht
1k计算 2k计算 1k计算下一个
2.3 线性多步法
? 2.3.1线性多步法基本原理
?基本原理,利用一个多项式去匹配变量若干
已知值和它们的导数值。
?设,时刻的 和
已知 ;
?预报：由 和 来计算
?校正：若 也已知，由它们来计算
t t tn n n k,,,? ? ?1 1? y y yn n n k,,,? ? ?1 1? ?,?,,?y y yn n n k? ? ?1 1?
y y yn n n k,,,? ? ?1 1? ?,?,,?y y yn n n k? ? ?1 1? y yn k n k? ?,?
kny? ?yn k?
线性多步法 (续 )
? 采用的多项式具有以下形式 (m阶 )
? 其中,是待定系数,
? 在 时刻，,可得到,
? （ 1）
y t d
t t
h
d
y t id
t t
h
id
m i
n k
i
i
i
i
m
i
m
m h i
n k
i
m i
h i
i
i
m
( )
? ( )
?
??
?
?
?
?
? ?
? ?
??
?
?
?
?
? ? ?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
??
? ?
?
?
00
1
1
1
1 1
1
di h tt kn ?? ??
tn k? 0?? y dn k? ? 0
?y dn k h? ? ? 1 1
线性多步法 (续 )
?由 和 确定,
需要 m+1个独立方程。该 m+1个方程可由以
下等式导出,
? （ 2）
y y yn n n k,,,? ? ?1 1? ?,?,,?y y yn n n k? ? ?1 1? ),2,1,0( mid
i ??
y t y d
t t
h
d j
y t y id
t t
h
id j
m n k j i
n k n k j
i
i
i
i
m
i
m
m n k j h i
n k n k j
i
m
i
h i
i
i
m
( )
? ( ) ?
? ?
??
?
?
?
?
? ?
? ? ?
??
?
?
?
?
? ? ?
?
?
?
?
?
?
?
? ?
? ? ?
??
? ?
? ? ?
?
?
?
?
??
? ?
00
1
1
1
1 1
1
线性多步法 (续 )
? 1、预报公式
? 令 m=2k-1，从（ 2）式得到如下方程组,
d d d d yp p p mp n k0 1 2 1? ? ? ? ? ? ??
d d d d yp p p m mp n k0 1 2 22 4 2? ? ? ? ? ? ??
d d d d yp p p m mp n k0 1 2 2 33 3 3? ? ? ? ? ? ??
?
d d md h yp p mp n k1 2 12? ? ? ? ? ? ?? ?
d d m d h yp p m m p n k1 2 1 22 2 2? ? ? ? ? ? ?? ? ?? ?
d d m d h yp p m m p n k1 2 1 32 3 3? ? ? ? ? ? ?? ? ?? ?
?
? 将其写成矩阵形式,（ 4）
? 其中上标 p表示预报。
? 其解为, （ 5）
V d Zp p p?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
??
??
??
??
??
??
?
?
?
?
?
?
n
kn
kn
kn
n
kn
kn
kn
m
k
k
k
k
m
m
m
m
m
m
yh
yh
yh
yh
y
y
y
y
d
d
d
d
d
d
d
d
kmkk
m
m
m
kkkk
?
?
?
?
?
?
?
?
?
??
?
?
?
?
??
?
?
?
3
2
1
3
2
1
2
1
1
2
1
0
12
12
12
32
32
32
3210
3233210
2232210
3210
1
33331
22221
11111
（ 3）
? ?d V Zp p 1 p? ?
其解为,? ?
d V Zp p 1 p? ?
由于 为常数阵，其逆存在,Z向量中的各元素为
已知值，因而 d向量的各元素值可计算得到，从而
由,得到下一时刻的预报值 。
Vp
y d y dn k n k h? ?? ? ?0 1 1,?
缺点：只有 是所需要的，其它元素的计算成
为多余,得不到 与 和
显式表达式
d d0 1,
y yn k n k? ?,? y y yn n n k,,,? ? ?1 1?
?,?,,?y y yn n n k? ? ?1 1?
?定义,（ m+1） ?1的列向量
? (6)
?定义辅助变量 (7)
?此式可改写为 (8)
?向量 的元素可划分为两个组
? (9)
T)0,,0,1( ??0e
? ?y n k? ?? ?e d e V Z0T p 0T p 1 p
? ? ? ?? p T ? ?e V0T p 1
? ?V ep T 0? p ?
?p
? ?? p p p k p p p k p Ta a a b b b? 1 2 1 2,,,,,,? ? ?
?例,k=3，则 (8)式为,
?可计算得到,
1 1 1 0 0 0
1 2 3 1 1 1
1 2 3 2 2 2 2 3
1 2 3 3 3 2 3 3
1 2 3 4 4 2 4 3
1 2 3 5 5 2 5 3
1
0
0
0
0
0
2 2
3 3 2 2
4 4 3 3
5 5 4 4
1
2
3
1
2
3
? ?
? ?
? ?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
a
a
a
b
b
b
p
p
p
p
p
p
? ? ? ?? p p p p p p p Ta a a b b b? ? ? ? ? ?1 2 3 1 2 3 18 9 10 9 18 3,,,,,,,,,,? ?
? 只依赖于 k，即先前 和 的个
数，而与它们的数值无关。这样，可以预
先求解（ 8）式得到。
?从而得到 的显式表达式,
?p yn k j? ? ?yn k j? ?
? ?y n k p T? ? ?e d Z0T p p?
yn k?
y a y h b yn k jp n k j jp n k j
j
k
j
k
? ? ? ? ?
??
? ? ?? ?
11
?例：试推导用 预报 公式
?条件：已知
?由
y y yn k n k n k? ? ? ? ? ?1 1 2,?,? yn k?
? ?V ep T 0? p ?
V p ?
?
?
?
?
?
?
?
?
?
?
1 1 1
0 1 2
0 1 4
y y yn k n k n k? ? ? ? ? ?1 1 2,?,?? ? ? ?? p T p p pa b b?
1 1 2,,
? ?V p T ? p
p
p
p
a
b
b
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1 0 0
1 1 1
1 2 4
1
0
0
1
1
2
a b bp p p1 1 21 3 2 1 2? ? ? ?,/,/
? ?y y h y yn k n k n k n k? ? ? ? ? ? ?? ? ?1 12 1 23 ? ?
2、校正公式
?预报公式 ―― 显式公式，未包括 。
?校正,对该预报值应进行校正，即先预报得
到,然后再用此值推出
?由 和 以及 来预
报,
?可令 m=2k-1，从（ 2）式得到如下方程组,
?yn k?
~yn k? ?yn k?
y y yn n n k,,,? ? ?1 1?,~ kny ? ?,?,,?y y yn n n k? ? ?1 1?
?yn k?
kn
c yd
?? ~0
?将其写成矩阵形式,
d d d d yc c c mc n k0 1 2 1? ? ? ? ? ? ??
d d d d yc c c m mc n k0 1 2 22 4 2? ? ? ? ? ? ??
d d d d yc c c m mc n k0 1 2 2 33 3 3? ? ? ? ? ? ??
?
d d md h yp c mc n k1 2 12? ? ? ? ? ? ?? ?
d d m d h yc c m mc n k1 2 1 22 2 2? ? ? ? ? ? ?? ? ?? ?
d d m d h yc p m mc n k1 2 1 32 3 3? ? ? ? ? ? ?? ? ?? ?
?
V d Zc c c?
校正公式 (续 )
?其中上标 c表示校正，可得
? (14)
? ?d V Zc c 1 c? ?
1 0 0 0 0
1 1 1 1 1
1 2 2 2 2
1 3 3 3 3
1
0 1 2 3
0 1 2 2 3 2 2
0 1 2 3 3 3 3
2 3
2 3
2 3
2 1
2 1
0
1
2
3
1
2
?
?
?
?
?
?
?
?
?
?
?
?
m
m
m
m
m
k
k
k
m
k k k k
m
m
m
d
d
d
d
d
d
d
d
? ? ?
? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ?
? ?
? ?
? ?
? ?
? ?
y
y
y
y
y
h y
h y
h y
n k
n k
n k
n k
n
n k
n k
n k
1
2
3
1
2
3
?
?
?
?
?
?定义,为（ m+1） ?1的列向量，
上标 T表示转置。将 左乘（ 13）式可得,
? （ 15）
?定义 （ 16）
?可改写为 （ 17）
? （ 18）
e 1T ? (,,,,)0 1 0 0?
e1T
? ?? ? ?? ?h y n k? e d e V Z0T c 1T c 1 c
? ? ? ?? c T ? ?e V1T c 1
? ?V ec T 1? c ?
? ?? c c c kc c c kc Ta a a b b b? 0 1 1 2,,,,,,? ? ?
?例,k=3
?同样,只依赖于 k，即先前 和 的个
数，而与它们的数值无关。这样
? （ 19）
?从而 （ 20）
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
??
0
0
0
0
1
0
35255210
34244210
33233210
32222210
111210
000111
3
2
1
3
2
1
445
334
223
2
p
p
p
p
p
p
b
b
b
a
a
a
?c yn k j? ? ?yn k j? ?
? ?? ? ??h y n k c T? e d Z1T c c?
? ( ? )y h a y h b yn k jc n k j jc n k j
j
k
j
k
? ? ? ? ?
??
? ? ? ??1
10
?例：已知,预估,然后用
校正 。
预估
预报公式为
校正
校正公式为
y y yn k n k x k? ? ? ? ? ?1 2 3,,kny ?~
21,,~ ????? knknkn yyy ?yn k?
1 1 1
1 2 4
1 3 9
0
1
2
1
2
3
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ?
? ?
? ?
d
d
d
y
y
y
p
p
p
n k
n k
n k
? ? ? ? ? ?? p T p p pa a a? ? ?1 2 3 3 3 1,,,,
321 33~ ??????? ??? knknknkn yyyy
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
2
1
2
1
0
~
421
111
001
kn
kn
kn
c
c
c
y
y
y
d
d
d
? ? ? ? ? ?2123210,2,,,???? cccTc aaa?
)2~( 2211231 ?????? ????? knknknhkn yyyy?
2.3.2 线性多步法误差分析
? 为了便于分析，对预报公式和校正公式，定义统一的表达式,
? （ 21）
? -------显式预报 显式预报
? 时称为后向差分公式（ BDF）
? 同时均不等于 0时为隐式校正公式
? k称为公式的阶次 。假设变量各时间的 精确值 已经得到，将
其代入（ 21）式，可得,
? （ 22）
? ?
? ?
???? ??
k
i
k
i
ikniikni yhy
0 0
0???
00 ?? ?yn k? 00 ?? yn k?
021 ???? k??? ?
00,??
? ?
? ? ??
???k
i
k
i knikni
ihtyhihty
0 0
)()( ???
线性多步法误差分析（续）
?在 附近，将每个函数展开成泰勒级数,
? （ 23）
?对所有 i(i=0,1,2,…,k), 将（ 23）式代入（ 22）
式，合并同类项，可得
?
? （ 24）
tn k?
???? ????? ???? )()()()( !2 )(!1 2 knihknihknkn tytytyihty
???????? ????? ???? )()()()( !2 )(!1 2 knihknihknkn tytytyihty
????? ????? ???? )()()()( )(2210 knmmmknknkn tyChtyChtyhCtyC
线性多步法误差分析（续）
?其中
? （ 25）
?如 均为 0，则称为 p阶公式
? （ 26）
kC ???? ????? ?2100
? ?)()2()1( 10!0121!1111 kkkC ?????? ????????? ??
? ?)2()2()1( 21!112221!2122 kk kkC ?????? ????????? ??
? ?)2()2()1( 2221!213231!3133 kk kkC ?????? ????????? ??
?
? ?)2()2()1( 1211)!1( 121!1 kmmmkmmmmm kkC ?????? ??? ????????? ??
C C C C p0 1 2,,,,?
?? ????? ????
? ? ??? ?
)()()( )1(11
0 0 kn
p
p
pk
i
k
i knikni
tyChihtyhihty ??
线性多步法误差分析（续）
? 下面，如果我们能证明上一节推导出的公式若能满
足求 均为 0的条件，则就得出了这些公
式的截断误差满足（ 26）式。
? 以三阶公式为例，将（ 25）式与
相关表达式表示成右端为零，可得,
C C C C p0 1 2,,,,?
32103210,,,,,,,????????
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
???
???
0
0
0
0
0
0
0
0
37271703210
36261603210
35251503210
34241403210
33231303210
32221203210
11113210
00001111
3
2
1
0
3
2
1
0
666
555
44455
33344
22233
22
?
?
?
?
?
?
?
?
?先讨论预报公式，由于 和,这意
味着要将上述矩阵的第 1列移到等式的右边，
并去掉第 5列。为使矩阵成为方阵，将其最后
两行也去掉。
10 ??? 0
0??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
???
???
0
0
0
0
0
0
0
0
37271703210
36261603210
35251503210
34241403210
33231303210
32221203210
11113210
00001111
3
2
1
0
3
2
1
0
666
555
44455
33344
22233
22
?
?
?
?
?
?
?
?
0
1
22
2
3 3 2 2 2
3
4 4 3 3 3
0
5 5 4 4 4
1
5 5 5
2
6 6 6
3
1 1 1 1 0 0 0 0
0 1 2 3 1 1 1 1
0 1 2 3 0 2 1 2 2 2 3
0 1 2 3 0 3 1 3 2 3 3
0 1 2 3 0 4 1 4 2 4 3
0 1 2 3 0 5 1 5 2 5 3
0 1 2 3 0 6 1 6 2 6 3
0 1 2 3 0 7 1 7 2 7 3
a
a
a
a
b
b
b
b
轾轾
犏犏
犏犏
犏犏
犏犏
创
犏犏
犏犏
创
犏犏
ê犏
创
ê犏
ê犏
创
ê犏
ê犏
创 ê犏
ê犏
创 ê犏
臌 臌
0
0
0
0
0
0
0
0
轾
犏
犏
犏
犏
犏
犏
犏
=
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
ú 犏
臌其结果与用于推导预报公式的矩阵方程 完全一样。
? 对校正公式，可采用类似的办法，只是,这样要将第 5
列移到等式的右边；若假定,则可去掉第 4列。同样
为得到方阵，去掉最后两行，结果就是 3阶校正公式。
10 ???
?3 0?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
???
???
???
???
0
0
0
0
0
0
0
0
37271703210
36261603210
35251503210
34241403210
33231303210
32221203210
11113210
00001111
3
2
1
0
3
2
1
0
666
555
44455
33344
22233
22
?
?
?
?
?
?
?
?
这就表明，上一节导出的预报与校正公式的截断误差系数可以用下式来计算
? ?)2()2()1( 21!11211)!1( 111 kpppkppppp kkC ?????? ????????? ????? ??（ 27）
2.4 稳定性分析
?仿真方法选择的基本要求,仿真计算不改变
原系统的绝对稳定性。
?原系统是稳定的。观察欧拉法仿真递推公式
?故有 （ 28）
? yn (n=0,1,2,?)为它的一个仿真解,
???? iydtdy ???,/ Re ? ?? ? 0
),(1 nnnn ythfyy ???
nnn yhyy ???? 1
稳定性分析（续）
?设 为其准确解，即
? （ 29）
?用 (29)式减去 (28)式，可得,
?即 （ 30）
?特征方程为 （ 31）
?显然，为了使扰动序列 ?n不随 n增加而增长，
必须要求,
?我们称它所对应的域就是该算法的稳定域,
h? 2?1/??，即 h小于等于系统时间常数的两
倍。
nny ??
)()()( 11 nnnnnn yhyy ???? ????? ??
nnn h ???? ??? 1
0)1(1 ???? nn h ???
z h? ? ?( )1 0?
11 ?? h?
确定数值积分法稳定域的一般方法
? 测试方程,
? 数值积分公式 (33)
? 其中 是一个关于 高阶多项式函数,
? 则只有当 时，算法才稳定。
0 -1 -2
Im ? h
Re ? h
图 2.8 稳定域
???? iydtdy ???,/
nn yhpy ??? )(1 ?
p( )h? h?
1)( ??hp
