连续系统的计算机模拟

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

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

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

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

5、使t以At间隔增加。假定化学量在At时间步长内不变，而只能在At结束的瞬间发生变化，这样在每个At结束日的A(或B或C)的数量就可以从At开始时的数值由下式求出a(t：t)=a(t)da.t(2同样的方程b(t+At)和C(t+At),也可写出。假定模拟周期为T,可将T分成N个小的时间步长At,及T=N-At.2)K1和K2值在时间为零时，我们知道a(0)、b(0)和c(0)的初始值，从这些初始值及常数出发，就可以计算出a(b(c(使用这些值,a(2b(2c(2At时间内的化学量的变化值。At)=a(0)+K2?C(0)-Kia(0)-b(0)-At)=b(0)+K2C(0)-Kia(0)-b

6、(0)-At)=c(0)+2Ki-a(0)-b(0)-2K2c(0)又可计算系统的下一个状态，即At)=a(At)=b(At)=c(AtAtAt2At时的状态。At)+K2?c(At)-Kia(At)b(At)At)+K2c(At)-Kia(At)-b(At)At)+2Kia(At)-b(At)-2K2c(AtAtAt)At2.1表不图2.i化学反应例子模拟程序框图C语言编写，程序中的初值如下：使用2At时的系统状态，又可写出3At时的状态，依此类推，以At为间隔，进彳TN步计算，就可写出周期T的系统状态得到所期望的结果，这个过程可用图模拟程序使用Ki=0.008/g.minC(0)=0,T=

7、5mins,K20.002/min,a(0)=100g,/t=0.1min,b(0)=50gN=50其程序清单如下：#include<stdio.h>#include<string.h>floatk1,k2;staticfloatA53,B53,C53,delta,t;voidstrut(int);main()inti;A1=100.0;B1=50.0;C1=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

8、,%10.2f,%10.2fn",t,Ai,Bi,Ci);strut(i);return;voidstrut(inti)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,89.33,39.33,21.3450.40,86.52,36.52,26.9

9、560.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.48171.60,67.28,17.28,65.44181.70,66.36,16

10、.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,9.41,81.18302.90,58.98,8.98,82.04313.00

11、,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.29,89.43434.20,55.07,5.07,89.86444.30,

12、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、值解法一般计算公式的基本思想，模拟程序也清晰易懂，故人们乐意用它做为构造数值解法的入门例子。其一般解法如下：设给定微分方程：v=f(t,y)y=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(tn1)ynf(tn，yn)h=yn.i(2.5)将上式写成以下差分方程形式：yn+=yn+hfnn=1,2,3，(2.6)这就是欧拉公式。它是一个递推的差分方程，任何一个新的数值解y

14、n+1均是基于前一个数值解以及它的导数f(tn,yn)值求得的。只要名定初始条件yo及步长h就可根据f(t0,y0)算出yi的值，再以yi算出y2,如此逐步算出y3,y4,，直到满足所需计算的范围才停止计算。欧拉法的基本思路是把连续的时间t分割成等间隔的At,在这些离散的时刻算得函数值，根据这些值在函数图上可得到一条折线，所以欧拉法又叫折线法，其特点是分析方法简单，计算量小，但计算精度低(后面将讨论欧拉法与其它方法的比较)。下图为欧拉折线法的几何意义。Y f(t图2.2欧拉法的几何意义如果用梯形面积来代替一个小区间的曲线积分，就可克服以小矩形计算的缺点，提高精度，梯形法计算公式为hyn1-yn

15、2f(tn，yn)f(tn1,yn1)=yn2(fnfni)-,汴(2(2.7)上式为隐式公式，因为公式右端含有yn+1,这是未知的待求量，故梯形法不能自行启动运算,而要依赖于其它算法的帮助，比如说用欧拉公式法求出初值，算出y(yn中)的近似值yF书，然后将其带入微分方程，计算fn书的近似值fnl=f(tn书，yP书)，最后利用梯形公式反复迭代。如迭代一次后就认为求得了近似解，这就是改进的欧拉法，其公式为预估yP 书=yn +hf(tn,yn)(2-8)校正yC书=yn+2f(tn.Yn)+fP(tn+.yP+)(2-9)h上式中第一个为预估公式，第二个为校正公式。通常这种方法称为预估矫正法。

16、在校正公式中计算了两点的斜率，再求其平均值，故计算量比欧拉法要大些。2. 3数值积分的几个基本概念1 .单步法与多步法数值积分法都用递推公式求解，而递推公式是步进式的，当由前一时刻的数值yn就能计算出后一时刻的数值yn+1时，这种方法称为单步法，它是一种能自启动的算法，欧拉法是单步法.反之，如果求yn+1时要用到yn,yn-1,yn-2,等多个值，则称为多步法，由于多步法在一开始时要用其它方法计算该时刻前面的函数值，所以它不能自启动，以上所讲的改进的欧拉法就是一个多步法的例子。2 显式与隐式在递推公式中，计算yn+1时所用到的数据均已算出的计算公式称为显示公式.相反，在算式中隐含有末知量yn+

17、1的则称为隐含公式.这需用另一个显示公式估计一个值，再用隐式公式进行迭代运算，这就是预估-校正法，这种方法不能自启动.用单步法与显示公式在计算上比用多步法和隐式公式方便.但有时要考虑精度与稳定性时，则必须采用隐式公式运算.3 截断误差一般使用台劳级数来分析数值积分公式的精度.为简化分析，假定前一步彳#的结果yn是准确的，则用台劳级数求得tn+1处的精确解为1一1(2.10)ylnkyghg)”，(«)Phn与前面的递推公式比较，欧拉法是从以上精确解中只取前两项之和来近似计算yn+1的，由这种方法单独一步引进的附加误差通常称着局部截断误差，它是该方法给出的值与微分方程解之间的差值，故又

18、称为局部离散误差。欧拉法的局部截断误差为O(h2),则方法称为r阶的,y(tn1)-y(tn)=O(h2)不同的数值解法，其局部截断误差也不同。一般若截断误差为所以方法的阶数可以作为衡量方法精确度的一个重要标志。欧拉法是一阶精度的数值解法，而改进的欧拉法(梯形法)相当于取台劳级数的前三项，故其解为二阶精度。4 舍入误差由于计算机的字长有限，数字不能表示彳#完全精确，在计算过程中不可避免地会产生舍入误差.舍入误差与步长h成反比，如计算步长小，计算次数多则舍入误差大.产生舍入误差的因素较多，除了与计算机字长有关外，还与计算机所使用的数字系统，数的运算次序及计算f(t,y)所用的程序的精确度等因素有

19、关.5 数值解的稳定性以上对连续系统的模拟用到了差分方程，把初始值代入递推公式进行迭代运算，如果原系统是稳定的，由数值积分法求得的解也应该是稳定的.但是，在计算机逐次计算时，初始数据的误差及计算过程中的舍入误差对后面的计算结果会产生影响.而且计算步长选择得不合理时，有可能使模拟结果不稳定.以下我们简要地讨论这一问题.差分方程的解与微分方程的解类似，均可分为特解与通解两部分.与稳定性有关的为通解部分，它取决于该差分方程的特征根是否满足稳定性要求。例如，为了考虑欧拉法的稳定性，可用检验方程y=Xyo其中入为方程的特征根。对此方程有yn1=yn'ynh=(11h)yn要使该微分方程稳定，须使

20、下式成立|1+入h|<1有此可得|入h|<2或h<2T要保证用欧拉法进行的计算是稳定的，积分步长h必须小于系统时间常数的两倍.所以积分步长的选择就要慎重，不能随意取任何值。2.4龙格贝塔(Runge-Kutta)法仍以(2.3)方程为例.已知y(t0)=y0,假设我t0开始以h增长，ti=t0+h,t1时刻为2一一一»yi=y(M+h),在b附近展开成台老级数，保留h项，则有.h3fdy；:f2yi=y()f(t0,y0)h()h(2-ii)2二ydt二t上式括号内下标为0,表示括号中的函数用t=to,y=yo代入，以下均同。假设上式的解可写成如下形式yi=yo(b

21、iKib2K2)h(2.i2)其中,Ki=f(10,yo)K2=f(t0C2h,y°aiKih)对K2式右端得函数在t=t0,y=y。处展开成台劳级数，保留h项，可得到：-汗汗K2：f(t。y。)(C2aiKi-)0h二t二y将Ki与K2代入(2-i2)式得:f;:fy1=y°,b)hf(to,y°),bzhf(to,y°),(C2aiKi)ohftFy将上式与(2-12)式比较，可得系数ai,bi,C2,b2方程如下:bi+b2=1,b2c2=二2ib2ai二二、2以上三个方程，四个未知数，因此有无穷个解，若限定bi=b2则可得一个解：ai=C2=i将

22、它们代入(2-i2)式可得出一组计算公式hyi=y。(KJ&)-其中Ki=f(t0,y0),k2=f(t0+h,y0+kih),写成一般的递推形式如下：hy(tni)：yni二yn2(KiK2)(2.i3)其中，Ki=f(tn,yn),K2=f(tn+h,yn十K1h),式中只取h,h2两项，而将h3以上的高阶项略去，所以此递推公式的截断误差正比于h的三次方，计算过程中只取h和h2两项，故此方法被称为二阶龙格-库塔法。一般使用的为四阶龙格-库塔公式，在展开台劳级数时保留h,h2,h3,和h4项,其截断误差为h5,四阶龙格库塔公式为y(tni)：yni=ynh(Ki2K22K3K4)(2

23、.i4)6式中：Ki=f(tn,yn)K2=f(tn+、+Ki)rhhK3=f(tn-,yn-K2)K4=f(tnh.ynhK3)四阶龙格一库塔法公式精度较高，故其应用较为普遍。下面再给出几组系数值设a1=1C2C2C21212b1 =0bibi1412b2b2b2二1二1读者可以写出其相应的四阶龙格一库塔法公式。龙格一库塔法的特点如下:1 .龙格库塔法有多组a1,C2b1b2值，故可有多种龙格库塔公式，常使用的有:yn书=yn+K?hK=f(tn,yn)K2=f(tn+,yn+K1h)所有龙格-库塔公式都有以下几个特点：(1)在计算yn+1时只用到yn,而不直接用yn-1,yn-2等项，也就

24、是计算后一步时只用到前一步的计算结果，故称为单步法，单步法可以自启动，且使用的存储量也小。(2)在计算过程中，步长h可以改变，这要由计算精度来定，但是在一步中，为计算若干个系数Ki(称成为龙格-库塔系数)，必须使用同一个步长ho(3).不论是几阶龙格-库塔公式，计算式中均为两部分组成。第一部分为上一步的结果yn,第二部分为h=tn-tn-1中对他y)的积分。它是步长h乘上各点斜率的加权平均。例如上面计算的四阶公式中，取四点的斜率：K1、K2、K3和K4,然后对K2和K3取两份，K1和K4各取一份，进行加权平均。3 .龙格一库塔法的精度取决于步长h的大小及求解的方法。实践证明：为达到同样的精度，

25、四阶方法的步长h可以比二阶方法的h大10倍，而四阶方法的每步计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小，使用的CPU机时也少。一般工程中四阶龙格库塔公式已能达到要求的精度，故不再使用更高阶的公式。下表示出了f的计算次数与精度阶数的关系。每一步计算f的次数234567n=>8精度阶数234456n-2可见精度的阶数与计算函数值f的次数之间的关系并非等量增加。4 .若在展开成台劳级数时，只取h这一项，而将h2以及h2以上的高次项略去，可得yn1=yn,hf(tn,Yn)这是欧拉公式，因此欧拉公式也可看成时一阶龙格-库塔公式，她的截断误差正比h2,其精度最低。如果将h取得很小，能否

26、用欧拉公式得到很高的精度呢？从理论上讲是可以的，但是实际上由于计算机字长有限，在计算中有舍入误差。步长h越小，计算的步数越多，舍入误差的积累会越来越大，故用欧拉法时不能用很小的ho2.5四阶龙格一库塔法模拟程序及应用目前已有多种四阶龙格一库塔法模拟(仿真)软件包推出，虽然子程序稍有不同，但总的结构还是相同的。对于连续系统的龙格一库塔法的计算机程序，其大致结构如下图所示：图2.3龙格一库塔法程序简要框图图2.3中每一个程序模块的功能为：1 .主函数：主函数实现对整个模拟系统的控制，这是通过调用各个子程序实现的。2 .输入及初始化函数：主要对系统参数输入初值，例如模拟时间，积分步长，方程阶数等。这

27、可以从键盘输入，也可从数据文件输入。3 .运行模块：在运行子程序中，将反复调用数值积分和方程右函数模块，进行迭代计算,可以每计算一步便显示一次结果，也可以计算数步显示一次结果，屏幕的显示使用户可以随时监视系统运行的进程，以便人工干预。4 .输出子程序：按要求输出打印结果，可以打印，也可以绘图，视需要而定。四阶龙格一库塔法公式前面已经给出，现在再将它写成可编程的向量形式1yy(Kfko2K°KJyi,n1yi,ni1i2i3i46Kilhfi(tn,y1n,y2n,yNn)Ki2= hfi(tnh2，y1n112K22, ,yNn 2 Kn2)Ki3rh111=hfi(tn2,yn2K

28、12,丫2n万心,，yNn2KN2)Ki4=hfi(tnh)1n*13/2n七3,yNn*N3)式中，i=1,2,3,.N,N为方程阶数。令DT=h,Ki1，Ki2，Ki3和Ki4分别换成D,D(I),D(3)(I)和口(4),D(I)为调微分子程序求得，这些分别为变量的导数，这样上式变为：Y(I)=YS(I)+0.166667*2.0*(RK1(I)+RK3(I)+4.0*RK2(I)+D(I)*DT其中RK1(I)=D(I)*0.5DTRK2(I)=D(I)*0.5DTRK3(I)=D(I)*DT代入上式：Y(I)=YS(I)+0.166667*DT*D(I)+2.0*D(I)+2.0*D

29、(I)+D(I)即yn+1=yn+h(K1+2K2+2K3+K4)设有一个微分方程：y(t)+P1y(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程序清单如下：主程序：/*p23example*/#include"stdio.h"#include"conio.h"#include"mat

30、h.h"#defineNORDER3#defineNPARAM2floatyNORDER,dNORDER,pNPARAM,t;voidoutput(void)intj;printf("%9s","time");for(j=0;j<NORDER;j+)printf("y%d",j);putchar('n');voiddifq(void)d0=-p0*y0*y1;d1=p0*y0*y1-p1*y1;d2=p1*y1;voidrun(floatdt)inti;floatrk3NORDER,ysNORDER;d

31、ifq();for(i=0;i<NORDER;i+)rk0i=di*0.5*dt;ysi=yi;yi=ysi+rk0i;t+=0.5*dt;difq();for(i=0;i<NORDER;i+)rk1i=di*0.5*dt;yi=ysi+rk1i;difq();for(i=0;i<NORDER;i+)rk2i=di*dt;yi=ysi+rk2i;t+=0.5*dt;difq();for(i=0;i<NORDER;i+)yi=ysi+1.0/6*(2*(rk0i+rk2i)+4*rk1i+di*dt);voidmain(void)charc;dofloatdt,tmax,

32、co,tnext,tol;intj;printf("Thisprogramwillfindtherootof:n");printf("y0'(t)=-p0*y0*y1n");printf("&&y1'(t)=p0*y0*y1-p1*y1n");printf("&&y2'(t)=p1*y1n");for(j=0;j<NORDER;j+)printf("Now,pleaseinputthefirsty%d:",j);scanf("

33、;%f",&yj);for(j=0;j<NPARAM;j+)printf("Pleaseinputp%d:",j);scanf("%f",&pj);printf("Pleaseinputthetimeofsimulation:");scanf("%f",&tmax);printf("Pleaseinputthetimeofonestep:");scanf("%f",&dt);printf("Now,thisisther

34、esult:n");co=5*dt;tol=0.0001*co;tnext=0;output();for(t=0;t<=tmax+tol;)intk;if(fabs(tnext-t)<tol)tnext+=co;printf("%10.4f",t);for(k=0;k<NORDER;k+)printf("%10d",(int)(yk+0.5);putchar('n');run(dt);printf("nRunthisprogramagainY:");c=getche();putchar(&#

35、39;n');while(c='r'|c='Y'|c='y');程序中所用的变量和数组说明如下：Y(20)状态变量数组G(20)状态变量的一阶导数P(1)=P1=0.1FALSEP(20)存数系统参数，本例中仅一个参数，TMAX模拟时间DT积分步长，即前面的hNP打印点数NORDER系统阶数NPARAM系统参数个数OUTPUT打印输出控制变量，.TRUE.为打印为不打印INIT打印表头控制变量A(8)8个参数变量CO打印时间间隔NOUT已打印点数计数器TNEXT打印点处的时间值SS剩下的打印点数NLIST输出结果对以上程序编译，运行，得到

36、如下结果TIMEY(1)Y(2)0.00000.00000.00000.50000.12040.46761.00000.44500.80081.50000.88630.92652.00001,33320.82472.50001.67880.53109.50001.66660.833310.00001.66660.83332.6 变步长的龙格一库塔法以上提到，步长h的选择十分重要，h太大不能达到预期的精度要求，h太小则增加了模拟时间，且有可能影响计算精度，要克服这一问题，就要采用变步长方法，使计算量尽可能的小，且精度也合乎要求。步长的自动控制是根据每一步的误差的大小来实现的。为了得到每一步的局部

37、误差的估计值，可以采用两种不同阶次的递推公式（一般采用低一阶的RK公式，同时计算yn+1并取Ki相同，使中间结果1957 年 Merson （默(2.15)两者之差），要使计算量最少，可以选择RK系数，要求两个公式中的对两种RK公式都适用，则这两个公式的计算结果之差就作为误差估计。森）给出了一个四阶变步长的方法，称为龙格一库塔默森法。其四阶公式为yn+1=yn+h(K1+4K4+K5)/6K1=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+(

38、h/2)(K1-3K2+4k4)计算估计误差的三阶公式如下n+=yn+?(3K1-9K3+12K4)6其绝对误差为一hE=yn1-yn1(2K1-9K38K4-七)6这一算法是四阶精度，三阶误差，故称为RKM34法，目前使用较多。其稳定域和一般RK4法相近，缺点是计算量稍大，每步计算5次f值，除用绝对误差外，也可用相对误差，最大相对误差RE为RE = max(Eyn 1 - yn在每一步计算后，先计算误差，根据误差的大小来控制步长，控制策略如下：(1)当误差ERR大于预先设定的最大允许误差Emax时，则缩小步长，一般将步长减半,并以新步长重新计算以后再比较。(2)当误差ERR小于预先定的最小允

39、许误差Emm时，步长增加一倍，以新步长计算下(3)如步长小于某一下限hmin时不再减半，以免增加模拟时间，舍入误差过大。(4)如步长大于某上限hmin(指打印间隔)，则不再加大步长。龙格一库塔默森法虽然增加了计算量，但在微机日益普及的今天，这一方法是值得广泛采用的。采用这种方法的程序RKM清单如下：#include<stdio.h>constintn=2;inlinefloatABS(floatinput)if(input>=0.0f)returninput;elsereturn-input;voidmain()voidrkm(float*t,floaty2,floateps

40、,floath,int*locp);intlocp=1,i;floatt=0.0f,eps=1.0E-8f,h=0.1f,y2;y0=1.0f;y1=1.0f;printf("%-1.8f,%-1.8f,%-1.8fn",t,y0,y1);for(i=1;i<=20;i+)rkm(&t,y,eps,h,&locp);printf("%-1.8f,%-1.8f,%-1.8fn",t,y0,y1);voidrkm(float*t,floaty2,floateps,floath,int*locp)voideql(floatt,floaty

41、2,floatf2);floatw5n,hc,er;intloc=1,lk=1,i;hc=h/(float)(*locp);label1:eql(*t,y,&w20);for(i=0;i<=n-1;i+)w0i=w2i*hc/3.0f+yi;eql(hc/3.0f+*t,&w00,&w3。);for(i=0;i<=n-1;i+)w0i=(w2i+w3i)*hc/6.0f+yi;eql(hc/3.0f+*t,&w00,&w3。);for(i=0;i<=n-1;i+)w0i=(w2i+3.0f*w3i)*hc*0.125f+yi;eql(h

42、c*0.5f+*t,&w0,&w40);for(i=0;i<=n-1;i+)w0i=w2i*hc*0.5f-w3i*hc*1.5f+w4i*hc*2.0f+yi;eql(hc+*t,&w00,&w30);for(i=0;i<=n-1;i+)w1i=(w2i+w3i+4.0f*w4i)*hc/6.0f+yi;lk=1;for(i=0;i<=n-1;i+)er=ABS(w0i-w1i)*0.2f;if(ABS(w1i)>1.0f)er=er/ABS(w1i);if(er>=eps)hc=hc*0.5f;*locp=(*locp)*2;l

43、oc=2*loc;gotolabel_1;if(er*64.0f>eps)lk=0;*t=(*t)+hc;for(i=0;i<=n-1;i+)yi=w1i;loc+;if(loc>=*locp)return;if(lk=0)|(loc%2=1)gotolabel_1;hc=2.0f*hc;10c=loc/2;*locp=*locp/2;gotolabel_1;voideql(floatt,floaty2,floatf2)f0=1.0f/(y1-t);f1=1.0f-1.0f/y0;使用RKM程序计算微分方程组dy_1dty2-t蛆二1dty1当t=02时，y1It=1y21t

44、_0=1步长取0.1程序中的变量与数组为:N方程个数T模拟时间H积分区间Y(2)存放因变量初值与积分结果EPS允许的误差LOCP表示在为h的积分区间等分数W(N,5)工作单元，二维数组运行以上程序得到如下，表2.2引出了精确解与R法与RKM法的比较表2.2t精确解RK(H=0.1)MS(H=0.1,EPS=10e-8)01.000000001.000000001.000000001.000000001.000000001.000000000.21.221402761.2214022041.221402761.018730751.018731221.018730750.41.491824701.

45、491822941.491824701.070320041.07320811.070320040.61.822118801.82115591.822118801.148811631.148812581.148811630.82.225540932.718273902.225540931,249328961.249329991.249328961.367879441.367880491.367879442.07.389056107.389013487.389056112.135335282.135336042.135335282.7 线性多步法上几节我们讨论了单步法，即在计算yn+1值时，只要知

46、道前一步的yn和fn的值就可以进行计算，而每往后计算一步都可以使用前面已经算出的结果。多步法就不同了，在计算yn+1时，不仅要用到前一步yn和fn的值，还要使用yn+1,yn+2及fn-1,fn-2的值，一般而言，充分利用前面多步信息来计算yn+1,可期望加快模拟速度，也会使精度得到提高，这显然比单步法要好一些，这就是构成多步法计算的基本思想。线性多步法中以亚当姆斯(Adams)法使用得最为普遍，下面介绍这种方法。欧拉法在计算函数y时，在tn-tn-1中积分，用矩形公式求得。tn1y(tn1)=yntf(t,y)dttnuyn.f(tn,yn)(tn4tn)=ynf(tn,Yn*其过程和示意图

47、在前面已经给出。欧拉法计算误差比较大，为了提高计算精度，将矩形改为梯形，即在求积分时使用下式tn1htf(t,y)dt：-f(tn1,yn1)f(tn,yn)tn2=h(fn1fn)(2.16)2则整个数值积分公式变成h(2-17)y(tn1)yn1=yn-f(tn1,yn1)f(tn,yn)1在(2-17)式中，为求yn+1值时，要用到yn+1本身，而这是一个未知数，故此式称为二阶隐式亚当姆斯公式，也称为梯形公式。(1)h(0).yn=yn2f(tn1，yn1)f(tn，yn)Adams法的几何意义可用下图表示用梯形公式进行数值积分图2-4f,时由于要用到 yn*本身，那就要用迭代法求解，即

48、先猜出一个初值(0)y n+ ，y：) = yn。£ (tn 1,yn0)1) f (tn,yn)然后用下式求出yn1)1再用y，1求出yn21。)h (八yn1 7n , 2 f(tnd，yn1) , f(tn，yn4)(2.18)如此重复下去，直到ynn；与yn、1)的差很小为止，此时可以认为yn? = yn+，其截断误差为h2还有二阶亚当姆斯显式公式，其可表示如下h _y(tn 1) ： yn 1 -3f(tn,yn) - f(tn4,yn)(2-19)yn1 =yn hB 4f (tn1,yn1B0 f (tn,yn)B-f(tnA1,yn*1)(2 - 20)表2-1到出了

49、各阶亚当姆斯法计算公式中的系数值，这是一个多步公式，为了计算yn+i,不仅要知道yn,还要知道yn-1,yn-2,，因此这些值都要保存，阶数越搞，要用到的值越多，存储量也就大了。多步法的缺点是不能自启动，即开始时要用单步法，然后才用多步法。名称B1B。B1B2B3一阶显式01000二阶显式03/2-1/200三阶显式023/12-16/125/120四阶显式055/24-59/2437/24-9/24:一阶隐式10000二阶隐式1/21/2000三阶隐式5/128/12-1/1200四阶隐式19/2419/24-5/241/240表2.1亚当姆斯系数表在隐式多步法中，一般先用显式多步法计算处值

50、，然后用隐式多步法作一次校正计算，这样使计算过程得到简化。例如，用一阶显式计算初值，即yn01=ynhf(tn,yn)(2-21)yn1=ynhf(tn,yn)f(tnh),y:0:(2-22)2以上计算公式与二阶龙格一库塔法公式相同。这种在多步法公式中先用显式公式估计一个值的方法叫预估一校正法，它的优点是在精度相同(四阶)，步长h相同时，预估一校正法，每步计算两次右函数，而龙格一库塔法要计算四次右函数，因此这种方法应用较为普遍。四阶亚当姆斯预估一校正法公式为：预估yF1=ynT-(55fn-59fnd37fn_2-9fnJ3)24校正YC+=yn+?(9炉虫+19fn5fn+fnj)24亚当

51、姆斯的显示公式与隐式公式的特点如下：(1)相同阶数的隐式公式的系数值比显式公式小(一阶例外)，因而隐式公式比显式公式精确得多.(2)隐式公式的稳定性亦比显式的好。下表给出了Adams法的稳定区。对于同样阶次，隐式的稳定域比显式的要大.随着阶数的增加，Adams法的稳定域逐步减小.阶数1234显式-216/11-3/10-8-8-6-3(3)从计算上看，显式比隐示计算量小。隐示公式需要计算fn+1,通常需要用Adams显式公式对它提供一个首次近似值yn+1。这种将显式公式和隐示公式联合起来使用以改进稳定性和精度的方法就是上面介绍的预估一校正法。亚当姆斯法每步只需要计算一次f函数值，因而计算量较小

52、，这是线性多步法的共同优点。但它需要知道多个出发值才能计算，而微分方程只能提供一个初值，显然多步法不象单步法那样可以自启动，它必须用其他方法先获得所求时刻以前多步的解。为了保证出发值的计算精度，一般采用比较高阶的单步法，以较小步长计算。此外，多步法改变步长不能象单步法那样方便。下面介绍一个使用亚当姆斯法进行模拟的程序，此程序可用于单输入单输出和多输入多输出线性系统的模拟，在模拟时面向状态方程。线性定常系统的状态空间表达式包括两个矩阵方程，X(t) = AX(t)Bu(T)(Y(t)=CX(t)+Du(t),X(t0)=X(0)第一个方程由n个一阶微分方程组成，称为状态方程，第二个方程是一个线性代数方程，称为输出方程，式中X为n*l的状态向量，u为m*1的控制向量。丫为l*l的输出向量，A为n*n阶状态矩阵，由控制对象的参数决定，B为n*m的控制矩阵，C为l*n的输出矩阵，D为l*m和直接传输矩阵。该程序用的计算方法为预估一校正线性多步法，首先用亚当姆斯显示公式计算预估值，再用亚当姆斯隐式公式计算校正值，其计算公式为预估：Xn1=Xn'(55Xn-59*37Xn_2-9Xnj3)24校正：XC书=Xn+；h(9XnF+19Xn-5XnJXn/)24输出：Yni=CX；.1DVAdambm函数流程图如下:计算预报值