第
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的数量，它们的增加和减少的速度服从下列微分方程。
da ?K2C?K1abdtdb ?K2C?K1abdtdc ?2K1ab?2K2Cdt （2.1）
其中K1和K2 是比例常数（一般而言这些比例常数会随温度和压力发生变化，但在模拟过程中，为了简化模型，一般不允许其变化，故一律视为常数）。在给出常数K1 和K2 值以及A和B的数量（C=0）后，我们希望能确定有多少C物质产在出来，这种化学反应速率
的决定在化学工业上是有意义的。
模拟该系统的一个直接的方法是在t＝0时开始，使t以Δt间隔增加。假定化学量在Δt时间步长内不变，而只能在Δt结束的瞬间发生变化，这样在每个Δt结束时的A（或B或C）的数量就可以从Δt开始时的数值由下式求出
a(t??t)?a(t)?da(t)?t dt (2．2) 同样的方程b(t+Δt)和C(t+Δt),也可写出。
假定模拟周期为T,可将T分成N个小的时间步长Δt,及 T＝N·Δt
在时间为零时，我们知道a(0)、b(0)和c(0)的初始值,从这些初始值及常数K1和K2值出发，就可以计算出Δt时间内的化学量的变化值。
a(Δt)=a(0)+[K2?C(0)-K1a(0)·b(0)]·Δt b(Δt)=b(0)+[K2·C(0)-K1a(0)·b(0)]·Δt
c(Δt)=c(0)+[2K1·a(0)·b(0)-2K2 c(0)]·Δt 使用这些值，又可计算系统的下一个状态，即2Δt时的状态。
a(2Δt)=a(Δt)+[K2?c (Δt)-K1a(Δt)·b(Δt)]·Δt b(2Δt)=b(Δt)+[K2·c(Δt)-K1a(Δt )·b(Δt )]·Δt
c(2Δt)=c(Δt )+[2K1·a(Δt )·b(Δt )-2K2 c(Δt)]·Δt
使用2Δt时的系统状态，又可写出3Δt时的状态，依此类推，以Δt为间隔，进行N步计算，就可写出周期T的系统状态得到所期望的结果，这个过程可用图2.1表示。 输入K1，K2，a(0),b(0),c(0) T, Δt和N
I＝0
I＝I＋1
保存a(I, Δt),b(I, Δt) c(I, Δt) Y N I<=N? 打印a,b,c
模拟完毕 图2.1 化学反应例子模拟程序框图
模拟程序使用C语言编写，程序中的初值如下：
K1=0.008/g.min K2＝0.002/min, a(0)=100g, b(0)=50g C(0)=0, T=5mins, ⊿t=0.1min, N=50
其程序清单如下： #include
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(\
printf(\ 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
以上例子的算法就是数值积分法中最简单的一种方法叫作欧拉法, 从这个例子中我们
可以看出使用计算机解题的过程，由于欧拉法本身的特点，决定了其计算精度差, 现在几乎无人在实际工作中使用，但它导出简单， 能说明构造数值解法一般计算公式的基本思想, 模拟程序也清晰易懂，故人们乐意用它做为构造数值解法的入门例子。其一般解法如下：
设给定微分方程
?y?f(t,y)??y(0)?y0
(2.3)
在区间(tn,tn+1)上求积分,得
y(tn+1)=y(tn)+∫f(t,y)dt (2.4)
如积分间隔足够小，使得tn与tn+1之间的f (t,y)可近似的看成常数f (tn,yn), 就可以用矩形面积近似地代替在该区间上的曲线积分，于是在tn+1时的积分值为
y (tn?1) yn?f(tn,yn)h?yn?1 (2.5)
将上式写成以下差分方程形式:
yn?1?yn?h?fn n?1,2,3,? （2.6）
这就是欧拉公式。它是一个递推的差分方程，任何一个新的数值解yn+1均是基于前一个数值解以及它的导数f(tn,yn)值求得的。只要给定初始条件y0及步长h就可根据f(t0,y0)算出y1的值，再以y1算出y2，如此逐步算出y3, y4,…，直到满足所需计算的范围才停止计算。欧拉法的基本思路是把连续的时间t分割成等间隔的Δt, 在这些离散的时刻算得函数值，根据这些值在函数图上可得到一条折线，所以欧拉法又叫折线法，其特点是分析方法简单，计算量小，但计算精度低（后面将讨论欧拉法与其它方法的比较）。下图为欧拉折线法的几何意义。
Y f(tn,yn) f(tn+1,yn+1)
欧拉折线
h y=∫f(t,y(t))dt+c 0 tn-1 tn tn+1 tn+2 图2.2 欧拉法的几何意义
如果用梯形面积来代替一个小区间的曲线积分，就可克服以小矩形计算的缺点，提高精度，梯形法计算公式为
连续系统的计算机模拟
