系统仿真技术
第 6章 病态系统仿真
剡昌锋 刘军
兰州理工大学机电工程学院
6.1病态系统的定义
?系统中各环节的时间常数差异巨大。为保证
仿真计算的稳定性，由于仿真步长必须限制
在最小时间常数的数量级而选得很小，然而
仿真结束的时间则决定于系统中的最大时间
常数，若按满足稳定性要求所选择的步长进
行仿真，则不仅整个仿真所花费的时间非常
长，甚至由于计算的舍入误差而导致整个仿
真的失败。这就是所谓,病态系统仿真,问
题。
病态系统的定义（续）
? 病态系统定义, （ 1）
? 令
? （ 2）
? J称为系统的雅可比矩阵。
? 若 J的特征值全部具有负实部，且有：,
则该系统称为病态系统，在某些文献中也叫做刚性
系统 (Stiff)，而,
? 称为病态比，一般在 50以上。
),( tyFy ??
1
1
1
1
J
y
F
n
nn
n
y
f
y
f
y
f
y
f
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
?
ii ?? m a xm in ??
i
iS
?
?
min
max?
6.2 线性病态系统仿真
?对线性定常系统，我们可用如下状态方程进
行一般性描述,
? (3)
?1,增广矩阵法 将 u(t)作为系统的增广状态
?
?
?
?
?
?
?
?
??
)0()(
)()(
)()(
0 xtx
tCxty
tButAxx?
?
?
?
??
?
?
?
?
)0(~)(~
)(~C
~
=( t )
)(~
~~
0 xtx
txy
txAx?
)(~)(~)1(~ kxtkx ???
线性病态系统仿真 (续）
? 其中,
? 由于病态系统特征值相差倍数很大，必须
用加速收敛的方法计算该状态转移矩阵的
值。
? 2.蛙跳算法 基本思想,
? （ 1）考虑作用函数为阶跃函数，则增广状
态十分简单。
? 选择 q，使得,
e hA~~ =(T )?
hqhq nn 2,2 ??
线性病态系统仿真 (续）
?(3) 仿真计算时采用如下“蛙跳”方式,
??
??
)2()3(
)()2(
)0()0()()2(
)0()0()4(
)0()0()2(
)0()(
~
2
~
2
~
2
~
2
~
2
~
4
~
2
~
2
~
2
~~
~
11
qhxeqhx
qhxeqhx
xexeeqhxhx
xexeehx
xexeehx
xehx
hA
hA
hAhAhAn
hAhAhA
hAhAhA
hA
n
n
nnn
?
?
???
?
?
????
???
???
?
??
线性病态系统仿真 (续）
? 即在 qh以前采用加倍跳跃式计算，而在 qh以后每隔
qh计算一次。
? 优点是,h可以取得很小（可按最小时间常数考
虑），从而保证初始阶段的精度而计算量却不大，
而到 qh以后，小时间常数的作用完成，则加大步长
计算，从而加快仿真计算速度。
? 从 x(h)到 x(qh)都是以 x(0)为基础进行计算，所以误
差传播比较均匀（仅仅是状态转移矩阵的误差）。
6.3 非线性病态系统仿真
?一般非线性系统的仿真大多采用数值积分法。
而数值积分法一般又只具有有限的稳定域，
典型的如 龙格－库塔法,仿真步长限定在系
统最小时间常数的数量级，才能保证计算的
稳定性，而系统的过渡过程时间却决定于最
大时间常数，因而对病态系统来说计算量极
大，加上存在误差传播，仿真的精度甚至稳
定性也会受到影响。
6.3.1 吉尔 (Gear)法
? 6.3.1.1.Stiff稳定域
? Gear研究后发现，并不要求一定采用恒稳方法，
而只要具有所谓 Stiff稳定域就可以了。
? Stiff稳定域定义：对实际的物理系统，时间常数 μ
一般小于零。选择仿真步长 h若满足,
? 与
? 可保证仿真的稳定性，称该算法具有 stiff稳定域。
(I) Im 0Re
??
?
?
?
??
?
h
h
( I I ) Re ?? ??h
Stiff稳定域（续）
? 实际上具有 Stiff稳定域的方法
与恒稳方法只在近虚轴处有一
点差别，即如果系统中的极点
全部为实极点，那么无论选择
多大的步长，计算是恒稳的。
如果系统中有复极点（实部仍
为负数），只要步长的选择满
足上述条件，也能保证算法稳
定。
0
θ
-θ
Reμ h
I
II
-δ
图 6.1 Stiff稳定域
Stiff稳定域（续）
?Stiff域中 θ与 δ的确定,
?按病态系统的大特征值来选择步长,
?该特征值所对应的模态大约要经过 4倍左右时
间常数的时间才能有效地衰减掉，即
?,也就是
?这样，此时即使加大步长 h，也能保持计算的
稳定性。基于这一考虑，可设 δ＝ - 4。
m a xm a xm a x ??? j??
m a x/4 ??h 4max ??h
Stiff稳定域（续）
? 另一方面，考虑到系统特征值为复数，它所对应的
瞬态响应呈振荡型。一个振荡周期内至少计算 N个点。
最小振荡周期为,
? 其中 h为计算步长，若选择 N≥8，则有：,因
此可选 θ＝ π/4。
? 综上所述，如果选择某一种方法，其稳定域 θ≥π/4，
且 |δ|≤4，则从使用的角度来看，图 6.1所示的稳定域
与恒稳域没有差别，从而完全可以用于病态系统的
仿真。
NhT ?? m a xm i n /2 ??
4/ma x ?? ?h
6.3.1.2吉尔（ Gear）法的基本原理
? 设系统,
?满足 Stiff稳定域的多步法,Gear提出的用于
病态系统仿真的计算公式是,
? (1)
?
00 )( ),(/ yttytyfdtdyy ?????
?
?
???? ??
k
j
knjknjkn yhyay
1
0 ??
用于病态系统仿真的 Gear公式的系数表
α1 α2 α3 α4 α5 α6 β0
G-1 1 1
G-2 4/3 -1/3 2/3
G-3 18/11 -9/11 2/11 6/11
G-4 48/25 -36/25 16/25 -3/25 12/25
G-5 300/137 -300/137 200/137 -75/137 12/137 60/137
G-6 360/147 -450/147 400/147 -225/147 72/147 -10/147 60/147
吉尔（ Gear）法的基本原理 (续）
? 稳定域如图 6.2所示。
? 从图上可以看出，该方
法在 5阶以下（包括 5阶）
的稳定域满足 Stiff稳定
域的条件
（ θ>π/4,|δ|<4）,
? 大于 5阶时 |δ|>4，而且
还可能穿过负实轴。 -10 -8 -6 -4 -2 0 2 4 6 8 10
0
2
4
6
8
10
图 6.2 用于病态系统仿真的 G e a r 法
C
2
C
3
C
4
C
5
C
6
吉尔（ Gear）法的基本原理 (续）
? 在用 Gear法仿真非线性病态系统时，有以下三个基本问题需要解决,
? 1） 启动问题
? 上述 Gear法本质上是隐式多步法。
? 对于初值问题，困难：隐式方法一般用显式方法启动，即先进行预报，
然后通过迭代进行校正。如果迭代方法的收敛性不好，可能引起计算发散或计算量加大。
? 即使选择的迭代方法收敛性满足要求，显式多步法预报，仍然难以启
动，必须采用单步法启动，由于单步法不具有 Stiff稳定域，因而很难
保证计算的稳定性。
? 2） 变步长策略
? 非线性病态系统仿真往往采用变步长策略，如何适时地将步长调整到
合适长度，以同时满足仿真精度和速度的要求。
? 3） 加速迭代 为了提高计算效能，加速迭代也是非线性病态系统仿
真中重要问题。
6.3.1.3 单步多值法
? 以三阶为例，采用显式多步法进行预报，然后用隐
式法校正。
? 其显式预报的公式是,（ 2）
?三阶隐式 Gear公式校正,
? （ 3）
kkkkk hfyyyy 33 221123)0( 1 ????? ???
?
)1(
111
6
211
2
111
9
11
18)2(
1
)0(
111
6
211
2
111
9
11
18)1(
1
????
????
????
????
kkkkk
kkkkk
hfyyyy
hfyyyy
单步多值法（续）
?其中等式右边第 4项为导函数项，它是通过将
第 i次迭代所得到的 y的预报值代入导函数后
计算得到的。为便于程序实现，由（ 2）式，
并令,
24
5
14
23
2
11)1(
1
)1(
1
)0(
111
6
22
1
12
3
)1(
1
)0(
111
6)0(
1
)0(
111
6
211
2
111
9
11
18)1(
1
7
)(33-=
)(
??
?
?
?
????
?
???
????
????
?????
???
????
kkkkk
kkkkkk
kkk
kkkkk
yyyhff
hfhfhfyyy
hfhfy
hfyyyy
单步多值法（续）
?迭代的校正公式可表示成,
?
?
)(
)(
)(
)1(
1
)(
111
6)(
1
)1(
1
)0(
1
)1(
111
6)1(
1
)2(
1
)1(
1
)0(
111
6)0(
1
)1(
1
?
???
?
?
????
?
????
???
???
???
i
k
i
k
i
k
i
k
kkkk
kkkk
hfhfyy
hfhfyy
hfhfyy
单步多值法（续）
?对一般情形，令,
Tpkkkkkk yyyhfyY ),,,,,( 121 ????? ?
Tpkkkkkk yyyhfyY ),,,,,( 21111 ?????? ? ?
),(),(= 1)1( 11)( 1
)1(
1
)(
1
)(
1
?
?
???
?
???
?
??
k
i
kk
i
k
i
k
i
k
i
k
tyhftyhf
hfhfG
单步多值法（续）
? 其中 p表示 Gear法的阶次。对三阶预报－迭代校正
算法,
? 用矩阵的方法加以表示,
? （ 4）
? （ 5）
? 其中 B,C分别是相应的系数矩阵。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
1
4
5
2
11
4
23
2
1
2
3
1
)1(
1
)0(
1
0100
0001
7
33
k
k
k
k
k
k
k
k
y
y
hf
y
y
y
hf
y
)(
1
11
6
1
)1(
1
)(
1
1
)(
1
)1(
1
0
0
1 i
k
k
k
i
k
i
k
k
k
i
k
i
k
G
y
y
hf
y
y
y
hf
y
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
),.,.,k -,( i =YBY k)(k 121 0 1 ??
)121( 1111,.,.,k -,i=GCYY )i(k)i(k)i(k ???? ??
单步多值法（续）
?然而，对初值问题，上式是不能自启动的。
?解决方法 ------单步多值法：用高阶导数值来
取代前几步的 y及 f的值。
? 先定义一个向量，称之为 Nordsieck向量,
? （ 6）
? （ 7）
?需要确定 Z向量与 Y向量之间的关系。
),,,,( )(!)2(!2 2 Tpkphkhkkk yyyhy pZ ???
),,,,( )( 1!)2(!2111 2 Tpkphkhkkk yyyhy pZ ???? ? ??
单步多值法（续）
?采用多项式逼近，以三阶为例,
? （ 8）
?在 处，有
? （ 9）
6)(
)(62)(
)(3)(2)(
)()()()(
3
32
2
321
3
3
2
210
?
?
?
?
?
?
?
?
???
?????
???????
aty
ttaaty
ttattaaty
ttattattaaty
k
kk
kkk
???
??
?
单步多值法（续）
?同样，也可以得到,
?由（ 8）式及（ 9）式消去,整理后
可得到,
!3,!2,,3210 ayayayay kkkk ???? ???
3210,,,aaaa
?
?
?
?
?
?
?
?????
?????
?
?
??
??
24
1
12
13
4
3
!3
24
1
12
3
4
7
!2
3
2
2
kkkkk
h
kkkkk
h
kk
kk
yyhfyy
yyhfyy
hfyh
yy
???
??
?
单步多值法（续）
?写成矩阵形式，就是,
?简记为,（ 10）
? Q阵就是 Y向量与 Z向量之间的变换阵。对 p
阶 Gear法,Q阵为（ p+1）阶非奇异方阵。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
1
4
1
2
1
4
3
4
1
2
3
4
7
!3
!2
1
2
0010
0001
3
2
k
k
k
k
k
h
k
h
k
k
y
y
f
y
y
y
yh
y
???
??
?
kk QYZ ?
单步多值法（续）
?将（ 10）式代入（ 4）及（ 5）式，可得到单
步多值法的计算公式,
? (11)
?说明：（ 11）式中计算 时以 为自
变量，但由于 是一个标量，所以实际上
只是用到它们的第一个分量 。
)()( ZZ
ZZ
Z
)1(
1
)(
1
)(
1
)(
1
)(
1
)1(
1
1)0(
1
?
?
?
?
?
?
?
??
??
?
?
???
??
?
?
?
?
i
k
i
k
i
k
i
k
i
k
i
k
kk
hfhf
QC
QB Q
G
G
Z
)( 1ikG? )1( 1)( 1,??? ikik ZZ
)( 1ikG?
)1( 1)( 1,??? ikik yy
单步多值法（续）
? 若记 （ 12）
? 则（ 11）式可简写为,
? （ 13）
? 在三阶的情况下,P及 L的值如下,
1
??
?
?
? ?
QCL
QB QP
)(
1
)(
1
)1(
1
)0(
1
??
???
??
?
??
?
?
?
i
k
i
k
i
k
kk
LGZZ
PZZ
? ?TP 111116116 1=L
1000
3100
3210
1111
?
?
?
?
?
?
?
?
?
?
?
?
?
6.3.1.4 误差估计与控制
? 当用单值多步的 Gear法对非线性病态系统进行仿真
时，它要求从初值开始，必须依靠显式法来启动，
即先从 及 开始，按一阶公式计算，然后逐次升
阶，计算出
? 在这种升阶过程中必须满足误差要求，对 k阶多步
法，其截断误差为,
? (14)
? 其中 是 t=ξ点上 y(t)的（ p+1）阶导数值，而 ξ
是所讨论区间中的某一个点。
0y 0y?
)p(y,,y,y,y 0000 ???????
)()1(11 ????? ppp yhCE
)()1( ??py
误差估计（续）
?例如，对三阶 Gear法,ξ就是 区间上
的某一个点。若设,则第 k步的截断
误差为,(15)
?为了估计,首先要估计,已知,
?现在要用 的差分来近似估计,即,
? ?htht kk ??,2
)1()1( )( ?? ? pkp yy ?
)1(11 ???? pkppk yhCE
kE )1( ?pky
TpkkkT)p(k!ph)(k!hk.kk )z,,z,z()y,,y,yh,y( pZ ?? 122 2 ??
)(pky )1( ?pky
h
yyy )p(k)p(k)p(
k 1
1 ?? ??
误差估计（续）
? 两边同乘以
? 则可得,
? （ 16）
? 将（ 16）式代入 (15)式可得,
? （ 17）
? 仿真中一般采用的是相对误差，若要求每一步的相对
误差不大于,即 （ 18）
?其中 为到目前为止已出现过的 y的最大绝对值
（注：若 y的初值为 0，则应取 ＝ 1为宜。
!
1
p
hp?
p
k
p
k
p
k
p
k
p
p
k
ppp
k
p
k
p
k
p
z
yphyphphh yypyh
z= z=
!!!!
1
)(
1
)(
1)(
1
)()1(1
??
????
?
?
?
?
??
pkpk pCE z!1 ?? ?
?0
0
m ax
1
z! ???
? ypC
p
k
p
maxy
maxy
误差估计（续）
?若系统为微分方程组，状态变量 y为 N个
（ y(1),y(2),…, y(N)），相对误差可定义
为,
? （ 19）
?若,则本步计算结果有效，进入下一
步；如果不满足，则需要减小步长或者采取
其它措施。
2
1
m a x
1
)(
)(z!
?
?
?
?
?
?
?
?
? ?? N
j
p
kp
k jy
jpC
E
0kE ??
误差控制技术
? 误差控制：一是改变步长，其二是改变阶次。
? 无论是 或,均需要考虑变阶或（与）变
步长，即通过改变 p及 h以使 。
?变阶或（与）变步长的原则,
?首先考虑仅改变步长 h（阶次不变），设新的
步长为,且 。那么,应取多大为
宜呢？
0??kE 0kE ???
0kE ??
sh hRh ps ? Rp
误差控制技术（续）
?为简便起见，我们根据单变量表达式（ 15）
来分析，即,
? 即
?根据（ 16）式,
?故有,
0
m a x
1)1(
1 )( ???
??
?
y
hRyC pppkp
0
m a x
1)1(
11 ???
??
??
y
hyCR ppkpp
p
pkpkp pyh z!)1(1 ????
0
m a x
11 z! ?????
y
pCR pkpp
p
误差控制技术（续）
?可得到,
?考虑到误差仅仅是估计值，考虑经验系数
1/1.2，这样,表达式如（ 20）式,
? （ 20）
? 若改变步长的效果不理想，则考虑要提升仿
真方法的阶次。
1
1
0
z?
?
?
?
p
max
p
k
p
p
y
!pC
R
?
pR
1
1
0
z21
1
?
?
?
?
p
m a x
p
k
p
p
y
!pC
.
R ?
误差控制技术（续）
?设当前为 p阶，考虑用 p＋ 1阶计算。假设此
时步长 能使相对误差接近规定的要求，
即,
?采用二阶差分来近似,即,
hRh ps 1??
0
222
12 ???
???
??
m a x
)p(
k
pp
pp
y
yhRC
)p(ky 2?
h
yyy )p(k)p(k)p(
k
1
1
1
2
?
?
?
? ??
误差控制技术（续）
? 考虑到,则不难得到,
? 也考虑经验系数（这里取其为 1/1.4，），则可得到升
阶时 的表达式如（ 21）式,
? （ 21）
? 如果 仿真步长缩短到相当小而其误差仍然达不到要求，
则可能是因阶次过高而稳定域达不到要求的缘故 。
p
kp
p
k h
py z!=
1
)1( ?
?
?
pkppkpkp)p(k h !phh !py z =zz= 22112 ???? ????
Rp?1
2 2
2
0
1 z41
1
?
?
? ??
p
m a x
p
k
p
p
y!pC
.R
?
误差控制技术（续）
?若将当前仿真的阶次 p降低一阶，仿真步长 h
变为,此时要使仿真误差接近规定要
求，则,
?考虑经验系数为 1/1.3，则可得到,
? (22)
?式中,它是向量 的最后一个分量。
hRh ps 1??
0
1 ??? ?
m a x
)p(
k
pp
pp
y
yhRC
p
m a x
p
k
p
p
y
!pC
.
R
z31
1 0
1
??
?
)(
!
p
k
p
p
k yp
hz ?
kZ
误差控制技术（续）
?步长发生变化时，如何由当前的 产生在
下的
?已知
?只变步长,
Zk sh
ksZ
Tp
kkk
Tp
kphkhkkk zzzyyyhy
pZ ),,,(=),,,,( 1)(!)2(!2,2 ???
ZZ k
p
p
pT)p(
k!p
)hR()(
k!
)hR(
k
.
pkks
R
R
)y,,y,yhR,y(
p
pp
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
1
2
2
2
误差控制技术（续）
?降阶变步长：从 p阶降为 p-1阶，步长由 h变
为,这时,,用 表示,hRp 1? kZ ksZ Z Zkp ksp? ?1 1,
TpkkkTpkphkhkkp
k zzzyyyhy
pZ ),,,(=),,,,( 11)1(
)!1()2(!2
1 12 ??
?
? ?? ???
ZZ
p
k
p
p
pTp
kp
hR
k
hR
kpk
p
ks
R
R
yyyhRy
p
pp 1
1
1
1)1(
)!1(
)()2(
!2
)(
1
1
1
= ),,,,( 1
2
1 ?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ??
?
??
误差控制技术（续）
?升阶变步长：当由 p阶升到 p+1阶，步长由 h
变为,则这时,,用 表示,
?然而,是未知的，为此，必须由已知的
数据来估计。若采用差分法，即,
R hp?1 kZ ksZ Z Zkp ksp? ?1 1,
TpkpkkkTpkphpkp hRkhkkp
k zzzzyyyyhy
pppZ ),,,,(=),,,,( 11)1(
)!1()(!
)()2(
!2
1 12 ??
?
? ?? ???
Tpkp hRpkp hRkhRkpkp
ks yyyyhRy
pppppZ ),,,,,( )1(
)!1(
)()(
!
)()2(
!2
)(
1
1 1121 ?
??
? ???? ??
ykp( )?1
y y yhk p k
p
k
p
( )
( ) ( )
? ?? ?1 1
误差控制技术（续）
? 则得到 的表达式如下,z
ksp?1
ZZ
p
k
p
p
pTp
kp
hR
k
hR
kpk
p
ks
R
R
yyyhRy
p
pp 1
1
1
1)1(
)!1(
)()2(
!2
)(
1
1
1
=),,,,( 1
2
1 ?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ??
?
??
6.3.1.5 加速收敛问题
? Gear法用于病态系统进行仿真时，必须先用显式公
式启动，然后用隐式法进行校正，经过多次迭代，
达到适当精度后，再往前推进。显然，迭代的收敛
速度极大地影响着仿真效能。
? 三阶 Gear法为例
? (1)
? 其中 就是 的简写，它一般是 的非线
性函数，（ 1）式是一非线性方程。
)62918( 1211111 ???? ???? kkkkk hfyyyy
fk?1 f yk( )?1 yk?1
加速收敛问题 (续）
?迭代计算时，我们采用的迭代公式为,
? （ 2）
?其中 是向量 L的第一个分量，这称为
Picard迭代。对该迭代过程,
? （ 3）
)ff(hly
)ff(hyy
)i(
k
)i(
k
)i(
k
)i(
k
)i(
k
)i(
k
)i(
k
1
1101
1
1111
6
1
1
1
= ????
?
???
?
?
??
???
l0
?
?
?
?
?
?
?
??
?
?
??
?
?
?
??? ???
)1,2,=(
)0(
1+0
)(
1+
)(
1+01
1
1
1)+(
1+
?i
y
f
hl
y
f
hlyy
k
)i(
i
k
i
k
)i(
k
)i(
k
i
k
??
??
?
?
?
?
加速收敛问题 (续）
?为保证收敛，则要求,即,
?一般 大约为 1～ 3之间，可见,h不能太大，
否则迭代将不收敛或收敛速度十分缓慢。这
大大限制了 Gear法的有效性。
?牛顿迭代法,
?（ 1）式可改写为,
? （ 4）
l h fy0 1?? ? h
f
y l
?
?
?
1
0
1 0/l
Cyhfl
yyyyhfly
i
k
kkk
i
k
i
k
?
????
?
???
?
?
)(=
)2918()(
)(
10
2111
1)(
10
)1(
1
加速收敛问题 (续）
?简记为,（ 5）
?令 （ 6）
?显然，（ 5）式与 g(y)=0同解。由牛顿迭代法，
不难得到,
? （ 7）
cyhfly 0 ?? )(
cyhflyyg 0 ???? )()(
)i(yy
)i(
)i()i(
y
g
)y(gyy
?
?
???
?
???
?
??
1
?
?
加速收敛问题 (续）
?从 (6)式求得 的表达式，并将 代入
及（ 7）式，可得,
? (8)
?将 代入 (5)式，所得结果在代入 (8)式，
可得,
?
?
g
y y y i? ( )
?? gy
)i(yy
)i()i(
)i()i(
y
f
hl
c)y(hfly
yy
?
?
??
?
?
??
?
?
?
???
??
1 0
01
?
?
y y i? ( )
加速收敛问题 (续）
? (9)
)i(yy
)i()i(
)i()i(
y
fhl
c)y(hfl)c)y(hfl(yy
?
?
?
???
?
???
?
?
??????
1
0
0
1
01
?
?
? ?
1
()(
=
)(
0
)1()(
0)(
iyy
ii
i
y
fhl
yhfyhfl
y
?
?
??
?
?
??
?
?
?
?
?
?
?
1
1+ = (
0
0
G )i
yy
)i( l
y
fhl
y
)i(?
???
?
???
??
?
?
加速收敛问题 (续）
?多了一个“加速”因子,
?对微分方程组来说,y是 N维向量,(9)式成为
如下形式,
? (10)
?其中 Y,F,G为 N维向量,I为 N阶单位阵，
而 为 N× N的方阵，亦称为雅可比矩
阵,
1
1 0?
?
??
?
?? ?l h
f
y y y i
?
? ( )
Y (0
1
0
1 Gyy )i
YY
)i()i( lFhlI
)i(
?
?
?
??
?
??
? ?
?
??
?
????
?
?
? ? F Y/
加速收敛问题 (续）
?
? (11)
?采用牛顿法,
? (12)
?其中 Z为（ K＋ 1） × 1的向量，而 G为标量。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
N
NNN
N
N
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
Y
F
21
2
2
2
1
2
1
2
1
1
1
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
??
?
?
???
?
?
?
?
?
?
?
)i(
k
ik
)i(
k
)i(
k
k
)(
k
G
y
f
hlLZZ
PZZ
1
1
01
1
1
0
1
1
?
?
加速收敛问题 (续）
? 对于 y为 N维的情况，则有,
? (13)
? 式中的 Z为（ K＋ 1） × N的矩阵，而 G为 N× 1的向
量。
? 牛顿法能加快迭代过程的收敛，这是该方法的优
点，然而，因每次迭代时都必须计算雅可比阵，
还要计算矩阵求逆，计算量将增加较多。
? ?
?
?
?
??
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
???
?
?
?
?
?
?
?
Ti
1k
1
i
0
i
1k
1i
1k
k
0
1k
ZZ
ZZ
G
Y
F
hlILZ
P
)( )()()(
)(
6.3.1.6 病态性探测
? 系统并不是每一步均处于病态。从数值计算的角度
进行定义,病态性”和“非病态性”,
? 如果用非 Stiff法仿真，其步长限制是由于计算精
度原因引起的而不是由于稳定性引起的，则称系统
呈现,非病态性,，反之，若步长是因稳定性而受
到限制时，则系统呈现,病态性,。
? 当系统呈现病态性时应采用 Stiff法，而当系统呈
现非病态性时应采用非 Stiff法。
病态性探测 (续）
? 病态性探测的三种方法,
? 1） 稳定半径法
? 非 Stiff法的稳定域的形状接近半圆，其稳定域与实
轴的交点称为稳定半径（ ）。若记 ρ为系统的最
大特征值，为保证计算稳定，则要求,
? (1)
? 各种方法的稳定半径已知，而 ρ的值可由雅可比阵
的范数 来估计，即,
? 因此，在用 stiff方法对病态系统仿真过程中，雅可
比矩阵的值可以得到，则在下一步计算时，先用
（ 1）式进行判断病态性，再决定是否仍采用 stiff方
法仿真。
rk
krh ???
J ?? J
病态性探测 (续）
?2） 嵌入低阶大稳定域法
?如果当前采用的是非 Stiff法，而计算误差不
能满足要求，为判断其稳定性，可降低阶次，
如果此时误差减少，则说明当前已处于非稳
定状态（非 Stiff法的低阶的稳定域大于高阶
的稳定域）。因而下一步应采用 Stiff法计算。
病态性探测 (续）
? 3） 小稳定域试探法
? 嵌入低阶大稳定域法对于高阶的非 Stiff法是有效的，
然而，如果当前使用的已经是低阶方法，则很难再
嵌入一个更低阶的算法。小稳定域试探法采用与嵌
入低阶大稳定域法相反的思路来进行病态性探测。
? 若当前采用的是非 Stiff法，为探测下一步的病态
性，步长不变，用比其高一阶的方法进行探测，若
误差加大，说明当前已处于病态或病态边缘，因此
下一步要么减少步长，要么采用 Stiff方法计算。
第 6章 病态系统仿真
剡昌锋 刘军
兰州理工大学机电工程学院
6.1病态系统的定义
?系统中各环节的时间常数差异巨大。为保证
仿真计算的稳定性，由于仿真步长必须限制
在最小时间常数的数量级而选得很小，然而
仿真结束的时间则决定于系统中的最大时间
常数，若按满足稳定性要求所选择的步长进
行仿真，则不仅整个仿真所花费的时间非常
长，甚至由于计算的舍入误差而导致整个仿
真的失败。这就是所谓,病态系统仿真,问
题。
病态系统的定义（续）
? 病态系统定义, （ 1）
? 令
? （ 2）
? J称为系统的雅可比矩阵。
? 若 J的特征值全部具有负实部，且有：,
则该系统称为病态系统，在某些文献中也叫做刚性
系统 (Stiff)，而,
? 称为病态比，一般在 50以上。
),( tyFy ??
1
1
1
1
J
y
F
n
nn
n
y
f
y
f
y
f
y
f
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
?
ii ?? m a xm in ??
i
iS
?
?
min
max?
6.2 线性病态系统仿真
?对线性定常系统，我们可用如下状态方程进
行一般性描述,
? (3)
?1,增广矩阵法 将 u(t)作为系统的增广状态
?
?
?
?
?
?
?
?
??
)0()(
)()(
)()(
0 xtx
tCxty
tButAxx?
?
?
?
??
?
?
?
?
)0(~)(~
)(~C
~
=( t )
)(~
~~
0 xtx
txy
txAx?
)(~)(~)1(~ kxtkx ???
线性病态系统仿真 (续）
? 其中,
? 由于病态系统特征值相差倍数很大，必须
用加速收敛的方法计算该状态转移矩阵的
值。
? 2.蛙跳算法 基本思想,
? （ 1）考虑作用函数为阶跃函数，则增广状
态十分简单。
? 选择 q，使得,
e hA~~ =(T )?
hqhq nn 2,2 ??
线性病态系统仿真 (续）
?(3) 仿真计算时采用如下“蛙跳”方式,
??
??
)2()3(
)()2(
)0()0()()2(
)0()0()4(
)0()0()2(
)0()(
~
2
~
2
~
2
~
2
~
2
~
4
~
2
~
2
~
2
~~
~
11
qhxeqhx
qhxeqhx
xexeeqhxhx
xexeehx
xexeehx
xehx
hA
hA
hAhAhAn
hAhAhA
hAhAhA
hA
n
n
nnn
?
?
???
?
?
????
???
???
?
??
线性病态系统仿真 (续）
? 即在 qh以前采用加倍跳跃式计算，而在 qh以后每隔
qh计算一次。
? 优点是,h可以取得很小（可按最小时间常数考
虑），从而保证初始阶段的精度而计算量却不大，
而到 qh以后，小时间常数的作用完成，则加大步长
计算，从而加快仿真计算速度。
? 从 x(h)到 x(qh)都是以 x(0)为基础进行计算，所以误
差传播比较均匀（仅仅是状态转移矩阵的误差）。
6.3 非线性病态系统仿真
?一般非线性系统的仿真大多采用数值积分法。
而数值积分法一般又只具有有限的稳定域，
典型的如 龙格－库塔法,仿真步长限定在系
统最小时间常数的数量级，才能保证计算的
稳定性，而系统的过渡过程时间却决定于最
大时间常数，因而对病态系统来说计算量极
大，加上存在误差传播，仿真的精度甚至稳
定性也会受到影响。
6.3.1 吉尔 (Gear)法
? 6.3.1.1.Stiff稳定域
? Gear研究后发现，并不要求一定采用恒稳方法，
而只要具有所谓 Stiff稳定域就可以了。
? Stiff稳定域定义：对实际的物理系统，时间常数 μ
一般小于零。选择仿真步长 h若满足,
? 与
? 可保证仿真的稳定性，称该算法具有 stiff稳定域。
(I) Im 0Re
??
?
?
?
??
?
h
h
( I I ) Re ?? ??h
Stiff稳定域（续）
? 实际上具有 Stiff稳定域的方法
与恒稳方法只在近虚轴处有一
点差别，即如果系统中的极点
全部为实极点，那么无论选择
多大的步长，计算是恒稳的。
如果系统中有复极点（实部仍
为负数），只要步长的选择满
足上述条件，也能保证算法稳
定。
0
θ
-θ
Reμ h
I
II
-δ
图 6.1 Stiff稳定域
Stiff稳定域（续）
?Stiff域中 θ与 δ的确定,
?按病态系统的大特征值来选择步长,
?该特征值所对应的模态大约要经过 4倍左右时
间常数的时间才能有效地衰减掉，即
?,也就是
?这样，此时即使加大步长 h，也能保持计算的
稳定性。基于这一考虑，可设 δ＝ - 4。
m a xm a xm a x ??? j??
m a x/4 ??h 4max ??h
Stiff稳定域（续）
? 另一方面，考虑到系统特征值为复数，它所对应的
瞬态响应呈振荡型。一个振荡周期内至少计算 N个点。
最小振荡周期为,
? 其中 h为计算步长，若选择 N≥8，则有：,因
此可选 θ＝ π/4。
? 综上所述，如果选择某一种方法，其稳定域 θ≥π/4，
且 |δ|≤4，则从使用的角度来看，图 6.1所示的稳定域
与恒稳域没有差别，从而完全可以用于病态系统的
仿真。
NhT ?? m a xm i n /2 ??
4/ma x ?? ?h
6.3.1.2吉尔（ Gear）法的基本原理
? 设系统,
?满足 Stiff稳定域的多步法,Gear提出的用于
病态系统仿真的计算公式是,
? (1)
?
00 )( ),(/ yttytyfdtdyy ?????
?
?
???? ??
k
j
knjknjkn yhyay
1
0 ??
用于病态系统仿真的 Gear公式的系数表
α1 α2 α3 α4 α5 α6 β0
G-1 1 1
G-2 4/3 -1/3 2/3
G-3 18/11 -9/11 2/11 6/11
G-4 48/25 -36/25 16/25 -3/25 12/25
G-5 300/137 -300/137 200/137 -75/137 12/137 60/137
G-6 360/147 -450/147 400/147 -225/147 72/147 -10/147 60/147
吉尔（ Gear）法的基本原理 (续）
? 稳定域如图 6.2所示。
? 从图上可以看出，该方
法在 5阶以下（包括 5阶）
的稳定域满足 Stiff稳定
域的条件
（ θ>π/4,|δ|<4）,
? 大于 5阶时 |δ|>4，而且
还可能穿过负实轴。 -10 -8 -6 -4 -2 0 2 4 6 8 10
0
2
4
6
8
10
图 6.2 用于病态系统仿真的 G e a r 法
C
2
C
3
C
4
C
5
C
6
吉尔（ Gear）法的基本原理 (续）
? 在用 Gear法仿真非线性病态系统时，有以下三个基本问题需要解决,
? 1） 启动问题
? 上述 Gear法本质上是隐式多步法。
? 对于初值问题，困难：隐式方法一般用显式方法启动，即先进行预报，
然后通过迭代进行校正。如果迭代方法的收敛性不好，可能引起计算发散或计算量加大。
? 即使选择的迭代方法收敛性满足要求，显式多步法预报，仍然难以启
动，必须采用单步法启动，由于单步法不具有 Stiff稳定域，因而很难
保证计算的稳定性。
? 2） 变步长策略
? 非线性病态系统仿真往往采用变步长策略，如何适时地将步长调整到
合适长度，以同时满足仿真精度和速度的要求。
? 3） 加速迭代 为了提高计算效能，加速迭代也是非线性病态系统仿
真中重要问题。
6.3.1.3 单步多值法
? 以三阶为例，采用显式多步法进行预报，然后用隐
式法校正。
? 其显式预报的公式是,（ 2）
?三阶隐式 Gear公式校正,
? （ 3）
kkkkk hfyyyy 33 221123)0( 1 ????? ???
?
)1(
111
6
211
2
111
9
11
18)2(
1
)0(
111
6
211
2
111
9
11
18)1(
1
????
????
????
????
kkkkk
kkkkk
hfyyyy
hfyyyy
单步多值法（续）
?其中等式右边第 4项为导函数项，它是通过将
第 i次迭代所得到的 y的预报值代入导函数后
计算得到的。为便于程序实现，由（ 2）式，
并令,
24
5
14
23
2
11)1(
1
)1(
1
)0(
111
6
22
1
12
3
)1(
1
)0(
111
6)0(
1
)0(
111
6
211
2
111
9
11
18)1(
1
7
)(33-=
)(
??
?
?
?
????
?
???
????
????
?????
???
????
kkkkk
kkkkkk
kkk
kkkkk
yyyhff
hfhfhfyyy
hfhfy
hfyyyy
单步多值法（续）
?迭代的校正公式可表示成,
?
?
)(
)(
)(
)1(
1
)(
111
6)(
1
)1(
1
)0(
1
)1(
111
6)1(
1
)2(
1
)1(
1
)0(
111
6)0(
1
)1(
1
?
???
?
?
????
?
????
???
???
???
i
k
i
k
i
k
i
k
kkkk
kkkk
hfhfyy
hfhfyy
hfhfyy
单步多值法（续）
?对一般情形，令,
Tpkkkkkk yyyhfyY ),,,,,( 121 ????? ?
Tpkkkkkk yyyhfyY ),,,,,( 21111 ?????? ? ?
),(),(= 1)1( 11)( 1
)1(
1
)(
1
)(
1
?
?
???
?
???
?
??
k
i
kk
i
k
i
k
i
k
i
k
tyhftyhf
hfhfG
单步多值法（续）
? 其中 p表示 Gear法的阶次。对三阶预报－迭代校正
算法,
? 用矩阵的方法加以表示,
? （ 4）
? （ 5）
? 其中 B,C分别是相应的系数矩阵。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
1
4
5
2
11
4
23
2
1
2
3
1
)1(
1
)0(
1
0100
0001
7
33
k
k
k
k
k
k
k
k
y
y
hf
y
y
y
hf
y
)(
1
11
6
1
)1(
1
)(
1
1
)(
1
)1(
1
0
0
1 i
k
k
k
i
k
i
k
k
k
i
k
i
k
G
y
y
hf
y
y
y
hf
y
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
),.,.,k -,( i =YBY k)(k 121 0 1 ??
)121( 1111,.,.,k -,i=GCYY )i(k)i(k)i(k ???? ??
单步多值法（续）
?然而，对初值问题，上式是不能自启动的。
?解决方法 ------单步多值法：用高阶导数值来
取代前几步的 y及 f的值。
? 先定义一个向量，称之为 Nordsieck向量,
? （ 6）
? （ 7）
?需要确定 Z向量与 Y向量之间的关系。
),,,,( )(!)2(!2 2 Tpkphkhkkk yyyhy pZ ???
),,,,( )( 1!)2(!2111 2 Tpkphkhkkk yyyhy pZ ???? ? ??
单步多值法（续）
?采用多项式逼近，以三阶为例,
? （ 8）
?在 处，有
? （ 9）
6)(
)(62)(
)(3)(2)(
)()()()(
3
32
2
321
3
3
2
210
?
?
?
?
?
?
?
?
???
?????
???????
aty
ttaaty
ttattaaty
ttattattaaty
k
kk
kkk
???
??
?
单步多值法（续）
?同样，也可以得到,
?由（ 8）式及（ 9）式消去,整理后
可得到,
!3,!2,,3210 ayayayay kkkk ???? ???
3210,,,aaaa
?
?
?
?
?
?
?
?????
?????
?
?
??
??
24
1
12
13
4
3
!3
24
1
12
3
4
7
!2
3
2
2
kkkkk
h
kkkkk
h
kk
kk
yyhfyy
yyhfyy
hfyh
yy
???
??
?
单步多值法（续）
?写成矩阵形式，就是,
?简记为,（ 10）
? Q阵就是 Y向量与 Z向量之间的变换阵。对 p
阶 Gear法,Q阵为（ p+1）阶非奇异方阵。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
2
1
4
1
2
1
4
3
4
1
2
3
4
7
!3
!2
1
2
0010
0001
3
2
k
k
k
k
k
h
k
h
k
k
y
y
f
y
y
y
yh
y
???
??
?
kk QYZ ?
单步多值法（续）
?将（ 10）式代入（ 4）及（ 5）式，可得到单
步多值法的计算公式,
? (11)
?说明：（ 11）式中计算 时以 为自
变量，但由于 是一个标量，所以实际上
只是用到它们的第一个分量 。
)()( ZZ
ZZ
Z
)1(
1
)(
1
)(
1
)(
1
)(
1
)1(
1
1)0(
1
?
?
?
?
?
?
?
??
??
?
?
???
??
?
?
?
?
i
k
i
k
i
k
i
k
i
k
i
k
kk
hfhf
QC
QB Q
G
G
Z
)( 1ikG? )1( 1)( 1,??? ikik ZZ
)( 1ikG?
)1( 1)( 1,??? ikik yy
单步多值法（续）
? 若记 （ 12）
? 则（ 11）式可简写为,
? （ 13）
? 在三阶的情况下,P及 L的值如下,
1
??
?
?
? ?
QCL
QB QP
)(
1
)(
1
)1(
1
)0(
1
??
???
??
?
??
?
?
?
i
k
i
k
i
k
kk
LGZZ
PZZ
? ?TP 111116116 1=L
1000
3100
3210
1111
?
?
?
?
?
?
?
?
?
?
?
?
?
6.3.1.4 误差估计与控制
? 当用单值多步的 Gear法对非线性病态系统进行仿真
时，它要求从初值开始，必须依靠显式法来启动，
即先从 及 开始，按一阶公式计算，然后逐次升
阶，计算出
? 在这种升阶过程中必须满足误差要求，对 k阶多步
法，其截断误差为,
? (14)
? 其中 是 t=ξ点上 y(t)的（ p+1）阶导数值，而 ξ
是所讨论区间中的某一个点。
0y 0y?
)p(y,,y,y,y 0000 ???????
)()1(11 ????? ppp yhCE
)()1( ??py
误差估计（续）
?例如，对三阶 Gear法,ξ就是 区间上
的某一个点。若设,则第 k步的截断
误差为,(15)
?为了估计,首先要估计,已知,
?现在要用 的差分来近似估计,即,
? ?htht kk ??,2
)1()1( )( ?? ? pkp yy ?
)1(11 ???? pkppk yhCE
kE )1( ?pky
TpkkkT)p(k!ph)(k!hk.kk )z,,z,z()y,,y,yh,y( pZ ?? 122 2 ??
)(pky )1( ?pky
h
yyy )p(k)p(k)p(
k 1
1 ?? ??
误差估计（续）
? 两边同乘以
? 则可得,
? （ 16）
? 将（ 16）式代入 (15)式可得,
? （ 17）
? 仿真中一般采用的是相对误差，若要求每一步的相对
误差不大于,即 （ 18）
?其中 为到目前为止已出现过的 y的最大绝对值
（注：若 y的初值为 0，则应取 ＝ 1为宜。
!
1
p
hp?
p
k
p
k
p
k
p
k
p
p
k
ppp
k
p
k
p
k
p
z
yphyphphh yypyh
z= z=
!!!!
1
)(
1
)(
1)(
1
)()1(1
??
????
?
?
?
?
??
pkpk pCE z!1 ?? ?
?0
0
m ax
1
z! ???
? ypC
p
k
p
maxy
maxy
误差估计（续）
?若系统为微分方程组，状态变量 y为 N个
（ y(1),y(2),…, y(N)），相对误差可定义
为,
? （ 19）
?若,则本步计算结果有效，进入下一
步；如果不满足，则需要减小步长或者采取
其它措施。
2
1
m a x
1
)(
)(z!
?
?
?
?
?
?
?
?
? ?? N
j
p
kp
k jy
jpC
E
0kE ??
误差控制技术
? 误差控制：一是改变步长，其二是改变阶次。
? 无论是 或,均需要考虑变阶或（与）变
步长，即通过改变 p及 h以使 。
?变阶或（与）变步长的原则,
?首先考虑仅改变步长 h（阶次不变），设新的
步长为,且 。那么,应取多大为
宜呢？
0??kE 0kE ???
0kE ??
sh hRh ps ? Rp
误差控制技术（续）
?为简便起见，我们根据单变量表达式（ 15）
来分析，即,
? 即
?根据（ 16）式,
?故有,
0
m a x
1)1(
1 )( ???
??
?
y
hRyC pppkp
0
m a x
1)1(
11 ???
??
??
y
hyCR ppkpp
p
pkpkp pyh z!)1(1 ????
0
m a x
11 z! ?????
y
pCR pkpp
p
误差控制技术（续）
?可得到,
?考虑到误差仅仅是估计值，考虑经验系数
1/1.2，这样,表达式如（ 20）式,
? （ 20）
? 若改变步长的效果不理想，则考虑要提升仿
真方法的阶次。
1
1
0
z?
?
?
?
p
max
p
k
p
p
y
!pC
R
?
pR
1
1
0
z21
1
?
?
?
?
p
m a x
p
k
p
p
y
!pC
.
R ?
误差控制技术（续）
?设当前为 p阶，考虑用 p＋ 1阶计算。假设此
时步长 能使相对误差接近规定的要求，
即,
?采用二阶差分来近似,即,
hRh ps 1??
0
222
12 ???
???
??
m a x
)p(
k
pp
pp
y
yhRC
)p(ky 2?
h
yyy )p(k)p(k)p(
k
1
1
1
2
?
?
?
? ??
误差控制技术（续）
? 考虑到,则不难得到,
? 也考虑经验系数（这里取其为 1/1.4，），则可得到升
阶时 的表达式如（ 21）式,
? （ 21）
? 如果 仿真步长缩短到相当小而其误差仍然达不到要求，
则可能是因阶次过高而稳定域达不到要求的缘故 。
p
kp
p
k h
py z!=
1
)1( ?
?
?
pkppkpkp)p(k h !phh !py z =zz= 22112 ???? ????
Rp?1
2 2
2
0
1 z41
1
?
?
? ??
p
m a x
p
k
p
p
y!pC
.R
?
误差控制技术（续）
?若将当前仿真的阶次 p降低一阶，仿真步长 h
变为,此时要使仿真误差接近规定要
求，则,
?考虑经验系数为 1/1.3，则可得到,
? (22)
?式中,它是向量 的最后一个分量。
hRh ps 1??
0
1 ??? ?
m a x
)p(
k
pp
pp
y
yhRC
p
m a x
p
k
p
p
y
!pC
.
R
z31
1 0
1
??
?
)(
!
p
k
p
p
k yp
hz ?
kZ
误差控制技术（续）
?步长发生变化时，如何由当前的 产生在
下的
?已知
?只变步长,
Zk sh
ksZ
Tp
kkk
Tp
kphkhkkk zzzyyyhy
pZ ),,,(=),,,,( 1)(!)2(!2,2 ???
ZZ k
p
p
pT)p(
k!p
)hR()(
k!
)hR(
k
.
pkks
R
R
)y,,y,yhR,y(
p
pp
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
1
2
2
2
误差控制技术（续）
?降阶变步长：从 p阶降为 p-1阶，步长由 h变
为,这时,,用 表示,hRp 1? kZ ksZ Z Zkp ksp? ?1 1,
TpkkkTpkphkhkkp
k zzzyyyhy
pZ ),,,(=),,,,( 11)1(
)!1()2(!2
1 12 ??
?
? ?? ???
ZZ
p
k
p
p
pTp
kp
hR
k
hR
kpk
p
ks
R
R
yyyhRy
p
pp 1
1
1
1)1(
)!1(
)()2(
!2
)(
1
1
1
= ),,,,( 1
2
1 ?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ??
?
??
误差控制技术（续）
?升阶变步长：当由 p阶升到 p+1阶，步长由 h
变为,则这时,,用 表示,
?然而,是未知的，为此，必须由已知的
数据来估计。若采用差分法，即,
R hp?1 kZ ksZ Z Zkp ksp? ?1 1,
TpkpkkkTpkphpkp hRkhkkp
k zzzzyyyyhy
pppZ ),,,,(=),,,,( 11)1(
)!1()(!
)()2(
!2
1 12 ??
?
? ?? ???
Tpkp hRpkp hRkhRkpkp
ks yyyyhRy
pppppZ ),,,,,( )1(
)!1(
)()(
!
)()2(
!2
)(
1
1 1121 ?
??
? ???? ??
ykp( )?1
y y yhk p k
p
k
p
( )
( ) ( )
? ?? ?1 1
误差控制技术（续）
? 则得到 的表达式如下,z
ksp?1
ZZ
p
k
p
p
pTp
kp
hR
k
hR
kpk
p
ks
R
R
yyyhRy
p
pp 1
1
1
1)1(
)!1(
)()2(
!2
)(
1
1
1
=),,,,( 1
2
1 ?
?
?
??
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
? ??
?
??
6.3.1.5 加速收敛问题
? Gear法用于病态系统进行仿真时，必须先用显式公
式启动，然后用隐式法进行校正，经过多次迭代，
达到适当精度后，再往前推进。显然，迭代的收敛
速度极大地影响着仿真效能。
? 三阶 Gear法为例
? (1)
? 其中 就是 的简写，它一般是 的非线
性函数，（ 1）式是一非线性方程。
)62918( 1211111 ???? ???? kkkkk hfyyyy
fk?1 f yk( )?1 yk?1
加速收敛问题 (续）
?迭代计算时，我们采用的迭代公式为,
? （ 2）
?其中 是向量 L的第一个分量，这称为
Picard迭代。对该迭代过程,
? （ 3）
)ff(hly
)ff(hyy
)i(
k
)i(
k
)i(
k
)i(
k
)i(
k
)i(
k
)i(
k
1
1101
1
1111
6
1
1
1
= ????
?
???
?
?
??
???
l0
?
?
?
?
?
?
?
??
?
?
??
?
?
?
??? ???
)1,2,=(
)0(
1+0
)(
1+
)(
1+01
1
1
1)+(
1+
?i
y
f
hl
y
f
hlyy
k
)i(
i
k
i
k
)i(
k
)i(
k
i
k
??
??
?
?
?
?
加速收敛问题 (续）
?为保证收敛，则要求,即,
?一般 大约为 1～ 3之间，可见,h不能太大，
否则迭代将不收敛或收敛速度十分缓慢。这
大大限制了 Gear法的有效性。
?牛顿迭代法,
?（ 1）式可改写为,
? （ 4）
l h fy0 1?? ? h
f
y l
?
?
?
1
0
1 0/l
Cyhfl
yyyyhfly
i
k
kkk
i
k
i
k
?
????
?
???
?
?
)(=
)2918()(
)(
10
2111
1)(
10
)1(
1
加速收敛问题 (续）
?简记为,（ 5）
?令 （ 6）
?显然，（ 5）式与 g(y)=0同解。由牛顿迭代法，
不难得到,
? （ 7）
cyhfly 0 ?? )(
cyhflyyg 0 ???? )()(
)i(yy
)i(
)i()i(
y
g
)y(gyy
?
?
???
?
???
?
??
1
?
?
加速收敛问题 (续）
?从 (6)式求得 的表达式，并将 代入
及（ 7）式，可得,
? (8)
?将 代入 (5)式，所得结果在代入 (8)式，
可得,
?
?
g
y y y i? ( )
?? gy
)i(yy
)i()i(
)i()i(
y
f
hl
c)y(hfly
yy
?
?
??
?
?
??
?
?
?
???
??
1 0
01
?
?
y y i? ( )
加速收敛问题 (续）
? (9)
)i(yy
)i()i(
)i()i(
y
fhl
c)y(hfl)c)y(hfl(yy
?
?
?
???
?
???
?
?
??????
1
0
0
1
01
?
?
? ?
1
()(
=
)(
0
)1()(
0)(
iyy
ii
i
y
fhl
yhfyhfl
y
?
?
??
?
?
??
?
?
?
?
?
?
?
1
1+ = (
0
0
G )i
yy
)i( l
y
fhl
y
)i(?
???
?
???
??
?
?
加速收敛问题 (续）
?多了一个“加速”因子,
?对微分方程组来说,y是 N维向量,(9)式成为
如下形式,
? (10)
?其中 Y,F,G为 N维向量,I为 N阶单位阵，
而 为 N× N的方阵，亦称为雅可比矩
阵,
1
1 0?
?
??
?
?? ?l h
f
y y y i
?
? ( )
Y (0
1
0
1 Gyy )i
YY
)i()i( lFhlI
)i(
?
?
?
??
?
??
? ?
?
??
?
????
?
?
? ? F Y/
加速收敛问题 (续）
?
? (11)
?采用牛顿法,
? (12)
?其中 Z为（ K＋ 1） × 1的向量，而 G为标量。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
N
NNN
N
N
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
y
f
Y
F
21
2
2
2
1
2
1
2
1
1
1
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
????
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
??
?
?
???
?
?
?
?
?
?
?
)i(
k
ik
)i(
k
)i(
k
k
)(
k
G
y
f
hlLZZ
PZZ
1
1
01
1
1
0
1
1
?
?
加速收敛问题 (续）
? 对于 y为 N维的情况，则有,
? (13)
? 式中的 Z为（ K＋ 1） × N的矩阵，而 G为 N× 1的向
量。
? 牛顿法能加快迭代过程的收敛，这是该方法的优
点，然而，因每次迭代时都必须计算雅可比阵，
还要计算矩阵求逆，计算量将增加较多。
? ?
?
?
?
??
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
???
?
?
?
?
?
?
?
Ti
1k
1
i
0
i
1k
1i
1k
k
0
1k
ZZ
ZZ
G
Y
F
hlILZ
P
)( )()()(
)(
6.3.1.6 病态性探测
? 系统并不是每一步均处于病态。从数值计算的角度
进行定义,病态性”和“非病态性”,
? 如果用非 Stiff法仿真，其步长限制是由于计算精
度原因引起的而不是由于稳定性引起的，则称系统
呈现,非病态性,，反之，若步长是因稳定性而受
到限制时，则系统呈现,病态性,。
? 当系统呈现病态性时应采用 Stiff法，而当系统呈
现非病态性时应采用非 Stiff法。
病态性探测 (续）
? 病态性探测的三种方法,
? 1） 稳定半径法
? 非 Stiff法的稳定域的形状接近半圆，其稳定域与实
轴的交点称为稳定半径（ ）。若记 ρ为系统的最
大特征值，为保证计算稳定，则要求,
? (1)
? 各种方法的稳定半径已知，而 ρ的值可由雅可比阵
的范数 来估计，即,
? 因此，在用 stiff方法对病态系统仿真过程中，雅可
比矩阵的值可以得到，则在下一步计算时，先用
（ 1）式进行判断病态性，再决定是否仍采用 stiff方
法仿真。
rk
krh ???
J ?? J
病态性探测 (续）
?2） 嵌入低阶大稳定域法
?如果当前采用的是非 Stiff法，而计算误差不
能满足要求，为判断其稳定性，可降低阶次，
如果此时误差减少，则说明当前已处于非稳
定状态（非 Stiff法的低阶的稳定域大于高阶
的稳定域）。因而下一步应采用 Stiff法计算。
病态性探测 (续）
? 3） 小稳定域试探法
? 嵌入低阶大稳定域法对于高阶的非 Stiff法是有效的，
然而，如果当前使用的已经是低阶方法，则很难再
嵌入一个更低阶的算法。小稳定域试探法采用与嵌
入低阶大稳定域法相反的思路来进行病态性探测。
? 若当前采用的是非 Stiff法，为探测下一步的病态
性，步长不变，用比其高一阶的方法进行探测，若
误差加大，说明当前已处于病态或病态边缘，因此
下一步要么减少步长，要么采用 Stiff方法计算。
