2.2 连续系统仿真的数值积分法

从控制理论中可知，对于一个连续时间系统可以在时域、频域中描述其动态特性。然而，在工程实际和科学研究中所遇到的实际问题往往很复杂，在很多情况下都不可能给出描述系统动态特性的数学模型的解析表达式，多数只能用近似的数值方法求解。随着计算机硬件、软件的发展和数值理论的进展，微分方程的数值解方法已成为当今研究、分析、设计动态系统的一种有力工具。

2.2.1 基本的数值积分公式

一个n阶连续系统可以用n阶微分方程表示，也可以用n个一阶微分方程组来描述，而每个一阶微分方程的解实际上是一个积分，利用数字计算机来进行连续系统仿真，从本质上讲，就是要在数字计算机上构造出几个数字积分器，也就是让数字计算机进行数值积分运算。

1.基本思想

在连续系统仿真中，主要的数值计算工作是对考虑初值问题的方程进行求解

对式（2.11）中所示的初值问题的解y(t)是一连续变量t的函数。现在要以一系列离散时刻近似值y(t1),y(t2),…,y(tn)来代替，其中ti=t0+ih（h称为步长，是相邻两点之间的距离）。称点列yk=y(tk)(k=0, 1, …, n)为式（2.11）在点列tk(k=0, 1, …, n)处的数值解。可见，数值解法就是首先把一个连续变量问题，用某种离散化方法化成离散变量问题——近似的差分方程的初值问题，然后逐步计算出yk。采用不同的离散化方法可得出具有不同精度的微分方程的数值积分方法。

由于

若令

则有

y(tn+1)≈yn+Qn （2.14）

因此，主要问题是如何对Qn进行求解，即如何对f(y,t)进行近似积分。

2.欧拉法

欧拉法是最简单的数值积分法，虽然它的精度较差，但由于公式简单，而且具有明显的几何意义，有利于初学者在直观上学习数值y(tn)是怎样逼近微分方程的精确解y(t)的，所以在讨论微分方程初值问题的数值解时通常先讨论欧拉法。

若把式（2.11）在(ti,ti+1)区间上积分，可得

式（2.15）右端的积分一般是很难求出的，其几何意义为曲线f(t,y)在(ti,ti+1)区间上的面积。当(ti,ti+1)区间足够小时，可用矩形面积来近似代替，即

因此，式（2.15）可以近似为

y(ti+1)=y(ti)+hf(ti,y(ti)) （2.17）

写成递推算式为

yn+1=yn+hfn （n=0，1，2，…，N） （2.18）

该式为欧拉公式。其中，h=tn+1-tn，即步长；fn=f(tn,yn)为函数y(t)在时刻tn时的导数值。

因已知y(t0)=y0，则由式（2.18）可得欧拉法的计算过程，即

由式（2.19）可知，由前一点ti上的数值yi，可以求得后一点ti+1上的数值yi+1，这种算法称为单步法。又由于式（2.18）可以直接由微分方程式（2.11）的初值y0作为递推计算时的初值，而不需要其他信息，因此它又是一种自启动算法。

与欧拉法相类似的改进的欧拉方法，是用超前一时刻的yn-1来替代式（2.18）中的前一时刻的yn，数值积分公式为

yn+1=yn-1+2hfn （2.20）

3.梯形法（预估-校正法）

欧拉法比较简单，却给我们一些启发，即在求积分时用梯形面积来代替矩形面积，可得到

写成递推公式为

用梯形公式（2.22）来计算时，产生了一个新问题，即在计算yn+1时，需要知道fn+1，而fn+1=f(tn+1,yn+1)又依赖于yn+1本身，因此，该式是不能自启动的。要求解yn+1，通常采用递推法，即取

如果，y(0)n+1,y(1)n+1,…这个序列是收敛的，那么极限就存在，即k→∞时，这个序列趋于某一个极限值，因而可以用此极限值来作为yn+1。可以证明，如果步长h取得足够小，则上述序列必定收敛。

收敛问题解决了，还有一个问题，即在具体的计算过程中，要迭代多少次才认为已经求得准确的yn+1呢？显然，迭代次数越多，求得的解越准确，但计算工作量将增加。所以，通常只迭代一次，这样计算公式就成为

通常，这类公式称为预估-校正公式，式（2.24）为预估公式，由它预估yn+1的一个值；式（2.25）为校正公式，由它得出yn+1的校正值。欧拉法每计算一次只需对f调用一次，而预估-校正法由于加入了校正过程，计算量较欧拉法增加一倍，但提高了计算精度。

4.龙格-库塔法

对于式（2.11），设其精确解是充分光滑的，若从t0时刻向前跨出一步，到t1=t0+h，t1时刻的解为y1=y(t0+h)；在t0附近将其展开成泰勒级数，则有其中，括号后的下标0，表示括号中的函数将用t=t0,y=y0代入，以下均同。

可见欧拉法是将式（2.26）截去h2及以后各项而得到的一阶一步法，所以精度较低。如果将式（2.26）多取几项后再截断，就可以得到精度较高的高阶数值解，但直接使用泰勒展开式涉及高阶导数的计算。为了避免对高阶导数的计算，龙格-库塔（Runge-Kutta）法采用间接利用泰勒展开式的思路，即用在n个点上的函数值f的线性组合来代替f的导数，然后按泰勒展开式确定其中的系数，以提高算法的阶数。

对于式（2.26）只保留到h2项，并假设它可以写成如下形式，即

y1=y0+h(a1K1+a2K2) （2.27）

其中，K1=f(t0,y0), K2=f(t0+b1h,y0+b2K1h)。

对K2式的右端函数在t=t0, y=y0处展开成泰勒级数，保留到h项，可得

将K1，K2式代入式（2.27），则有

将式（2.29）与式（2.26）相比较，可得

可见，有四个未知参数a1，a2，b1和b2，但只有三个方程，因此有无穷多个解。若限定a1=a2，则可得其中的一个解，即

a1=a2=, b1=b2=1

将它们代入式（2.27），可得一组计算公式为

将式（2.31）写成一般的递推形式，为

由于式（2.32）只保留到泰勒展开式中的h2项，而将h3以后的高阶项全部省略去了，故称其为2阶龙格-库塔法。

根据上述原理，若在泰勒展开式中保留到h4项，可得各种4阶龙格-库塔法公式。其中，典型公式为

对式（2.32）和式（2.33）所示的2阶和4阶龙格-库塔法公式可做如下的直观解释：由于y(t)满足式（2.32），已知在(tn,yn)处y(t)的导数为K1=f(tn,yn)，另外，根据线性外推的原理，即假定在(tn,tn+h)这个区间中y(t)的导数不变，则可得y(t)在t=tn+h的估计值yn+hk1。令K2为这个估计值处的导数，因此式（2.32）所示的2阶龙格-库塔法公式就是按(K1+K2)/2这个导数线性外推而得的yn+1。式（2.33）所示的4阶龙格-库塔法公式则是4点导数的加权平均，第2点和第3点是以h/2为步长，故精度较高。

根据上述解释，各种龙格-库塔法可写成如下的一般形式，即

其中，各系数应满足以下关系，即

由于龙格-库塔法有很多优点，故在许多仿真程序包中，它是最基本的算法之一。龙格-库塔法的主要特点有以下几个方面。

（1）若在展开泰勒级数时，截去h2及以后各项可得欧拉公式，因此，可以将欧拉公式看做是1阶龙格-库塔法公式，它是精度最低的一种数值积分公式。而对于1～4阶龙格-库塔法，每积分一步所需右端函数的计算次数等于阶数，而对于大于4阶的方法，右端函数的计算次数要大于阶数，使积分工作量增加，所以通常使用4阶或4阶以下的方法。

（2）龙格-库塔法属于外推单步法，可以自启动。步长h在整个计算过程中并不要求固定不变，可以根据精度要求进行改变。但是在同一步中，为计算若干个系数，则必须用同一个步长h。不论几阶龙格-库塔法，计算公式总是由两部分组成的，第一部分为上一步的结果yn，第二部分是在tn～tn+1中对f(t,y)的积分，它是步长h乘上各点斜率的加权平均。

（3）龙格-库塔法的精度取决于步长的大小和解决的方法。许多计算实例表明：为达到相同的精度，4阶方法的步长可以比2阶方法的步长大10倍，而4阶方法的每步计算量仅比2阶方法大一倍，所以总的计算量仍比2阶方法小。因此，一般系统进行数字仿真时常用4阶龙格-库塔法公式。

（4）在式（2.34）所表示的龙格-库塔法公式中，并未给出截断阶次的限制。一般来讲，为了减少计算量，总是希望尽量减少每步中计算右端函数f的次数（即减少计算K的次数）。可以证明，1阶龙格-库塔法公式至少要计算1次右端函数，即式（2.34）中的s（称为级数）至少等于1，而2阶公式的smin=2，3阶公式的smin=3，4阶公式的smin=4，5阶公式的smin=6，6阶公式的smin=7。

5.ADAMS线性多步法

以上所述的欧拉法和龙格-库塔法均为单步法。在计算中只要知道yn，f(tn,yn)的值即可递推算出yn+1。也就是说，根据初始条件可以递推计算出相继各时刻的y值，所以这种方法都可以自启动。与单步法相对应，还有一类数值积分方法，它的递推公式为

由于计算yn+1时，不仅要求知道yn及fn=f(tn,yn)，而且还要知道y及f(t,y)在tn, tn-1, tn-2,…各时刻的值，因此称为多步法，显然该计算公式不能自启动，并且在计算过程中占用的内存较大，但可以提高计算精度和速度。

由式（2.36）可知，若b-1=0，则计算yn+1公式中的右端全部已知，因此计算公式为显式，称为显式多步法。而若b-1≠0，则右端函数中含有fn+1=f(tn+1,yn+1)，因此计算公式为隐式，称为隐式多步法。

（1）ADAMS显式多步法。

ADAMS显式多步法是利用一个插值多项式来近似代替f(t,y(t))。如果在tn点以前的k个节点上，用多项式pk,n(t)近似表示f(t,y(t))，则k称为多项式的阶数。根据牛顿后插公式

其中，

并设t-tn=sh，用Pk,n(t)近似代替式（2.15）中的f(t,y)，经过简单地推导，可得ADAMS多步法的计算公式为

其中，

在式（2.39）中，当k=1时，可得欧拉公式为

yn+1=yn+hfn （2.41）

当k=2时，得到2阶亚当斯公式为

当k=3时，可得3阶亚当斯公式为

由式（2.42）和式（2.43）可以看出，如果在tn点已知yn，fn，fn-1，fn-2，那么以后求得的yn+1的值是k-1步以前各导数的线性组合，各导数都以显式出现在式（2.42）或式（2.43）中，所以称为显式线性多步法。

（2）ADAMS隐式多步法。

根据插值理论可以得出，插值节点的选择对精度有直接的影响。同样阶数的内插公式比外插公式更为精确。牛顿前插公式为

其中，∇fn为向前插分算子，定义为

用牛顿前插公式近似代替牛顿后插公式，可以得到ADAMS隐式多步法的计算公式为

其中，系数βi的值参见表2.1。

表2.1 隐式多步法系数表

如果将ADAMS方法的显式公式与隐式公式联合使用，前者提供预测值，后者将预测值加以校正，可得到预测-校正法。常用的4阶ADAMS预测-校正法的计算公式为

式（2.47）为预测公式，式（2.48）为校正公式。计算步骤：

（1）利用单步法计算预测公式（2.47）中的附加值fn-3, fn-2, fn-1,fn；

（2）计算预测值ypn+1；

（3）计算fpn+1=f(tn+1,ypn+1)；

（4）计算校正值ycn+1。

2.2.2 误差分析、稳定性分析与步长控制

2.2.2.1 误差分析

由于数值解法是以一系列离散时刻近似值y1, y2, …, yn来逼近微分方程的解y(t)，理论上讲，当n→∞时，数值解yn收敛于精确解y(t)。但是，积分步长不可能无限小，数值解与精确值之间必然存在误差。

数值计算所产生的误差有两种：截断误差和舍入误差。截断误差与采用的计算方法有关，而舍入误差则由计算机的字长所决定。下面以欧拉法为例进行介绍。

1.截断误差

对于微分方程，在不考虑舍入误差情况下，其误差称为截断误差，也称局部截断误差，它是因积分算法的固有局限性造成的误差，将其记为e(t0)，根据定义，可得

将yc在(t0+h)点进行泰勒级数展开，得

取式（2.50）右端函数的前两项，而将后面各项截断，即得到欧拉公式。余下各项可得欧拉公式的局部截断误差为

由式（2.51）可知，e(t0)与h2成正比，Ο(h2)是h→∞时与h2同阶的无穷小量，表明若步距缩短一半，局部截断误差便只有原误差值的1/4。

另外，以t=0开始继续到t=tn所积累的误差称为整体误差。一般情况，整体误差比局部误差大，其值不易估计。局部截断误差的精确值也是难以找到的，但可根据不同步距得到的不同结果来进行估计，当所取步距对于计算结果的影响达到足够小时，便认为其步距已满足要求。

2.舍入误差

舍入误差是由于计算机进行计算时，数字的位数有限所引起的误差。随着计算终止时间和积分法阶次的增加，舍入误差会增大，而且舍入误差还会因积分步距的减小而更趋严重。这是因为舍入误差在计算机进行每一步运算时都会产生，步距越小，运算次数越多，舍入误差自然会越大。舍入误差的积累值是难以精确预测的，一般认为它与h-1成正比。

最后得到的欧拉法总误差可表示为

由式（2.52）可以看出，步长h增加，截断误差Ο1(h2)增加，而舍入误差Ο2(h-1)减小。反之，截断误差Ο1(h2)减小，而舍入误差Ο2(h-1)增加。

2.2.2.2 稳定性分析

在利用数值积分法进行仿真时，常常会发生这样的现象：本来是稳定的系统，但仿真结果却得出不稳定的结论。这种现象通常是由于计算步长选得太大而造成的，即当步长选得过大时，数值积分方法可能会使各种误差恶性发展，以致引起计算不稳定。

研究数值计算方法要涉及算法的收敛性和稳定性。收敛性反映了差分方程的截断误差对计算结果的影响，它不考虑由于具体的计算机字长的限制而带来的舍入误差；稳定性反映了初值误差或某一计算步骤引入的舍入误差对计算结果的影响，它与积分步长密切相关。只有既收敛又稳定的差分方程才有实用价值。对于常用的基本数值积分方法，它们的收敛性已得到了证明，因此下面主要讨论积分方法的稳定性。

1.测试方程

如果差分方程的精确解为yn，计算所得的近似解为，则积累误差rn=yn-，对于允许的误差限ε>0，当|rn|<ε时，认为≈yn。误差rn的变化规律直接影响差分方程的稳定性，而使初值误差（或称为扰动）逐步减小的差分方程是稳定的。

给y(t)一个小的扰动δy，则

(y+δy)′=f(t,y+δy) （2.53）

将式（2.53）展开，得

由于δy很小，可略去高阶无穷小量，得

(δy)′=δy·f′(t,y) （2.55）

式（2.55）称为方程y′=f(t,y)的第一变分方程。

若f(t,y)是y的线性函数，f′(t,y)=A(t)，于是

(δy)′=A(t)δy （2.56）

若f(t,y)不是y的线性函数，由于δy很小，仍可认为式（2.55）成立。

为了方便，仅讨论∂f/∂y=A的情况（取足够小的区域，总可以如此），即

(δy)′=A·δy （2.57）

它的解为

可见，当A<0时，方程稳定；当A>0时，方程不稳定。

在A<0的情况下，用式（2.57）来讨论各种差分方程的稳定性问题。以y代替δy，可得

称式（2.58）为测试方程。

2.数值积分方法的稳定域

在实际的数值计算中，给出的初值y(t0)=y0不一定很准确（即有初始误差）；同时，由于计算机字长的有限，在计算中会有舍入误差；另外，对应于一定的步长h，还有截断误差。所有这些误差都会在逐步的计算中传播下去，对以后的计算结果产生影响。对于一种数值积分法，如果计算结果对初值误差或计算误差不敏感，就可以说该计算方法是稳定的，否则是不稳定的。对于不稳定的算法，误差会恶性发展，以致计算失败。

为了考察欧拉法的稳定性，研究测试方程=Ay，A<0（即方程是稳定的）。对此方程用欧拉法进行计算，得递推计算公式为

yn+1=yn+Ahyn=(1+Ah)yn （2.59）

假设yn（n=0, 1, 2, …）为该方程的一个解，另外，设yn+εn是一个受扰解，即

yn+1+εn+1=(yn+εn)+Ah(yn+εn)=(1+Ah)(yn+εn) （2.60）

将式（2.60）减去式（2.59），得

εn+1=εn+Ahεn=(1+Ah)εn （2.61）

显然，为了使扰动εn不随n的增加而增加，必须要求|1+Ah|<1，即该方程的特征根在单位圆内，|1+Ah|<1所对应的区域称为稳定域，即欧拉法的稳定区域为-2<Ah<0。如果所选步长不满足-2<Ah<0的条件，尽管原系统微分方程是稳定的，利用差分方程式（2.59）进行数值计算时会产生很大的误差，从而得到不稳定的数值解。这种对积分步长有限制的数值积分方法称为条件稳定积分法。

对于其他数值积分法，可以用同样的方法求它们的稳定域。一般方法：设系统方程为式（2.58）所示的测试方程，而数值积分公式为

yn+1=yn+p(Ah)·yn （2.62）

则只有当|p(Ah)|<1时，算法才稳定，由此可得出数值积分方法的稳定域。

表2.2给出了各阶龙格-库塔公式的稳定条件。

表2.2 龙格-库塔公式的稳定条件

对各种基本数值积分方法进行稳定性分析可知，除隐式1阶、2阶ADAMS线性多步法为恒稳法外，其他方法都是条件稳定的，这些条件稳定的方法中步长h的选择应限制在系统最小时间常数的数量级。对龙格-库塔法，阶数增大，稳定域略有增大；而对于ADAMS法来讲，阶数大，稳定域反而缩小。

2.2.2.3 步长控制

一个实用的仿真程序必须将步长的自动控制作为必要的手段，因为要求用户给出合适的仿真步长往往是困难的，更何况为保证仿真过程满足一定的精度而又要求计算量尽可能小，仿真步长是需要不断改变的。

一般来讲，积分步长的选择应遵循以下两个原则。

（1）首先要求保证计算的稳定性。例如，当采用4阶龙格-库塔法时，就要求步长小于系统最小时间常数的2.78倍。

（2）要求有一定的计算精度。采用数值积分方法时，有两种计算误差：截断误差和舍入误差，通常后者比较小，可忽略不计。因此，主要误差是截断误差，它将随步长h的增加而增加，使截断误差超过允许值。

以上两条都是要求限制步长h不能过大，但并不是说步长越小越好，因为步长越小，计算量越大，所以最好的选择是在保证计算稳定性及计算精度的前提下，选择最大步长。

对步长的控制可以通过对截断误差的估计来自动改变和控制步长，即每积分一步，都设法估计出计算误差，若误差在允许范围之内，则该步计算结果有效，并设法调整下一步的步长（一般是略微放大）；若估计出的误差大于允许误差，则该步计算结果无效，设法减小步长，重新进行该步的积分。

2.2.3 积分方法的选择

连续系统一般用微分方程（组）来描述，这些微分方程（组）的数值解包含某些积分过程。在进行连续系统仿真时，积分方法的选择十分重要，但究竟如何进行选择，没有一个具体的办法，下面列出几点需要考虑的因素。

1.精度问题

数值积分的精度主要由三类误差影响：截断误差、舍入误差和积累误差。其中截断误差与积分公式的阶次有关，阶数越高，截断误差越小，一般减小步长可减小每一步的截断误差。舍入误差则与计算机字长有关，字长越长，舍入误差越小。积累误差是上述两类误差累积的结果，它与积分时间长短有关，一般积分步长越小，则积累误差越大（在一定积分时间下），所以在一定积分方法条件下，若从总误差考虑必须有一个最佳步长值。

2.计算速度问题

计算速度决定于计算的步数及每一步积分所需时间，而每步的计算时间与积分方法有关，它主要决定于计算导数的次数。例如，由于4阶龙格-库塔法每一步中要计算4次导数，所以每步花费时间较多。为了加快计算速度，在积分方法已定的条件下，应在保证一定精度下，尽量以选用大步长，以缩短积分时间。

3.数值解稳定性问题

数值积分方法实际上是将微分方程化成差分方程进行求解，对于一个稳定的微分方程组，经过变换得到的差分方程不一定是稳定的，不同的积分方法具有不同的稳定域，所以在选择积分方法时，要考虑选择具有较好稳定性的积分公式。

4.自启动问题

在应用积分法求解时，若可以直接从微分方程的初值开始，就属于自启动，如单步法中的欧拉法和龙格-库塔法。而ADAMS线性多步法除了初始值以外，还需要以前时刻的值，所以无法自启动，而需要用另一种自启动后，再应用多步法积分。

2.2.4 刚性系统的仿真算法

在实践中，在研究液压系统、电子网络、控制系统时，常常会碰见这样的情形：一个高阶系统中常有不同的时间常数相互作用。以惯性导航为例，修正回路时间常数大，稳定回路时间常数小。导弹等航行器的运动也是如此，质点运动较慢，偏航与俯仰运动较快。所有这些现象，主要是由于系统模型中的一些小参数，如小的时间常数、小质量等的存在而引起的。描述这些现象的微分方程，在数学上常常称为刚性方程，这种系统就称为刚性系统。

1.刚性系统及其特点

一个刚性系统可以这样描述，对于n阶微分方程组

其中，y=[y1 y2 … yn]T, f=[f1 f2 … fn]T ,η=[η1 η2 … ηn]T,为n维向量。该方程组雅可比矩阵∂f/∂y的特征值为λi（i=1, 2, …, n），特征值的实部表示系统衰减的速率，其最大特征值与最小特征值实部之比，即

作为系统刚性程度的度量。当ρ≫1时，系统为刚性系统，或称为Stiff系统。

现以一个3阶微分方程组为例，讨论刚性系统及它在数值求解上的特点，即

其中，y=[y1 y2 y3]T；A为3×3维矩阵，设

满足初始条件y0=[2 1 2]T的解为

其中，λ1=-120, λ2=-50, λ3=-0.1为A的3个特征值。

于是，方程的解中有对应于λ2或λ3的一个衰减很快的分量，它们的权重很快变得很小。然而，从数值解稳定性的要求出发，步长又受最小时间常数的限制。例如，对于欧拉法，应保证|λih|<2，即|λ1h|<2，因此最大步长只能为1/60。虽然对于对应λ1的解分量，很快就没有实际价值了，但是在整个积分区间，由于受到绝对稳定的限制，步长又必须取得很小。因此，当一个系统的矩阵的特征值范围变化很大时，数值求解会引起很大的困难。

对刚性系统进行数字仿真，其最大的困难：要求步长限制在最小时间常数（相当于最大特征值实部的倒数）的数量级，而系统的过渡过程的时间却又决定于最大时间常数（相当于最小特征值实部的倒数）。因而占用大量的计算时间，有时会由于计算的舍入误差而导致计算失败。

刚性系统在实践中的普遍性和重要性已得到广泛的重视，这种方程的数值解已成为常微分方程数值求解研究中的重点。

微分方程数值解中对步长的限制并不是物理本质上的原因，它仅是为保证数值解稳定性而采取的措施。而刚性方程求解数值解时，要解决稳定性和计算次数的矛盾，就应放宽对步长大小的限制。到目前为止，已提出的解刚性方程的数值方法，基本上分为显式公式、隐式公式和预测-校正型。

（1）显式公式常用雷纳尔法，其着眼点：在保证稳定的前提下，尽可能地扩大稳定区域。这一方法的优点是，它是显式的，便于程序设计。对一般条件的方程，它就还原为4阶龙格-库塔方法，而对刚性方程它又有增加稳定性的好处。

（2）隐式1、2阶公式是恒稳的，因而适合于解描述刚性系统的方程组，但这种方法每计算一步都需要进行迭代，计算工作量大，在工程上使用有一定的困难。因而，在应用中常使用Rosenbrock提出的半隐式龙格-库塔法。

（3）预测-校正型中常用的解刚性方程的方法是吉尔（Gear）法。吉尔法首先引进了刚性系统稳定性的概念，它可以满足稳定性，而降低对步长的限制；它又是一种通用的方法，不仅适用于解刚性方程组，而且也适用于解非刚性方程组。

2.Stiff稳定域

一般的数值积分方法，如龙格-库塔法、隐式或显式ADAMS法，除隐式ADAMA法中的1阶、2阶法为恒稳法外，其他都是条件稳定的。对于刚性系统，若采用条件稳定法，由于要求步长限制在最小时间常数的数量级，而系统的过渡过程的时间却决定于最大时间常数，因此计算量极大。如果用恒稳方法，虽然从计算稳定性方面考虑可以采用大步长，但由于只有1阶和2阶方法，截断误差较大，因此也不适用于刚性系统的仿真。

实际上，为仿真刚性系统不一定要求用恒稳法，而只要求算法要具有如图2.1所示的稳定域，这种稳定域称为Stiff稳定域。

定义：一种数值方法，如果它是收敛的，且存在正常数α和θ，使该方法在区域R1和R2中是稳定的，则称为是Stiff稳定的。其中

3.求解刚性方程的吉尔法

吉尔法首先推出了满足Stiff稳定域的多步法。试求解常微分方程初值问题，其递推方程的形式为

式中有2p+3个系数（αi,bi）待定，但只有n+1个方程，为此，设定p=K-1，且令b0=b1=…=bk+1=0，则式（2.69）可改写为

yn+1=a0yn+a1yn-1+…+ak-1yn-k+1+hb-1fn+1 （2.70）

而约束条件改写为

由此可得到各阶隐式吉尔法的系数，参见表2.3。

表2.3 各阶隐式吉尔法的系数

隐式吉尔法的稳定域符合Stiff稳定域的条件，因此隐式吉尔法适用于对刚性系统的仿真。同样也可以推出显式吉尔法，但它们的稳定域不满足Stiff稳定域的特性，在此不再赘述。

2.2.5 实时系统的仿真算法

1.实时仿真的要求

假设在计算机上仿真的连续系统由非线性常微分方程描述为

其中，计算机的输入序列是u，由实物经计算机的输入接口输入给计算机，u=u(kh), k=0, 1, 2,…。

计算机从时刻kh开始，根据用户所采用的不同仿真算法，利用y(kh), u(kh)和kh时刻以前的数值，计算出yk+1=y(kh+h)。很明显，实时仿真算法的第一个要求是计算机求解式（2.72）一步解的实际时间必须等于或小于h秒，以便与实物的运行时间同步。

对实时仿真算法的第二个要求是，计算机在kh时刻求解方程时，不能要求从实物上取得kh时刻的值，即计算机的输入也必须满足实时仿真的条件。例如，对于龙格-库塔法，前面所介绍的计算公式是不适用于实时仿真的了。

如果对式（2.72）采用2阶龙格-库塔法进行求解，其递推公式可写为

由于算法在一个仿真步中要计算两次右端函数，所以可假定在h/2时间内计算机计算一次右端函数。即：要计算在tn+1=tn+h时刻的输出yn+1，计算机在tn～(tn+h)/2时间内根据yn和tn时刻的外部输入值un来计算k1，而由于只有在tn+1=tn+h时计算机才能真正获得外部输入值un+1，才具备计算k2的条件，所以到tn+1+h/2时刻才能计算出yn+1，并输出到外部设备，也就是说，计算机输出就要滞后半个计算步距。可见，式（2.73）不适合于实时仿真。其他公式的分析同上。

2.实时仿真的要求实时仿真算法

为了适用于实时仿真计算，一般采用下面的方法计算。

（1）选择ADAMS多步法。因为在这类算法中，为计算yn+1，只要求知道tn及tn以前的各类右端函数值，tn以前的各类右端函数值可事先存储在内存中，而tn时刻的各类右端函数值和外部输入均可在tn+h这段时间内计算出来，或由外部设备输入给计算机，所以yn+1不会被延迟。如果用隐式算法，则可用显式法计算预估值。

（2）合理地选择龙格-库塔法计算公式中的系数，使之能适用于实时仿真。例如，可将2阶龙格-库塔法公式改写为适用于实时仿真的计算公式

使得在(tn+h)/2～tn+1时间内，根据(tn+h)/2时刻的外部输入值就可计算k2，而无须获得tn+1时刻的外部输入值un+1，从而可避免输出的延迟。