系统仿真技术
第 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,… )
?则可认为两模型等价。
),,( tuyfy ??
)(? ntu
)(? nty nt
)()(? nn tutu ? )()(? nn tyty ?
0)()(?)( ??? nnnu tutute 0)()(?)( ??? nnny tytyte
u(t)
h
y(t)
-
+
图 2.1 相似原理
原连续模型 ),,( tuyfy ??
仿真模型 ),?,?(? ntuyfy ??
)(? 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
? 截断误差正比于
),,( tuyfy ??
00 )( yty ?
00 )( yty ? )( 00 ytf,
??? tt dtytfyty 0 ),()( 0
y y t y t f t y1 1 0 0 0? ? ? ?( ) ( )?,
)()()( 111 nnnnnnn ytfttytyy,????? ???
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项,
可得,
?代入,则有,
y y a k a k h1 0 1 1 2 2? ? ?( )
),( 001 ytfk ? k f t b h y b k h2 0 1 0 2 1? ? ?( ),
k2
hyfkbtfbytfk t 0)(),( 121002 ???????
])(),([),( 012100200101 hyfkbtfbytfhaythfayy t?????????
龙格 -库塔法基本原理(续)
?将 (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
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
?
?
?
?
?
?
?
?
其结果与用于推导预报公式的矩阵方程 完全一样。
? 对校正公式,可采用类似的办法,只是,这样要将第 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