《连续系统的计算机模拟》由会员分享，可在线阅读，更多相关《连续系统的计算机模拟（35页珍藏版）》请在金锄头文库上搜索。
1、第2章连续系统的计算机模拟本章讨论连续系统的模拟技术，由于这类系统中状态随时间连续动态地变化，常常具有 一定的规律，故可用一些数学方程来描述，这些方程就是系统的数学模型，通常以微分方程、 代数方程为多见。下面将介绍利用数值积分法对连续系统进行数字模拟的基本原理和具体方 法，并给出数值积分法中几个常用的算法以及实现这些算法的计算程序，最后介绍两个建模实例。数值积分法不仅方法种类多，而且有较强的理论性， 本章由浅入深地介绍几种常见的数值解法。主要为单步法中的四阶龙格-库塔法与默森法和多步法中的亚当斯法。使用数字计算机对连续系统进行模拟，首先必须将连续系统离散化，并将它转化为差分方程，以建立所谓的模
2、拟系统的数学模型。 描述连续系统动态特征的数学模型是多种多样的， 除微分方程外，还有传递函数、结构图及状态方程等，由于篇幅所限，本书不讨论后两种方法。建立系统的模拟模型之后，就要选择计算机语言（也叫算法语言）编写系统模拟程序， 在计算机上运行，将结果保留在数据文件中以待传输和处理。由于模拟的目的不同， 可以选用不同的模拟模型和算法，其特点是运算精度高，对于不同的计算机，字长一般在16位-72位之间，也可采用浮点运算和双精度运算，其精度一般可达千万分之一到百万分之一。当然结果的精度与所选的算法有关。这可以根据实际需要选择机器、算法和模拟的步长。 数字计算机储存容量大，可进行各种运算，在以前认为是
3、不可能解决的问题，利用数字计算机都可容易地或有可能得到解决。本章介绍的方法适应性较强，应用也十分广泛。数字计算机上还 有各种功能的软件包（即一些子程序），用户可以稍加修改或不经修改就可以用于自己的模拟程序中，解决自己的实际问题，使用非常方便。2 . 1 欧拉（Euler ）法在讨论连续系统的计算机模拟之前， 让我们先看一个化学反应的例子， 通过这个例子我 们可以看到怎样使用数字计算机模拟一个实际问题，虽然介绍的是欧拉法， 但是分析问题的思路同样适用与其它数值积分法。当两种物质A和B放到一起产生化学反应时，产生第三种物质C, 一般一克A与一克B结合产生2克的C物质，形成C的速率与A和B的数量乘积
4、成正比，同样 C也可分解 为A和B, C的分解速率正比与 C的数量，即在任何时刻，如果a,b,和c是化学物质 A,B,和C的数量，即在任何时刻，如果 a,b,和c是化学物质A,B,C的数量，它们的增加和减少的 速度服从下列微分方程。da dt db dtK2CK2CK1abK1ab(2.1 )2K ab 2K C dt 12其中K1和K2是比例常数（一般而言这些比例常数会随温度和压力发生变化，但在模拟过程中，为了简化模型，一般不允许其变化，故一律视为常数）。在给出常数 K1和K2值以及A和B的数量（C=0）后，我们希望能确定有多少 C物质产在出来，这种化学反应速率的决定在化学工业上是有意义的。
5、模拟该系统的一个直接的方法是在t=0时开始，使t以At间隔增加。假定化学量在At时间步长内不变，而只能在A t结束的瞬间发生变化，这样在每个At结束日的A(或B或C)的数量就可以从A t开始时的数值由下式求出a(t t) a(t) ddt) t(2同样的方程 b(t+ At)和C(t+At),也可写出。假定模拟周期为T,可将T分成N个小的时间步长A t,及T= N At出发，就可以计算出At时间内的化学量的变化值。.2)Ki和K2值在时间为零时，我们知道 a(0)、b(0)和c(0)的初始值，从这些初始值及常数a(At)=a(0)+K2? C(0)-Kia(0) - b(0) -Atb(At)
6、=b(0)+K2 C(0)-Kia(0) - b(0) -Atc(At)=c(0)+2K 1 - a(0)- b(0)-2K 2 c(0)A t使用这些值，又可计算系统的下一个状态，即 2At时的状态。a(2A t)=a(A t)+K 2?c ( At)-Kia(At)- b( At) A tb(2A t)=b(At)+K 2c( At)-K ia( At ) b( At ) Atc(2A t)=c( A t )+2K 1 a( A t ) b( A t )-2K 2 c( At) A t模拟程序使用K=0.008/g.minC(0)=0,使用2At时的系统状态，又可写出3At时的状态，依此类
7、推，以A t为间隔，进彳T N步计算，就可写出周期 T的系统状态得到所期望的结果，这个过程可用图2.1表示。图2.1化学反应例子模拟程序框图C语言编写，程序中的初值如下：K2=0.002/min,T=5mins,a(0)=100g, /t=0.1min,b(0)=50gN=50其程序清单如下：#include #include float k1,k2;static float A53,B53,C53,delta,t;void strut(int);main()int i;A1=100.0;B1=50.0;C1=0.0;t=0;delta=0.1;k1=0.008;k2=0.002;for(i=
8、1;i53;i+)printf(%2d,i);printf(%10.2f,%10.2f,%10.2f,%10.2fn,t,Ai,Bi,Ci);strut(i);return;void strut(int i)Ai+1=Ai+(k2*Ci-k1*Ai*Bi)*delta;Bi+1=Bi+(k2*Ci-k1*Ai*Bi)*delta;Ci+1=Ci+2.0*(k1*Ai*Bi-k2*Ci)*delta;t=t+delta;运行此程序，输出模拟结果如下：10.00,100.00,50.00,0.0020.10,96.00,46.00,8.0030.20,92.47,42.47,15.0640.30,
9、89.33,39.33,21.3450.40,86.52,36.52,26.9560.50,84.00,34.00,32.0070.60,81.72,31.72,36.5580.70,79.66,29.66,40.6990.80,77.77,27.77,44.45100.90,76.05,26.05,47.89111.00,74.48,24.48,51.04121.10,73.03,23.03,53.94131.20,71.70,21.70,56.61141.30,70.46,20.46,59.07151.40,69.32,19.32,61.36161.50,68.26,18.26,63.48
10、171.60,67.28,17.28,65.44181.70,66.36,16.36,67.28191.80,65.51,15.51,68.99201.90,64.71,14.71,70.59212.00,63.96,13.96,72.08222.10,63.26,13.26,73.48232.20,62.60,12.60,74.79242.30,61.99,11.99,76.03252.40,61.41,11.41,77.18262.50,60.86,10.86,78.27272.60,60.35,10.35,79.30282.70,59.87,9.87,80.27292.80,59.41,
11、9.41,81.18302.90,58.98,8.98,82.04313.00,58.57,8.57,82.86323.10,58.19,8.19,83.63333.20,57.82,7.82,84.36343.30,57.48,7.48,85.05353.40,57.15,7.15,85.70363.50,56.84,6.84,86.32373.60,56.55,6.55,86.91383.70,56.27,6.27,87.46393.80,56.00,6.00,87.99403.90,55.75,5.75,88.50414.00,55.51,5.51,88.97424.10,55.29,5
12、.29,89.43434.20,55.07,5.07,89.86444.30,54.86,4.86,90.27454.40,54.67,4.67,90.66464.50,54.48,4.48,91.03474.60,54.31,4.31,91.39484.70,54.14,4.14,91.73494.80,53.98,3.98,92.05504.90,53.82,3.82,92.35515.00,53.68,3.68,92.65525.10,53.54,3.54,92.93从这个例子中我们以上例子的算法就是数值积分法中最简单的一种方法叫作欧拉法可以看出使用计算机解题的过程，由于欧拉法本身的特点
13、，决定了其计算精度差，现在几乎无人在实际工作中使用，但它导出简单，能说明构造数值解法一般计算公式的基本思想，模拟程序也清晰易懂，故人们乐意用它做为构造数值解法的入门例子。其一般解法如下： 设给定微分方程y f(t,y) y(0) Vo (2.3) 在区间(t n, t n+1)上求积分,得y(tn+1)=y(t n)+ f f(t,y)dt(2.4)如积分间隔足够小，使得 tn与tn+1之间的f (t,y)可近似的看成常数f (tn,yn),就可以用矩形面积近似地代替在该区间上的曲线积分，于是在t n+1时的积分值为y(tn 1) yn f(tn，yn)h 九 1(2.5)将上式写成以下差分方程形式：ynlyn h fnn 1,2,3,*.6)这就是欧拉公式。它是一个递推的差分方程，任何一个新的数值解 yn+1均是基于前一个数值 解以及它的导数f(tn，yn)值求得的。只要名定初始条件yo及步长h就可根据f(to，yo)算出y1的值，再以y1算出y2,如此逐步算出 y3, y4,，直到满足所需计算的范围才停止计算。欧拉 法的
