连续系统的计算机模拟.doc

《连续系统的计算机模拟.doc》由会员分享，可在线阅读，更多相关《连续系统的计算机模拟.doc（36页珍藏版）》请在三一办公上搜索。

1、 第2章 连续系统的计算机模拟 本章讨论连续系统的模拟技术,由于这类系统中状态随时间连续动态地变化，常常具有一定的规律，故可用一些数学方程来描述，这些方程就是系统的数学模型，通常以微分方程、代数方程为多见。下面将介绍利用数值积分法对连续系统进行数字模拟的基本原理和具体方法，并给出数值积分法中几个常用的算法以及实现这些算法的计算程序，最后介绍两个建模实例。 数值积分法不仅方法种类多，而且有较强的理论性，本章由浅入深地介绍几种常见的数值解法。主要为单步法中的四阶龙格-库塔法与默森法和多步法中的亚当斯法。 使用数字计算机对连续系统进行模拟，首先必须将连续系统离散化，并将它转化为差分方程，以建立所谓的

2、模拟系统的数学模型。描述连续系统动态特征的数学模型是多种多样的，除微分方程外，还有传递函数、结构图及状态方程等，由于篇幅所限，本书不讨论后两种方法。 建立系统的模拟模型之后，就要选择计算机语言（也叫算法语言）编写系统模拟程序，在计算机上运行，将结果保留在数据文件中以待传输和处理。由于模拟的目的不同，可以选用不同的模拟模型和算法，其特点是运算精度高，对于不同的计算机，字长一般在16位-72位之间，也可采用浮点运算和双精度运算，其精度一般可达千万分之一到百万分之一。当然结果的精度与所选的算法有关。这可以根据实际需要选择机器、算法和模拟的步长。数字计算机储存容量大，可进行各种运算，在以前认为是不可能

3、解决的问题，利用数字计算机都可容易地或有可能得到解决。本章介绍的方法适应性较强，应用也十分广泛。数字计算机上还有各种功能的软件包（即一些子程序），用户可以稍加修改或不经修改就可以用于自己的模拟程序中，解决自己的实际问题，使用非常方便。 21 欧拉（Euler）法 在讨论连续系统的计算机模拟之前，让我们先看一个化学反应的例子，通过这个例子我们可以看到怎样使用数字计算机模拟一个实际问题，虽然介绍的是欧拉法，但是分析问题的思路同样适用与其它数值积分法。 当两种物质A和B放到一起产生化学反应时，产生第三种物质C，一般一克A与一克B结合产生2克的C物质，形成C 的速率与A 和B的数量乘积成正比，同样C也

4、可分解为A和B，C的分解速率正比与C的数量，即在任何时刻，如果a,b,和c是化学物质A,B,和C的数量，即在任何时刻，如果a,b,和c是化学物质A,B,C的数量，它们的增加和减少的速度服从下列微分方程。 （2.1） 其中K1和K2 是比例常数（一般而言这些比例常数会随温度和压力发生变化，但在模拟过程中，为了简化模型，一般不允许其变化，故一律视为常数）。在给出常数K1 和K2 值以及A和B的数量（C=0）后，我们希望能确定有多少C物质产在出来，这种化学反应速率的决定在化学工业上是有意义的。 模拟该系统的一个直接的方法是在t0时开始，使t以t间隔增加。假定化学量在t时间步长内不变，而只能在t结束的

5、瞬间发生变化，这样在每个t结束时的A（或B或C）的数量就可以从t开始时的数值由下式求出 (22)同样的方程b(t+t)和C(t+t),也可写出。假定模拟周期为T,可将T分成N个小的时间步长t,及 TNt 在时间为零时，我们知道a(0)、b(0)和c(0)的初始值,从这些初始值及常数K1和K2值出发，就可以计算出t时间内的化学量的变化值。 a(t)=a(0)+K2C(0)-K1a(0)b(0)t b(t)=b(0)+K2C(0)-K1a(0)b(0)t c(t)=c(0)+2K1a(0)b(0)-2K2 c(0)t 使用这些值，又可计算系统的下一个状态，即2t时的状态。 a(2t)=a(t)+K

6、2c (t)-K1a(t)b(t)t b(2t)=b(t)+K2c(t)-K1a(t )b(t )t c(2t)=c(t )+2K1a(t )b(t )-2K2 c(t)t 使用2t时的系统状态，又可写出3t时的状态，依此类推，以t为间隔，进行N步计算，就可写出周期T的系统状态得到所期望的结果，这个过程可用图2.1表示。 输入K1，K2，a(0),b(0),c(0) T, t和N I0 II1 保存a(I, t),b(I, t) c(I, t) Y N I=N? 打印a,b,c 模拟完毕 图2.1 化学反应例子模拟程序框图 模拟程序使用C语言编写，程序中的初值如下：K1=0.008/g.min

7、 K20.002/min, a(0)=100g, b(0)=50g C(0)=0, T=5mins, t=0.1min, N=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=1;i53;i+) printf(%2d,i); printf(%10.2f,%10.2f,%10.2f,%10.2fn,t,Ai

8、,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; 运行此程序，输出模拟结果如下： 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

9、.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.4510 0.90, 76.05, 26.05, 47.8911 1.00, 74.48, 24.48, 51.0412 1.10, 73.03, 23.03, 53.9413 1.20, 71.70, 21.70, 56.6114 1.30, 70.46, 20.46, 59.0715 1.40, 69.32, 19.32, 61.3616 1.50, 6

10、8.26, 18.26, 63.4817 1.60, 67.28, 17.28, 65.4418 1.70, 66.36, 16.36, 67.2819 1.80, 65.51, 15.51, 68.9920 1.90, 64.71, 14.71, 70.5921 2.00, 63.96, 13.96, 72.0822 2.10, 63.26, 13.26, 73.4823 2.20, 62.60, 12.60, 74.7924 2.30, 61.99, 11.99, 76.0325 2.40, 61.41, 11.41, 77.1826 2.50, 60.86, 10.86, 78.2727

11、 2.60, 60.35, 10.35, 79.3028 2.70, 59.87, 9.87, 80.2729 2.80, 59.41, 9.41, 81.1830 2.90, 58.98, 8.98, 82.0431 3.00, 58.57, 8.57, 82.8632 3.10, 58.19, 8.19, 83.6333 3.20, 57.82, 7.82, 84.3634 3.30, 57.48, 7.48, 85.0535 3.40, 57.15, 7.15, 85.7036 3.50, 56.84, 6.84, 86.3237 3.60, 56.55, 6.55, 86.9138 3

12、.70, 56.27, 6.27, 87.4639 3.80, 56.00, 6.00, 87.9940 3.90, 55.75, 5.75, 88.5041 4.00, 55.51, 5.51, 88.9742 4.10, 55.29, 5.29, 89.4343 4.20, 55.07, 5.07, 89.8644 4.30, 54.86, 4.86, 90.2745 4.40, 54.67, 4.67, 90.6646 4.50, 54.48, 4.48, 91.0347 4.60, 54.31, 4.31, 91.3948 4.70, 54.14, 4.14, 91.7349 4.80

13、, 53.98, 3.98, 92.0550 4.90, 53.82, 3.82, 92.3551 5.00, 53.68, 3.68, 92.6552 5.10, 53.54, 3.54, 92.93 以上例子的算法就是数值积分法中最简单的一种方法叫作欧拉法, 从这个例子中我们可以看出使用计算机解题的过程，由于欧拉法本身的特点，决定了其计算精度差, 现在几乎无人在实际工作中使用，但它导出简单， 能说明构造数值解法一般计算公式的基本思想, 模拟程序也清晰易懂，故人们乐意用它做为构造数值解法的入门例子。其一般解法如下：设给定微分方程(2.3)在区间(tn,tn+1)上求积分,得 y(tn+1)=

14、y(tn)+f(t,y)dt (2.4)如积分间隔足够小，使得tn与tn+1之间的f (t,y)可近似的看成常数f (tn,yn), 就可以用矩形面积近似地代替在该区间上的曲线积分，于是在tn+1时的积分值为 (2.5) 将上式写成以下差分方程形式: （2.6）这就是欧拉公式。它是一个递推的差分方程，任何一个新的数值解yn+1均是基于前一个数值解以及它的导数f(tn,yn)值求得的。只要给定初始条件y0及步长h就可根据f(t0,y0)算出y1的值，再以y1算出y2，如此逐步算出y3, y4,，直到满足所需计算的范围才停止计算。欧拉法的基本思路是把连续的时间t分割成等间隔的t, 在这些离散的时刻

15、算得函数值，根据这些值在函数图上可得到一条折线，所以欧拉法又叫折线法，其特点是分析方法简单，计算量小，但计算精度低（后面将讨论欧拉法与其它方法的比较）。下图为欧拉折线法的几何意义。 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 欧拉法的几何意义如果用梯形面积来代替一个小区间的曲线积分，就可克服以小矩形计算的缺点，提高精度，梯形法计算公式为(2.7)上式为隐式公式，因为公式右端含有yn+1, 这是未知的待求量，故梯形法不能自行启动运算，而要依赖于其它算法的帮助，比如说用欧拉公式法求出初值，算出的近似值

16、，然后将其带入微分方程，计算的近似值，最后利用梯形公式反复迭代。如迭代一次后就认为求得了近似解，这就是改进的欧拉法，其公式为 预估 （2-8） 校正 （2-9）上式中第一个为预估公式，第二个为校正公式。通常这种方法称为预估矫正法。在校正公式中计算了两点的斜率，再求其平均值，故计算量比欧拉法要大些。23 数值积分的几个基本概念 1. 单步法与多步法 数值积分法都用递推公式求解, 而递推公式是步进式的, 当由前一时刻的数值yn就能计算出后一时刻的数值yn+1时, 这种方法称为单步法, 它是一种能自启动的算法, 欧拉法是单步法. 反之, 如果求yn+1时要用到 yn, yn-1, yn-2 ,等多个

17、值, 则称为多步法,由于多步法在一开始时要用其它方法计算该时刻前面的函数值, 所以它不能自启动,以上所讲的改进的欧拉法就是一个多步法的例子。 2 显式与隐式 在递推公式中, 计算yn+1时所用到的数据均已算出的计算公式称为显示公式. 相反,在算式中隐含有末知量yn+1 的则称为隐含公式. 这需用另一个显示公式估计一个值, 再用隐式公式进行迭代运算, 这就是预估-校正法, 这种方法不能自启动. 用单步法与显示公式在计算上比用多步法和隐式公式方便. 但有时要考虑精度与稳定性时, 则必须采用隐式公式运算. 3 截断误差 一般使用台劳级数来分析数值积分公式的精度. 为简化分析, 假定前一步得的结果yn

18、是准确的, 则用台劳级数求得tn+1处的精确解为 (2.10)与前面的递推公式比较, 欧拉法是从以上精确解中只取前两 项之和来近似计算yn+1的 , 由这种方法单独一步引进的附加误差通常称着局部截断误差，它是该方法给出的值与微分方程解之间的差值，故又称为局部离散误差。欧拉法的局部截断误差为 不同的数值解法，其局部截断误差也不同。一般若截断误差为，则方法称为r阶的，所以方法的阶数可以作为衡量方法精确度的一个重要标志。欧拉法是一阶精度的数值解法，而改进的欧拉法(梯形法)相当于取台劳级数的前三项，故其解为二阶精度。 4 舍入误差 由于计算机的字长有限, 数字不能表示得完全精确, 在计算过程中不可避免

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

20、, 均可分为特解与通解两部分. 与稳定性有关的为通解部分, 它取决于该差分方程的特征根是否满足稳定性要求。例如, 为了考虑欧拉法的稳定性，可用检验方程y=y。其中为方程的特征根。对此方程有 要使该微分方程稳定，须使下式成立 |1+h|1有此可得 |h|2或h8 精度阶数234456n-2可见精度的阶数与计算函数值f 的次数之间的关系并非等量增加。 其精度最低。 如果将h取得很小，能否用欧拉公式得到很高的精度呢？从理论上讲是可以的，但是实际上由于计算机字长有限，在计算中有舍入误差。步长h越小，计算的步数越多，舍入误差的积累会越来越大，故用欧拉法时不能用很小的h。 2.5 四阶龙格库塔法模拟程序及

21、应用目前已有多种四阶龙格库塔法模拟（仿真）软件包推出，虽然子程序稍有不同，但总的结构还是相同的。对于连续系统的龙格库塔法的计算机程序，其大致结构如下图所示： 主 函 数输 入 及 运 行 输 出 初 始 化 微分方程 数值积分 显 示 打印 绘图 右函数 图2.3 龙格库塔法程序简要框图 图2.3中每一个程序模块的功能为：1. 主函数：主函数实现对整个模拟系统的控制, 这是通过调用各个子程序实现的。 2. 输入及初始化函数：主要对系统参数输入初值，例如模拟时间，积分步长，方程阶数等。这可以从键盘输入，也可从数据文件输入。3运行模块: 在运行子程序中, 将反复调用数值积分和方程右函数模块,进行迭

22、代计算，可以每计算一步便显示一次结果，也可以计算数步显示一次结果，屏幕的显示使用户可以随时监视系统运行的进程，以便人工干预。4输出子程序：按要求输出打印结果，可以打印，也可以绘图，视需要而定。四阶龙格库塔法公式前面已经给出, 现在再将它写成可编程的向量形式式中，i =1,2,3,.N, N为方程阶数。程序求得，这些分别为变量的导数，这样上式变为：Y(I)=YS(I)+0.166667* 2.0*(RK1（I）+RK3(I)+4.0*RK2(I)+D(I)*DT 其中 RK1(I)=D(I)*0.5DT RK2(I)=D(I)*0.5DT RK3(I)=D(I)*DT代入上式： Y(I) =YS

23、(I) +0.166667*DT*D (I)+2.0*D (I)+2.0*D (I)+D (I) 即 yn+1 =yn +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*/#inclu

24、de stdio.h#include conio.h#include math.h#define NORDER 3#define NPARAM 2float yNORDER,dNORDER,pNPARAM,t;void output(void)int j; printf(%9s,time); for(j=0;jNORDER;j+) printf( y%d,j); putchar(n);void difq(void) d0=-p0*y0*y1; d1=p0*y0*y1-p1*y1; d2=p1*y1; void run(float dt) int i; float rk3NORDER,ysNOR

25、DER; difq(); for(i=0;iNORDER;i+) rk0i=di*0.5*dt; ysi=yi; yi=ysi+rk0i; t+=0.5*dt; difq(); for(i=0;iNORDER;i+) rk1i=di*0.5*dt; yi=ysi+rk1i; difq(); for(i=0;iNORDER;i+) rk2i=di*dt; yi=ysi+rk2i; t+=0.5*dt; difq(); for(i=0;iNORDER;i+) yi=ysi+1.0/6*(2*(rk0i+rk2i)+4*rk1i+di*dt); void main(void)char c; do f

26、loat dt,tmax,co,tnext,tol; int j; printf(This program will find the root of:n); printf( y0(t)=-p0*y0*y1n); printf(&y1(t)=p0*y0*y1-p1*y1n); printf(&y2(t)=p1*y1n); for(j=0;jNORDER;j+) printf(Now,please input the first y%d:,j); scanf(%f,&yj); for(j=0;jNPARAM;j+) printf(Please input p%d:,j); scanf(%f,&p

27、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;kNORDER;k+

28、)printf(%10d,(int)(yk+0.5); putchar(n); run(dt); printf(nRun this program againY:); c=getche(); putchar(n); while(c=r|c=Y|c=y);程序中所用的变量和数组说明如下：Y(20) 状态变量数组G(20) 状态变量的一阶导数P(20) 存数系统参数，本例中仅一个参数，P(1)=P1=0.1TMAX 模拟时间DT 积分步长，即前面的hNP 打印点数NORDER 系统阶数NPARAM 系统参数个数OUTPUT 打印输出控制变量，.TRUE.为打印 FALSE 为不打印INIT 打印表

29、头控制变量A(8) 8个参数变量CO 打印时间间隔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的选择十分

30、重要，h太大不能达到预期的精度要求，h太小则增加了模拟时间，且有可能影响计算精度，要克服这一问题，就要采用变步长方法，使计算量尽可能的小，且精度也合乎要求。步长的自动控制是根据每一步的误差的大小来实现的。为了得到每一步的局部误差的估计值，可以采用两种不同阶次的递推公式（一般采用低一阶的RK公式，同时计算yn+1并取两者之差），要使计算量最少，可以选择系数，要求两个公式中的Ki相同，使中间结果对两种RK公式都适用，则这两个公式的计算结果之差就作为误差估计。1957年Merson（默森）给出了一个四阶变步长的方法，称为龙格库塔默森法。 其四阶公式为 yn+1 =yn +h (K1+4K4+K5)/6 (2.15) K1 =f (tn, yn )K2 =f (tn +h/3, yn +(h/3)*K1)K3 =f (tn+h/3, yn+(h/6)(K1+K2)K4 =f (tn+h/3, yn+(h/8)(K1+3K3) K5 =f (tn+h, yn+ (h/2)(K1-3K2+4k4)计算估计误差的三阶公式如下 其绝对误差为这一算法是四阶精度，三阶误差，故称为RKM34法，目前使用较多。其稳定域和一