系统仿真技术
第 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方法计算。