第3章连续系统的数字仿真通用算法

第3章连续系统的数字仿真通用算法欧拉(Euler)法数值积分的几个概念龙格-库塔(Runge-Kutta)法建模实例一、欧拉(Euler)法[例3-1]当两种物质A和B放到一起产生化学反应时，产生第三种物质C，一般1克A与1克B结合产生2克的C物质，形成C的速率与A和B的数量乘积成正比，同样C也可分解为A和B，C的分解速率正比于C的数量，即在任何时刻，如果a，b，c分别是化学物质A，B，C的数量，它们的增加和减少的速度服从下列微分方程。（3-1）要求：在给出常数K1和K2值以及A和B的数量(C=0)后，我们希望能确定有多少C物质产生出来，这种化学反应速率的决定在化学工业上是有意义的。模拟该系统的一个直接的方法是在t=0时开始，使t以

t间隔增加。假定化学量在

t时间步长内不变，而只能在

t结束的瞬间发生变，这样在每个

t结束时的A(或B或C)的数量就可以从

t开始时的数值求出：（3-2）假定模拟周期为T，可将T分成N个小的时间步长t，即：T=N·t在时间为0时，我们知道a(0),b(0)和c(0)的实始值，从这些初始值及常数K1和K2值出发，就可以计算出t时间内的化学量的变化值。输入K1,K2,a(0),b(0),c(0)T,t和NI=0I=I+1计算保存a(I,t),b(I,t),c(I,t)I<=N?打印a,b,c模拟结束YN作业：已知k1=0.008/(g.min),k2=0.002/min,a(0)=100g,b(0=50g,c(0)=0,T=5min,t=0.1min,N=50,试模拟该问题，并输出计算结果。欧拉法的一般解法设给定微分方程（3-3）在区间[tn,tn+1]上求积分，得（3-4）（3-5）（3-6）如果用梯形面积代替则较好地解决欧拉缺点（3-7）（3-8）（3-9）预估校正二、数值积分的几个概念（一）单步法与多步法单步法：当由前一时刻的数值yn就能计算出后一时刻的数值yn+1时，这种方法称为单步法，它是一种能自启动的算法，欧拉法是单步法。多步法：如果求yn+1时要用到yn,yn-1,yn-2,…等多个值，则称为多步法，多步法在一开始时要用其他方法计算该时刻前面的函数值，所以它个能自启动的。（二）显式与隐式在递推公式中，计算yn+1时所用到的数据均巳算出的计算公式称为显示公式。相反，在算式中隐含有未知量yn+1的则称为隐含公式，这需用另—个显示公式估计一个值，再用隐式公式进行迭代运算，这就是预估-校正法，这种方法不能自启动。（三）截断误差一般位用台劳级数来分析数值积分公式的精度。为简化分析，假定前一步所得的结果yn是准确的，则用台劳级数求得tn+1处的精确解为：（3-10）（三）截断误差与前面的递推公式比较、欧拉法是从以上精确解中只取前两项之和来近似计算儿yn+1的，由这种方法单独一步引进的附加误差通常称着局部截断误差，它是该方法给出的值与微分方程解之间的差值，故又称为局部离散误差。欧拉法的局部截断误差为：y(tn+1)-y(tn)=0(h2)不同的数值解法，其局部截断误差也不同。一般若截断误差为0(hr+1)，则方法称为r阶的所以方法的阶数可以作为衡量方法精确度的一个重要标志。欧拉法是一阶精度的数值解法，而改进的欧拉法(梯形法)相当于取台劳级数的前三项，故其解为二阶精度。（四）合入误差由于计算机的字长有限，数字不能表示得完全精确，在计算过程中不可避免地会产生舍入误差。舍入误差与步长h成反比，如计算步长小，计算次数多则舍入误差大。产生舍入误差的因素较多，除了与计算机字长有关外，还与计算机所使用的数字系统，数的运算次序及计算f(t,y)所用的程序的精确度等因素有关。（五）数值解的稳定性三、龙格-库塔(Runge-Kutta)法仍以(3-3)方程为例，已知y(t0)＝y0，假设从t0开始以h增长，t1＝t0+h，t1时刻为y1=y(t0+h)，在t0附近展开成台劳级数，保留h2项，则有：（3-11）假设上式的解可写成如下形式:（3-12）对K2式右端的函数在t=t0,y=y0处展开成台劳级数,保留h项,可得到:将K1与K2代入(3-12)式,得:将上式与(3-12)式比较,可得系数a1,b1,c2,b2方程如下:若b1=b2将它们代入(3-12)可得出一组计算公式:写成一般的递推形式如下:（3-13）式中只取h,h2两项,而将h3以上的高阶项略去,所以此递推公式的截断误差正比于h的三次方,计算过程中只取h和h2两项,故此方法被称为二阶龙格—库塔法.一般使用的为四阶龙格—库塔公式,在展开台劳级数时保留h,h2,h3和h4项,其截断误差为h5,四阶龙格—库塔公式为:下面给出几组系数值（3-14）龙格—库塔法的特点(1)龙格—库塔法有多组a1,c2,b1,b2值,故可有多种龙格—库塔公式,常使用的有:（3-15）（3-16）二阶四阶(2)所有龙格—库塔公式都有以下特点①在计算yn+1时只用到yn，而不直接用yn-1，yn-2等项，也就是计算后—步时只用到前一步的计算结果，故称为单步法，单步法可以自启动，且使用的存储量也小。②在计算过程中，步氏h可以改变，这要由计算精度来定，但是在每一步中，为计算若干个系数K(称成为龙格一库塔系数)，必须使用同一个步长h。③不论是几阶龙格--库塔公式，计算式中均为两部分组成。第一部分为上一步的结果yn，第二部分为h=tn-tn-1中对f(t,y)的积分。它是步h乘上各点斜率的加权平均。例如上面计算的四阶公式中，取四点的斜率：K1，K2，K3和K4，然后对K2和K3取两份，K1和K4各取一份，进行加权平均。(3)龙格--库塔法的精度取决于步长h的大小及求解的方法。实践证明：为达到同样的精度，四阶方法的步长h可以比二阶方法的h大10倍，而四阶方法的每步计算量仅比二阶方法大一倍，所以总的计算星仍比二阶方法小，使用的CPU机时也少。一般工程中四阶龙格--库塔公式已能达到要求的精度，故不再使用更高阶的公式。表3-1示出了f的计算次数与精度阶数的关系。表3-1计算f的次数与精度阶数的关系第一步计算f的次数234567

8精度阶数234456n-2(4)若在展开成台劳级数时，只取h这一项，而将h2以及h2以上的高次项略去，可得yn+1=yn+hf(tn,yn)这是欧拉公式，因此欧拉公式也可看成时一阶龙格--库塔公式，它的截断误差正比h2，其精度最低。如果将h取得很小，能否用欧拉公式得到很高的精度呢?从理论是讲是可以的，但是实际上由于计算机字长有限，在计算中有合入误关。步长h越小，计算的步数越多，合入误差的积累会越来越大，故用欧拉法时不能用很小的h。四、建模实例传染病的传播问题在发生自然灾害的年份内，常常伴有各种瘟疫或传染病流行，被传染的人数与哪些因素有关?如何预报传染病高潮的