连续系统的计算机模拟

第
2章 连续系统的计算机模拟
本章讨论连续系统的模拟技术,由于这类系统中状态随时间连续动态地变化，常常具有一定的规律，故可用一些数学方程来描述，这些方程就是系统的数学模型，通常以微分方程、代数方程为多见。

下面将介绍利用数值积分法对连续系统进行数字模拟的基本原理和具体方法，并给出数值积分法中几个常用的算法以及实现这些算法的计算程序，最后介绍两个建模实例。

数值积分法不仅方法种类多，而且有较强的理论性，本章由浅入深地介绍几种常见的数值解法。

主要为单步法中的四阶龙格--库塔法与默森法和多步法中的亚当斯法。

使用数字计算机对连续系统进行模拟，首先必须将连续系统离散化，并将它转化为差分方程，以建立所谓的模拟系统的数学模型。

描述连续系统动态特征的数学模型是多种多样的，除微分方程外，还有传递函数、结构图及状态方程等，由于篇幅所限，本书不讨论后两种方法。

建立系统的模拟模型之后，就要选择计算机语言（也叫算法语言）编写系统模拟程序，在计算机上运行，将结果保留在数据文件中以待传输和处理。

由于模拟的目的不同，可以选用不同的模拟模型和算法，其特点是运算精度高，对于不同的计算机，字长一般在16位---72位之间，也可采用浮点运算和双精度运算，其精度一般可达千万分之一到百万分之一。

当然结果的精度与所选的算法有关。

这可以根据实际需要选择机器、算法和模拟的步长。

数字计算机储存容量大，可进行各种运算，在以前认为是不可能解决的问题，利用数字计算机都可容易地或有可能得到解决。

本章介绍的方法适应性较强，应用也十分广泛。

数字计算机上还有各种功能的软件包（即一些子程序），用户可以稍加修改或不经修改就可以用于自己的模拟程序中，解决自己的实际问题，使用非常方便。

2．1 欧拉（Euler ）法
在讨论连续系统的计算机模拟之前，让我们先看一个化学反应的例子，通过这个例子我们可以看到怎样使用数字计算机模拟一个实际问题，虽然介绍的是欧拉法，但是分析问题的思路同样适用与其它数值积分法。

当两种物质A 和B 放到一起产生化学反应时，产生第三种物质C ，一般一克A 与一克B 结合产生2克的C 物质，形成C 的速率与A 和B 的数量乘积成正比，同样C 也可分解为A 和B ，C 的分解速率正比与C 的数量，即在任何时刻，如果a,b,和c 是化学物质A,B,和C 的数量，即在任何时刻，如果a,b,和c 是化学物质A,B,C 的数量，它们的增加和减少的速度服从下列微分方程。

C K ab K dt dc ab K C K dt
db ab K C K dt
da 21121222 -=-=-= （2.1）
其中K 1和K 2 是比例常数（一般而言这些比例常数会随温度和压力发生变化，但在模拟过程中，为了简化模型，一般不允许其变化，故一律视为常数）。

在给出常数K 1 和K 2 值以及A 和B 的数量（C=0）后，我们希望能确定有多少C 物质产在出来，这种化学反应速率
的决定在化学工业上是有意义的。

模拟该系统的一个直接的方法是在t ＝0时开始，使t 以Δt 间隔增加。

假定化学量在Δt 时间步长内不变，而只能在Δt 结束的瞬间发生变化，这样在每个Δt 结束时的A （或B 或C ）的数量就可以从Δt 开始时的数值由下式求出
(2．2) 同样的方程b(t+Δt)和C(t+Δt),也可写出。

假定模拟周期为T,可将T 分成N 个小的时间步长Δt,及
T ＝N ·Δt
在时间为零时，我们知道a(0)、b(0)和c(0)的初始值,从这些初始值及常数K 1和K 2值出发，就可以计算出Δt 时间内的化学量的变化值。

a(Δt)=a(0)+[K 2•C(0)-K 1a(0)·b(0)]·Δt
b(Δt)=b(0)+[K 2·C(0)-K 1a(0)·b(0)]·Δt
c(Δt)=c(0)+[2K 1·a(0)·b(0)-2K 2 c(0)]·Δt
使用这些值，又可计算系统的下一个状态，即2Δt 时的状态。

a(2Δt)=a(Δt)+[K 2•c (Δt)-K 1a(Δt)·b(Δt)]·Δt
b(2Δt)=b(Δt)+[K 2·c(Δt)-K 1a(Δt )·b(Δt )]·Δt
c(2Δt)=c(Δt )+[2K 1·a(Δt )·b(Δt )-2K 2 c(Δt)]·Δt
使用2Δt 时的系统状态，又可写出3Δt 时的状态，依此类推，以Δt 为间隔，进行N 步计算，就可写出周期T 的系统状态得到所期望的结果，这个过程可用图2.1表示。

图2.1 化学反应例子模拟程序框图
模拟程序使用C 语言编写，程序中的初值如下：
K 1=0.008/g.min K 2＝0.002/min, a(0)=100g, b(0)=50g C(0)=0, T=5mins, ⊿t=0.1min, N=50 )()()(t dt
t da t a t t a ∆+=∆+
其程序清单如下：
#include <stdio.h>
#include <string.h>
float k1,k2;
static float A[53],B[53],C[53],delta,t;
void strut(int);
main()
{
int i;
A[1]=100.0;
B[1]=50.0;
C[1]=0.0;
t=0;
delta=0.1;
k1=0.008;
k2=0.002;
for(i=1;i<53;i++)
{
printf("%2d",i);
printf("%10.2f,%10.2f,%10.2f,%10.2f\n",t,A[i],B[i],C[i]); strut(i);
}
return;
}
void strut(int i)
{
A[i+1]=A[i]+(k2*C[i]-k1*A[i]*B[i])*delta;
B[i+1]=B[i]+(k2*C[i]-k1*A[i]*B[i])*delta;
C[i+1]=C[i]+2.0*(k1*A[i]*B[i]-k2*C[i])*delta;
t=t+delta;
}
运行此程序，输出模拟结果如下：
1 0.00, 100.00, 50.00, 0.00
2 0.10, 96.00, 46.00, 8.00
3 0.20, 92.47, 42.47, 15.06
4 0.30, 89.33, 39.33, 21.34
5 0.40, 86.52, 36.52, 26.95
6 0.50, 84.00, 34.00, 32.00
7 0.60, 81.72, 31.72, 36.55
8 0.70, 79.66, 29.66, 40.69
9 0.80, 77.77, 27.77, 44.45
10 0.90, 76.05, 26.05, 47.89
11 1.00, 74.48, 24.48, 51.04
12 1.10, 73.03, 23.03, 53.94
13 1.20, 71.70, 21.70, 56.61
14 1.30, 70.46, 20.46, 59.07
15 1.40, 69.32, 19.32, 61.36
16 1.50, 68.26, 18.26, 63.48
17 1.60, 67.28, 17.28, 65.44
18 1.70, 66.36, 16.36, 67.28
19 1.80, 65.51, 15.51, 68.99
20 1.90, 64.71, 14.71, 70.59
21 2.00, 63.96, 13.96, 72.08
22 2.10, 63.26, 13.26, 73.48
23 2.20, 62.60, 12.60, 74.79
24 2.30, 61.99, 11.99, 76.03
25 2.40, 61.41, 11.41, 77.18
26 2.50, 60.86, 10.86, 78.27
27 2.60, 60.35, 10.35, 79.30
28 2.70, 59.87, 9.87, 80.27
29 2.80, 59.41, 9.41, 81.18
30 2.90, 58.98, 8.98, 82.04
31 3.00, 58.57, 8.57, 82.86
32 3.10, 58.19, 8.19, 83.63
33 3.20, 57.82, 7.82, 84.36
34 3.30, 57.48, 7.48, 85.05
35 3.40, 57.15, 7.15, 85.70
36 3.50, 56.84, 6.84, 86.32
37 3.60, 56.55, 6.55, 86.91
38 3.70, 56.27, 6.27, 87.46
39 3.80, 56.00, 6.00, 87.99
40 3.90, 55.75, 5.75, 88.50
41 4.00, 55.51, 5.51, 88.97
42 4.10, 55.29, 5.29, 89.43
43 4.20, 55.07, 5.07, 89.86
44 4.30, 54.86, 4.86, 90.27
45 4.40, 54.67, 4.67, 90.66
46 4.50, 54.48, 4.48, 91.03
47 4.60, 54.31, 4.31, 91.39
48 4.70, 54.14, 4.14, 91.73
49 4.80, 53.98, 3.98, 92.05
50 4.90, 53.82, 3.82, 92.35
51 5.00, 53.68, 3.68, 92.65
52 5.10, 53.54, 3.54, 92.93
以上例子的算法就是数值积分法中最简单的一种方法叫作欧拉法, 从这个例子中我们
可以看出使用计算机解题的过程，由于欧拉法本身的特点，决定了其计算精度差, 现在几乎无人在实际工作中使用，但它导出简单， 能说明构造数值解法一般计算公式的基本思想, 模拟程序也清晰易懂，故人们乐意用它做为构造数值解法的入门例子。

其一般解法如下：
设给定微分方程
(2.3)
在区间(t n ,t n+1)上求积分,得
y(t n+1)=y(t n )+∫f(t,y)dt (2.4)
如积分间隔足够小，使得t n 与t n+1之间的f (t,y)可近似的看成常数f (t n ,y n ), 就可以用矩形面积近似地代替在该区间上的曲线积分，于是在t n+1时的积分值为
1),()1( +=+ +n n n n n y h y t f y t y (2.5)
将上式写成以下差分方程形式:
1,2,3,n 1=⋅+=+n n n f h y y （2.6） 这就是欧拉公式。

它是一个递推的差分方程，任何一个新的数值解y n+1均是基于前一个数值解以及它的导数f(t n ,y n )值求得的。

只要给定初始条件y 0及步长h 就可根据f(t 0,y 0)算出y 1的值，再以y 1算出y 2，如此逐步算出y 3, y 4,…，直到满足所需计算的范围才停止计算。

欧拉法的基本思路是把连续的时间t 分割成等间隔的Δt, 在这些离散的时刻算得函数值，根据这些值在函数图上可得到一条折线，所以欧拉法又叫折线法，其特点是分析方法简单，计算量小，但计算精度低（后面将讨论欧拉法与其它方法的比较）。

下图为欧拉折线法的几何意义。

n-1 n n+1 n+2
图2.2 欧拉法的几何意义
如果用梯形面积来代替一个小区间的曲线积分，就可克服以小矩形计算的缺点，提高精度，梯形法计算公式为
⎩⎨⎧==0
)0(),(y y y t f y
- )(2y )],(),([2
1n 111++++++=++=n n n n n n n n f f h y t f y t f h y y (2.7)
上式为隐式公式，因为公式右端含有y n+1, 这是未知的待求量，故梯形法不能自行启动运算，而要依赖于其它算法的帮助，比如说用欧拉公式法求出初值，算出)(1+n y y 的近似值F
n y 1+，
然后将其带入微分方程，计算1+n f 的近似值),(111P n n F n y t f f +++=，最后利用梯形公式反复迭代。

如迭代一次后就认为求得了近似解，这就是改进的欧拉法，其公式为
预估 ),(1n n n P n y t hf y y +=+ （2-8）
校正 []
),(),(2111P n n P n n n C n y t f y t f h y y +++++= （2-9） 上式中第一个为预估公式，第二个为校正公式。

通常这种方法称为预估矫正法。

在校正公式中计算了两点的斜率，再求其平均值，故计算量比欧拉法要大些。

2．3 数值积分的几个基本概念
1. 单步法与多步法
数值积分法都用递推公式求解, 而递推公式是步进式的, 当由前一时刻的数值y n 就能计算出后一时刻的数值y n+1时, 这种方法称为单步法, 它是一种能自启动的算法, 欧拉法是单步法. 反之, 如果求y n+1时要用到 y n , y n-1, y n-2 ,…等多个值, 则称为多步法,由于多步法在一开始时要用其它方法计算该时刻前面的函数值, 所以它不能自启动,以上所讲的改进的欧拉法就是一个多步法的例子。

2 显式与隐式
在递推公式中, 计算y n+1时所用到的数据均已算出的计算公式称为显示公式. 相反,在算式中隐含有末知量y n+1 的则称为隐含公式. 这需用另一个显示公式估计一个值, 再用隐式公式进行迭代运算, 这就是预估----校正法, 这种方法不能自启动.
用单步法与显示公式在计算上比用多步法和隐式公式方便. 但有时要考虑精度与稳定性时, 则必须采用隐式公式运算.
3 截断误差
一般使用台劳级数来分析数值积分公式的精度. 为简化分析, 假定前一步得的结果y n 是准确的, 则用台劳级数求得t n+1处的精确解为
)()(!1)(!21)()()(1)(2+++++=++r n r r n n n n h O t y h r t y h t hy t y h t y (2.10)
与前面的递推公式比较, 欧拉法是从以上精确解中只取前两 项之和来近似计算y n+1的 , 由这种方法单独一步引进的附加误差通常称着局部截断误差，它是该方法给出的值与微分方程解之间的差值，故又称为局部离散误差。

欧拉法的局部截断误差为
)()()(21h O t y t y n n =-+
不同的数值解法，其局部截断误差也不同。

一般若截断误差为)(2
h O ，则方法称为r 阶的，
所以方法的阶数可以作为衡量方法精确度的一个重要标志。

欧拉法是一阶精度的数值解法，而改进的欧拉法(梯形法)相当于取台劳级数的前三项，故其解为二阶精度。

4 舍入误差
由于计算机的字长有限, 数字不能表示得完全精确, 在计算过程中不可避免地会产生舍入误差. 舍入误差与步长 h 成反比, 如计算步长小, 计算次数多则舍入误差大. 产生舍入误差的因素较多, 除了与计算机字长有关外, 还与计算机所使用的数字系统, 数的运算次序及计算f(t,y)所用的程序的精确度等因素有关.
5 数值解的稳定性
以上对连续系统的模拟用到了差分方程, 把初始值代入递推公式进行迭代运算, 如果原系统是稳定的, 由数值积分法求得的解也应该是稳定的. 但是, 在计算机逐次计算时, 初始数据的误差及计算过程中的舍入误差对后面的计算结果会产生影响. 而且计算步长选择得不合理时, 有可能使模拟结果不稳定. 以下我们简要地讨论这一问题.
差分方程的解与微分方程的解类似, 均可分为特解与通解两部分. 与稳定性有关的为通解部分, 它取决于该差分方程的特征根是否满足稳定性要求。

例如, 为了考虑欧拉法的稳定性，可用检验方程y=λy 。

其中λ为方程的特征根。

对此方程有
要使该微分方程稳定，须使下式成立
|1+λh|<1
有此可得
|λh|<2或h<2T
要保证用欧拉法进行的计算是稳定的, 积分步长 h 必须小于系统时间常数的两倍. 所以积分步长的选择就要慎重, 不能随意取任何值。

2．4 龙格--贝塔(Runge--Kutta) 法
仍以(2.3)方程为例.已知y(t 0)=y 0,假设我t 0 开始以h 增长, t 1=t 0+h ，t1 时刻为 )(01h t y y +=,在0t 附近展开成台老级数，保留2h 项，则有 20001)(2),(h t
f dt dy y f h h y t f y y ∂∂+∂∂++= (2-11) 上式括号内下标为0，表示括号中的函数用t=t 0,y=y 0代入，以下均同。

假设上式的解可写成如下形式
其中, K 1=f(t 0,y 0)
对K 2式右端得函数在t=t 0, y=y 0处展开成台劳级数，保留h 项,可得到:
将K 1 与K 2代入 ( 2-12)式得:
n
n n n y h h y y y ) 1( 1λλ+=+=+)
,(110202h K a y h C t f K ++=h y
f K a t f C y t f K 0112002)()(∂∂+∂∂++≈
(2.12) )(221101h K b K b y y ++=
将上式与（2-12）式比较，可得系数a 1,b 1,c 2,b 2方程如下：
以上三个方程,四个未知数,因此有无穷个解,若限定 b 1=b 2 则可得一个解:
1a 21 2121====C b b
将它们代入(2-12)式可得出一组计算公式
- )(22101K K h y y ++=
其中 K1=f(t0,y0),k2=f(t0+h,y0+k1h),写成一般的递推形式如下：
式中：
),( )2
,2( )2
,2( )
,( 3423121hK y h t f K K h y h t f K K h y h t f K y t f K n n n n n n n n ++=++=++==
四阶龙格—库塔法公式精度较高，故其应用较为普遍。

下面再给出几组系数值
⎪⎪⎪⎩
⎪⎪⎪⎨⎧===+21
211122221a b c b b b h y
f K a t f C y t f h b y t hf b y y 0112
00200101)(),([),(∂∂+∂∂+++=(2.13) )(2
)(2111K K h y y t y n n n ++=≈++二阶龙格－库塔法。

两项，故此方法被称为和取只的三次方，计算过程中正比于此递推公式的截断误差上的高阶项略去，所以以
两项，而将式中只取其中，232121,),,(),,(h h h h h h h K y h t f K y t f K n n n n ++==四阶龙格库塔公式为
，其截断误差为项和级数时保留库塔公式，在展开台劳一般使用的为四阶龙格,,,, 5432h h h h h -(2.14) )22(6
)(432111K K K K h y y t y n n n ++++=≈++
21b 21b 1C 1b 4
1b 21C 1b 0b 2
1C 1
2122122121==========a 设 读者可以写出其相应的四阶龙格—库塔法公式。

龙格—库塔法的特点如下：
1.龙格－库塔法有多组a 1,c 2 b 1 b 2值，故可有多种龙格－库塔公式，常使用的有：
所有龙格－库塔公式都有以下几个特点：
（1）在计算y n ＋1时只用到y n ，而不直接用y n －1，y n －2等项，也就是计算后一步时 只用到前一步的计算结果，故称为单步法，单步法可以自启动，且使用的存储量也小。

（2）在计算过程中，步长h 可以改变，这要由计算精度来定，但是在一步中，为计算 若干个系数K i (称成为龙格－库塔系数)，必须使用同一个步长h 。

(3). 不论是几阶龙格－库塔公式，计算式中均为两部分组成。

第一部分为上一步的结果y n ，第二部分为h=t n--t n －1中对f(t,y)的积分。

它是步长h 乘上各点斜率的加权平均。

例如上面计算的四阶公式中，取四点的斜率：K 1、K 2、K 3和K 4, 然后对K 2和K 3取两份，K 1和K 4各取一份，进行加权平均。

3.龙格－库塔法的精度取决于步长h 的大小及求解的方法。

实践证明：为达到同样的精度，四阶方法的步长h 可以比二阶方法的h 大10倍，而四阶方法的每步计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小，使用的CPU 机时也少。

一般工程中四阶龙格－库塔公式已能达到要求的精度，故不再使用更高阶的公式。

下表示出了f 的计算次数与精度阶数的关系。

可见精度的阶数与计算函数值f 的次数之间的关系并非等量增加。

其精度最低。

如果将h 取得很小，能否用欧拉公式得到很高的精度呢？从理论上讲是可以的，但是可得
以上的高次项略去以及而将这一项只取若在展开成台劳级数时,h h ,h ,.422)
,(1n n n n y t hf y y +=+,h ,,2她的截断误差正比库塔公式时一阶龙格因此欧拉公式也可看成这是欧拉公式-⎪⎪⎩
⎪⎪⎨⎧++==+=+)2,2(
),(12121h K y h t f K y t f K h K y y n n n n n n
实际上由于计算机字长有限，在计算中有舍入误差。

步长h 越小，计算的步数越多，舍入误差的积累会越来越大，故用欧拉法时不能用很小的h 。

2.5 四阶龙格—库塔法模拟程序及应用
目前已有多种四阶龙格—库塔法模拟（仿真）软件包推出，虽然子程序稍有不同，但总的结构还是相同的。

对于连续系统的龙格—库塔法的计算机程序，其大致结构如下图所示：
主 函 数
输 入 及 运 行 输 出
初 始 化
微分方程 数值积分 显 示 打印 绘图
右函数
图2.3 龙格—库塔法程序简要框图
图2.3中每一个程序模块的功能为：
1. 主函数：主函数实现对整个模拟系统的控制, 这是通过调用各个子程序实现的。

2. 输入及初始化函数：主要对系统参数输入初值，例如模拟时间，积分步长，方程阶数等。

这可以从键盘输入，也可从数据文件输入。

3．运行模块: 在运行子程序中, 将反复调用数值积分和方程右函数模块,进行迭代计算，可以每计算一步便显示一次结果，也可以计算数步显示一次结果，屏幕的显示使用户可以随时监视系统运行的进程，以便人工干预。

4．输出子程序：按要求输出打印结果，可以打印，也可以绘图，视需要而定。

四阶龙格—库塔法公式前面已经给出, 现在再将它写成可编程的向量形式
),,,,()2
1,21,21,2()2
1,,21,21,2()
,,,,()2(6
1323213142,2221213222211122114321,1,N Nn n n n i i N Nn n n n i i N Nn n n n i i Nn n n n i i i i i i n i n i K y K y K y h t hf K K y K y K y h t hf K K y K y K y h t hf K y y y t hf K K K K K y y ++++=++++=++++==++++=+
式中，i =1,2,3,...N, N 为方程阶数。

为调微分子和分别换成和令)(),()(),(),(,,,)4()3()2()1(4321I D I D I D I D I D K K K K h DT i i i i = 程序求得，这些分别为变量的导数，这样上式变为：
Y(I)=YS(I)+0.166667* 2.0*(RK 1（I ）+RK 3(I))+4.0*RK 2(I)+D(I)*DT
其中RK1(I)=D(I)*0.5DT
RK2(I)=D(I)*0.5DT
RK3(I)=D(I)*DT
代入上式：
Y(I) =YS(I) +0.166667*DT*[D (I)+2.0*D (I)+2.0*D (I)+D (I)] 即y n+1 =y n +h (K1+2K2+2K3+K4)
设有一个微分方程：
y(t) +P1 y(t) +y(t) =1.0
令y1(t), y2(t)为两各个状态变量，且
y1(t) =y(t),y1(t) =y2(t)
代入原方程得：y2(t）=-y1(t)-P1y1(t)+1.0
其中P1为常数，且P1 = 0.1
初始条件：y1(0)=0, y2(0)=0
模拟参数初值：
模拟时间为50.0
积分步长为0.1
打印点数为101
程序清单如下：
主程序：
/******p23 example******/
#include "stdio.h"
#include "conio.h"
#include "math.h"
#define NORDER 3
#define NPARAM 2
float y[NORDER],d[NORDER],p[NPARAM],t;
void output(void)
{int j;
printf("%9s","time");
for(j=0;j<NORDER;j++)
printf(" y[%d]",j);
putchar('\n');
}
void difq(void)
{d[0]=-p[0]*y[0]*y[1];
d[1]=p[0]*y[0]*y[1]-p[1]*y[1];
d[2]=p[1]*y[1];
}
void run(float dt)
{int i;
float rk[3][NORDER],ys[NORDER];
difq();
for(i=0;i<NORDER;i++)
{rk[0][i]=d[i]*0.5*dt;
ys[i]=y[i];
y[i]=ys[i]+rk[0][i];
}
t+=0.5*dt;
difq();
for(i=0;i<NORDER;i++)
{rk[1][i]=d[i]*0.5*dt;
y[i]=ys[i]+rk[1][i];
}
difq();
for(i=0;i<NORDER;i++)
{rk[2][i]=d[i]*dt;
y[i]=ys[i]+rk[2][i];
}
t+=0.5*dt;
difq();
for(i=0;i<NORDER;i++)
y[i]=ys[i]+1.0/6*(2*(rk[0][i]+rk[2][i])+4*rk[1][i]+d[i]*dt); }
void main(void)
{char c;
do
{float dt,tmax,co,tnext,tol;
int j;
printf("This program will find the root of:\n");
printf(" {y[0]'(t)=-p[0]*y[0]*y[1]\n");
printf("&&{y[1]'(t)=p[0]*y[0]*y[1]-p[1]*y[1]\n");
printf("&&{y[2]'(t)=p[1]*y[1]\n");
for(j=0;j<NORDER;j++)
{printf("Now,please input the first y[%d]:",j);
scanf("%f",&y[j]);
}
for(j=0;j<NPARAM;j++)
{printf("Please input p[%d]:",j);
scanf("%f",&p[j]);
}
printf("Please input the time of simulation:");
scanf("%f",&tmax);
printf("Please input the time of one step:");
scanf("%f",&dt);
printf("Now,this is the result:\n");
co=5*dt;
tol=0.0001*co;
tnext=0;
output();
for (t=0;t<=tmax+tol;)
{int k;
if(fabs(tnext-t)<tol)
{tnext+=co;
printf("%10.4f",t);
for(k=0;k<NORDER;k++)
printf("%10d",(int)(y[k]+0.5));
putchar('\n');
}
run(dt);
}
printf("\nRun this program again[Y]:");
c=getche();
putchar('\n');
}
while(c=='\r'||c=='Y'||c=='y');
}
程序中所用的变量和数组说明如下：
Y(20) 状态变量数组
G(20) 状态变量的一阶导数
P(20) 存数系统参数，本例中仅一个参数，P(1)=P1=0.1 TMAX 模拟时间
DT 积分步长，即前面的h
NP 打印点数
NORDER 系统阶数
NPARAM 系统参数个数
OUTPUT 打印输出控制变量，.TRUE. 为打印FALSE
为不打印
INIT 打印表头控制变量
A(8) 8个参数变量
C O 打印时间间隔
NOUT 已打印点数计数器
TNEXT 打印点处的时间值
SS 剩下的打印点数
NLIST 输出结果
对以上程序编译，运行，得到如下结果
TIME Y （1） Y （2）
0.0000 0.0000 0.0000
0.5000 0.1204 0.4676
1.0000 0.4450 0.8008
1.5000 0.8863 0.9265
2.0000 1,3332 0.8247
2.5000 1.6788 0.5310
9.5000
1.6666 0.8333 10.0000
1.6666 0.8333
2.6 变步长的龙格—库塔法
以上提到，步长h 的选择十分重要，h 太大不能达到预期的精度要求，h 太小则增加了模拟时间，且有可能影响计算精度，要克服这一问题，就要采用变步长方法，使计算量尽可能的小，且精度也合乎要求。

步长的自动控制是根据每一步的误差的大小来实现的。

为了得到每一步的局部误差的估计值，可以采用两种不同阶次的递推公式（一般采用低一阶的RK 公式，同时计算y n+1并取两者之差），要使计算量最少，可以选择ＲＫ系数，要求两个公式中的K i 相同，使中间结果对两种RK 公式都适用，则这两个公式的计算结果之差就作为误差估计。

1957年Merson （默森）给出了一个四阶变步长的方法，称为龙格—库塔—默森法。

其四阶公式为
y n+1 =y n +h (K 1+4K 4+K 5)/6 (2.15)
K 1 =f (t n , y n )
K 2 =f (t n +h/3, y n +(h/3)*K 1)
K 3 =f (t n +h/3, y n +(h/6)(K 1+K 2))
K 4 =f (t n +h/3, y n +(h/8)(K 1+3K 3))
K 5 =f (t n +h, y n + (h/2)(K 1-3K 2+4k 4))
计算估计误差的三阶公式如下 )1293(6y ~ 4311n K K K h y n +-+=+
其绝对误差为
)892(6y ~E 543111n K K K K h y n -+-=-=++
这一算法是四阶精度，三阶误差，故称为RKM3—4法，目前使用较多。

其稳定域和一般RK4法相近，缺点是计算量稍大，每步计算5次f 值，除用绝对误差外，也可用相对误差，最大相对误差RE 为
)
1y max(RE 12n +-+=++n n y y E
在每一步计算后，先计算误差，根据误差的大小来控制步长，控制策略如下：
（1）当误差ERR 大于预先设定的最大允许误差E max 时，则缩小步长，一般将步长减半，并以新步长重新计算以后再比较。

（2）当误差ERR 小于预先定的最小允许误差E min 时，步长增加一倍，以新步长计算下一步并计算误差。

（3）如步长小于某一下限h min 时不再减半，以免增加模拟时间，舍入误差过大。

（4）如步长大于某上限h min （指打印间隔），则不再加大步长。

龙格—库塔—默森法虽然增加了计算量，但在微机日益普及的今天，这一方法是值得广泛采用的。

采用这种方法的程序RKM 清单如下：
#include <stdio.h>
const int n=2;
inline float ABS(float input)
{if (input>=0.0f) return input;
else return -input;}
void main()
{void rkm(float *t,float y[2],float eps,float h,int *locp);
int locp=1,i;
float t=0.0f,eps=1.0E-8f,h=0.1f,y[2];
y[0]=1.0f; y[1]=1.0f;
printf("%-1.8f,%-1.8f,%-1.8f\n",t,y[0],y[1]);
for (i=1;i<=20;i++)
{rkm(&t,y,eps,h,&locp);
printf("%-1.8f,%-1.8f,%-1.8f\n",t,y[0],y[1]);}
}
void rkm(float *t,float y[2],float eps,float h,int *locp)
{void eql(float t,float y[2],float f[2]);
float w[5][n],hc,er;
int loc=1,lk=1,i;
hc=h/(float)(*locp);
label_1:
eql(*t,y,&w[2][0]);
for (i=0;i<=n-1;i++) w[0][i]=w[2][i]*hc/3.0f+y[i];
eql(hc/3.0f+*t,&w[0][0],&w[3][0]);
for (i=0;i<=n-1;i++) w[0][i]=(w[2][i]+w[3][i])*hc/6.0f+y[i];
eql(hc/3.0f+*t,&w[0][0],&w[3][0]);
for (i=0;i<=n-1;i++) w[0][i]=(w[2][i]+3.0f*w[3][i])*hc*0.125f+y[i];
eql(hc*0.5f+*t,&w[0][0],&w[4][0]);
for (i=0;i<=n-1;i++) w[0][i]=w[2][i]*hc*0.5f-w[3][i]*hc*1.5f+w[4][i]*hc*2.0f+y[i];
eql(hc+*t,&w[0][0],&w[3][0]);
for (i=0;i<=n-1;i++) w[1][i]=(w[2][i]+w[3][i]+4.0f*w[4][i])*hc/6.0f+y[i];
lk=1;
for (i=0;i<=n-1;i++)
{{er=ABS(w[0][i]-w[1][i])*0.2f;
if (ABS(w[1][i])>1.0f) er=er/ABS(w[1][i]);
if (er>=eps)
{hc=hc*0.5f;
*locp=(*locp)*2;
loc=2*loc;
goto label_1;}
if (er*64.0f>eps) lk=0;}}
*t=(*t)+hc;
for (i=0;i<=n-1;i++) y[i]=w[1][i];
loc++;
if (loc>=*locp) return;
if ((lk==0)||(loc%2==1)) goto label_1;
hc=2.0f*hc;
loc=loc/2;
*locp=*locp/2;
goto label_1;
}
void eql(float t,float y[2],float f[2])
{f[0]=1.0f/(y[1]-t);
f[1]=1.0f-1.0f/y[0];
}
使用RKM 程序计算微分方程组
1
1
1
1
2
21+-=-=y dt dy t
y dt dy
当2~0=t 时，1|01=-t y 1|02=-t y
步长取0.1
程序中的变量与数组为：
N 方程个数
T 模拟时间
H 积分区间
Y(2) 存放因变量初值与积分结果
EPS 允许的误差
LOCP 表示在为h 的积分区间等分数
W(N,5) 工作单元，二维数组
运行以上程序得到如下，表2.2引出了精确解与R 法与RKM 法的比较
表2.2
t 精确解 RK(H=0.1) MS(H=0.1,EPS=10e-8)
0 1.00000000 1.00000000 1.00000000
1.00000000 1.00000000 1.00000000
0.2 1.22140276 1.221402204 1.22140276
1.01873075 1.01873122 1.01873075
0.4 1.49182470 1.49182294 1.49182470
1.07032004 1.0732081 1.07032004
0.6 1.82211880 1.8211559 1.82211880
1.14881163 1.14881258 1.14881163
0.8 2.22554093 2.71827390 2.22554093
1,24932896 1.24932999 1.24932896
1.36787944 1.36788049 1.36787944
2.0 7.38905610 7.38901348 7.38905611
2.13533528 2.13533604 2.13533528
2.7 线性多步法
上几节我们讨论了单步法，即在计算y n+1值时，只要知道前一步的y n 和f n 的值就 可以进行计算，而每往后计算一步都可以使用前面已经算出的结果。

多步法就不同了，在计算y n+1时，不仅要用到前一步y n 和f n 的值, 还要使用y n+1, y n+2…及f n-1, f n-2…的值，一般而言，充分利用前面多步信息来计算y n+1，可期望加快模拟速度，也会使精度得到提高，这显然比单步法要好一些，这就是构成多步法计算的基本思想。

线性多步法中以亚当姆斯（Adams ）法使用得最为普遍，下面介绍这种方法。

欧拉法在计算函数y(t)时，在t n -t n-1中积分，用矩形公式求得。

其过程和示意图在前面已经给出。

欧拉法计算误差比较大，为了提高计算精度，将矩形改为梯形, 即在求积分时使用下式
h
y t f y t t y t f y dt
y t f y t y n n n n n n n n t t n n n n ),( )
)(,( ),()(111
+=-+=+=-+⎰+(2.16) )(2 )],(),([2
),(1111
n n n n n n t t f f h y t f y t f h dt y t f n n +=+≈+++⎰+
则整个数值积分公式变成
[]),(),(2
)(1111n n n n n n n y t f y t f h y y t y ++=≈++++ (2-17) 在(2-17)式中，为求y n+1值时，要用到y n+1本身，而这是一个未知数，故此式称为二阶隐式亚当姆斯公式，也称为梯形公式。

)],(),([2 )0(11)1(1n n n n n n y t f y t f h y y ++=+++
Adams 法的几何意义可用下图表示
f f n
f
f 0=f(t 0
t 0 t 1 t n t n+1
图2-4 用梯形公式进行数值积分
在求1+n y ，时由于要用到1+n y 本身，那就要用迭代法求解，即先猜出一个初值)
0(1+n y ，然后用下式求出)1(1+n y
再用)1(1+n y 求出)2(1+n y
[]
),(),(21)1(11)2(1-+++++=n n n n n n y t f y t f h y y (2.18) 如此重复下去，直到)(1n n y +与)1(1-+n n y 的差很小为止，此时可以认为1)(1++=n n n y y ，其截断误差为
2h
还有二阶亚当姆斯显式公式，其可表示如下
)],(),([2
)0(11)1(1n n n n n n y t f y t f h y y ++=+++19)-(2 )],(),(3[2
)(1111--++-+≈n n n n n n y t f y t f h y t y 20)-(2 )],(),(,([11101111+-+--++-++++=k n k n k n n n n n n y t f B y t f B y t f B h y y
表2-1到出了各阶亚当姆斯法计算公式中的系数值，这是一个多步公式，为了计算y n+1,不仅要知道yn,还要知道yn-1,yn-2,…,因此这些值都要保存，阶数越搞，要用到的值越多，
表2.1亚当姆斯系数表
在隐式多步法中，一般先用显式多步法计算处值，然后用隐式多步法作一次校正计算，这样使计算过程得到简化。

例如，用一阶显式计算初值，即
以上计算公式与二阶龙格—库塔法公式相同。

这种在多步法公式中先用显式公式估计一个值的方法 叫预估—校正法，它的优点是在精度相同（四阶），步长h 相同时，预估—校正法，每步计算两次右函数，而龙格—库塔法 要计算四次右函数，因此这种方法应用较为普遍。

四阶亚当姆斯预估—校正法公式为：
预估 )9375955(243211---+-+-+
=n n n n n F n f f f f h y y 校正 )5199(24
2111--+++-++=n n n F n n C n f f f f h y y 亚当姆斯的显示公式与隐式公式的特点如下：
（1）相同阶数的隐式公式的系数值比显式公式小（一阶例外），因而隐式公式比显式 公式精确得多．
（2）隐式公式的稳定性亦比显式的好。

下表给出了Adams 法的稳定区。

对于同样阶次，隐式的稳定域比显式的要大．随着阶数的增加，Adams 法的稳定域逐步减小．
（3）从计算上看，显式比隐示计算量小。

隐示公式需要计算f n+1,通常需要用Adams 显式公式对它提供一个首次近似值y n+1。

这种将显式公式和隐示公式联合起来使用以改进稳定性和精度的方法就是上面介绍的预估—校正法。

亚当姆斯法每步只需要计算一次 f 函数值，因而计算量较小，这是线性多步法的共同优点。

但它需要知道多个出发值才能计算，而微分方程只能提供一个初值，显然多步法不象单步法那样可以自启动，它必须用其他方法先获得所求时刻以前多步的解。

为了保证出发值的计算精度，一般采用比较高阶的单步法，以较小步长计算。

此外，多步法改变步长不能象21)
-(2 ),()0(1n n n n y t hf y y +=+22)-(2 ),(),([2
)0(11+++++=n n n n n n y h t f y t f h y y
单步法那样方便。

下面介绍一个使用亚当姆斯法进行模拟的程序，此程序可用于单输入单输出和多输入多输出线性系统的模拟，在模拟时面向状态方程。

线性定常系统的状态空间表达式包括两个矩阵方程，
第一个方程由n 个一阶微分方程组成，称为状态方程，第二个方程是一个线性代数方程，称为输出方程，式中X(t)为n ＊l 的状态向量，u(t)为m ＊1的控制向量。

Y(t)为l ＊l 的输出向量，A 为n ＊n 阶状态矩阵，由控制对象的参数决定，B 为n ＊m 的控制矩阵，C 为l*n 的输出矩阵，D 为l ＊m 和直接传输矩阵。

该程序用的计算方法为预估—校正线性多步法，首先用亚当姆斯显示公式计算预估值，再用亚当姆斯隐式公式计算校正值，其计算公式为
预估： )9375955(243211---+-+-+
=n n n n n F n X X X X h X X 校正： )5199(24
2111--+++-++=n n n F n n C n X X X X h X X 输出： DV CX Y C n n +=++11
Adambm 函数流程图如下：
⎪⎩
⎪⎨⎧=+=+=)0()()
()()()()()(0X t X t Du t CX t Y T Bu t AX t X
图2.5 Adambm函数流程图
程序中一些主要数组和变量如下：
X n x维状态向量
Y n y维输出向量
A n x *n y维系数矩阵
B n x *n u维输入系数矩阵
C n y *n x维与状态有关的输出矩阵
D n y *n u维输入有关的输出矩阵
U n u维控制输入
X(20),DX(20) 分别存放状态变量及导数的数组
U(6),Y(6) 分别存放输入与输出变量的数组
WORK(4,20) 存放状态变量四步导数的数组
TMAX,DT,PRDT 为模拟时间总长度，积分步长，打印步长
NX,NU,NR分别为状态变量、输入变量和输出变量的维数
NXM,NUM,NRM 分别为程序中允许的状态变量、输入变量和输出变量的最大维数。

MM,NP 为打印间隔数和打印点数。

模拟程序清单如下：
主程序：
#include "stdio.h"
#include "stdlib.h"
#include "iostream.h"
void adambm(float[][20],float[][6],float[][20],float[][6],float*,
int,int,int,float,float*,float *,float&,int&);。