连续系统的计算机模拟.docx
《连续系统的计算机模拟.docx》由会员分享，可在线阅读，更多相关《连续系统的计算机模拟.docx（52页珍藏版）》请在冰豆网上搜索。
连续系统的计算机模拟
第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的数量，它们的增加和减少的速度服从下列微分方程。
（2.1）
其中K1和K2是比例常数（一般而言这些比例常数会随温度和压力发生变化，但在模拟过程中，为了简化模型，一般不允许其变化，故一律视为常数）。
在给出常数K1和K2值以及A和B的数量（C=0）后，我们希望能确定有多少C物质产在出来，这种化学反应速率的决定在化学工业上是有意义的。
模拟该系统的一个直接的方法是在t＝0时开始，使t以Δt间隔增加。
假定化学量在Δt时间步长内不变，而只能在Δt结束的瞬间发生变化，这样在每个Δt结束时的A（或B或C）的数量就可以从Δt开始时的数值由下式求出
（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）-2K2c（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）-2K2c（Δ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）
YN
I<=N?
打印a,b,c
模拟完毕
图2.1化学反应例子模拟程序框图
模拟程序使用C语言编写，程序中的初值如下：
K1=0.008/g.minK2＝0.002/min,a（0）=100g,b（0）=50gC（0）=0,T=5mins,⊿t=0.1min,N=50
其程序清单如下：
#include
#include
floatk1,k2;
staticfloatA[53],B[53],C[53],delta,t;
voidstrut（int）;
main（）
{
inti;
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（"%2d",i）;
printf（"%10.2f,%10.2f,%10.2f,%10.2f\n",t,A[i],B[i],C[i]）;
strut（i）;
}
return;
}
voidstrut（inti）
{
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;
}
运行此程序，输出模拟结果如下：
10.00,100.00,50.00,0.00
20.10,96.00,46.00,8.00
30.20,92.47,42.47,15.06
40.30,89.33,39.33,21.34
50.40,86.52,36.52,26.95
60.50,84.00,34.00,32.00
70.60,81.72,31.72,36.55
80.70,79.66,29.66,40.69
90.80,77.77,27.77,44.45
100.90,76.05,26.05,47.89
111.00,74.48,24.48,51.04
121.10,73.03,23.03,53.94
131.20,71.70,21.70,56.61
141.30,70.46,20.46,59.07
151.40,69.32,19.32,61.36
161.50,68.26,18.26,63.48
171.60,67.28,17.28,65.44
181.70,66.36,16.36,67.28
191.80,65.51,15.51,68.99
201.90,64.71,14.71,70.59
212.00,63.96,13.96,72.08
222.10,63.26,13.26,73.48
232.20,62.60,12.60,74.79
242.30,61.99,11.99,76.03
252.40,61.41,11.41,77.18
262.50,60.86,10.86,78.27
272.60,60.35,10.35,79.30
282.70,59.87,9.87,80.27
292.80,59.41,9.41,81.18
302.90,58.98,8.98,82.04
313.00,58.57,8.57,82.86
323.10,58.19,8.19,83.63
333.20,57.82,7.82,84.36
343.30,57.48,7.48,85.05
353.40,57.15,7.15,85.70
363.50,56.84,6.84,86.32
373.60,56.55,6.55,86.91
383.70,56.27,6.27,87.46
393.80,56.00,6.00,87.99
403.90,55.75,5.75,88.50
414.00,55.51,5.51,88.97
424.10,55.29,5.29,89.43
434.20,55.07,5.07,89.86
444.30,54.86,4.86,90.27
454.40,54.67,4.67,90.66
464.50,54.48,4.48,91.03
474.60,54.31,4.31,91.39
484.70,54.14,4.14,91.73
494.80,53.98,3.98,92.05
504.90,53.82,3.82,92.35
515.00,53.68,3.68,92.65
525.10,53.54,3.54,92.93
以上例子的算法就是数值积分法中最简单的一种方法叫作欧拉法,从这个例子中我们可以看出使用计算机解题的过程，由于欧拉法本身的特点，决定了其计算精度差,现在几乎无人在实际工作中使用，但它导出简单，能说明构造数值解法一般计算公式的基本思想,模拟程序也清晰易懂，故人们乐意用它做为构造数值解法的入门例子。
其一般解法如下：
设给定微分方程
（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时的积分值为
（2.5）
将上式写成以下差分方程形式:
（2.6）
这就是欧拉公式。
它是一个递推的差分方程，任何一个新的数值解yn+1均是基于前一个数值解以及它的导数f（tn,yn）值求得的。
只要给定初始条件y0及步长h就可根据f（t0,y0）算出y1的值，再以y1算出y2，如此逐步算出y3,y4,…，直到满足所需计算的范围才停止计算。
欧拉法的基本思路是把连续的时间t分割成等间隔的Δt,在这些离散的时刻算得函数值，根据这些值在函数图上可得到一条折线，所以欧拉法又叫折线法，其特点是分析方法简单，计算量小，但计算精度低（后面将讨论欧拉法与其它方法的比较）。
下图为欧拉折线法的几何意义。
Yf（tn,yn）f（tn+1,yn+1）
欧拉折线
hy=∫f（t,y（t））dt+c
0tn-1tntn+1tn+2
图2.2欧拉法的几何意义
如果用梯形面积来代替一个小区间的曲线积分，就可克服以小矩形计算的缺点，提高精度，梯形法计算公式为
（2.7）
上式为隐式公式，因为公式右端含有yn+1,这是未知的待求量，故梯形法不能自行启动运算，而要依赖于其它算法的帮助，比如说用欧拉公式法求出初值，算出
的近似值
，然后将其带入微分方程，计算
的近似值
，最后利用梯形公式反复迭代。
如迭代一次后就认为求得了近似解，这就是改进的欧拉法，其公式为
预估
（2-8）
校正
（2-9）
上式中第一个为预估公式，第二个为校正公式。
通常这种方法称为预估矫正法。
在校正公式中计算了两点的斜率，再求其平均值，故计算量比欧拉法要大些。
2．3数值积分的几个基本概念
1.单步法与多步法
数值积分法都用递推公式求解,而递推公式是步进式的,当由前一时刻的数值yn就能计算出后一时刻的数值yn+1时,这种方法称为单步法,它是一种能自启动的算法,欧拉法是单步法.反之,如果求yn+1时要用到yn,yn-1,yn-2,…等多个值,则称为多步法,由于多步法在一开始时要用其它方法计算该时刻前面的函数值,所以它不能自启动,以上所讲的改进的欧拉法就是一个多步法的例子。
2显式与隐式
在递推公式中,计算yn+1时所用到的数据均已算出的计算公式称为显示公式.相反,在算式中隐含有末知量yn+1的则称为隐含公式.这需用另一个显示公式估计一个值,再用隐式公式进行迭代运算,这就是预估----校正法,这种方法不能自启动.
用单步法与显示公式在计算上比用多步法和隐式公式方便.但有时要考虑精度与稳定性时,则必须采用隐式公式运算.
3截断误差
一般使用台劳级数来分析数值积分公式的精度.为简化分析,假定前一步得的结果yn是准确的,则用台劳级数求得tn+1处的精确解为
（2.10）
与前面的递推公式比较,欧拉法是从以上精确解中只取前两项之和来近似计算yn+1的,由这种方法单独一步引进的附加误差通常称着局部截断误差，它是该方法给出的值与微分方程解之间的差值，故又称为局部离散误差。
欧拉法的局部截断误差为
不同的数值解法，其局部截断误差也不同。
一般若截断误差为
，则方法称为r阶的，所以方法的阶数可以作为衡量方法精确度的一个重要标志。
欧拉法是一阶精度的数值解法，而改进的欧拉法（梯形法）相当于取台劳级数的前三项，故其解为二阶精度。
4舍入误差
由于计算机的字长有限,数字不能表示得完全精确,在计算过程中不可避免地会产生舍入误差.舍入误差与步长h成反比,如计算步长小,计算次数多则舍入误差大.产生舍入误差的因素较多,除了与计算机字长有关外,还与计算机所使用的数字系统,数的运算次序及计算f（t,y）所用的程序的精确度等因素有关.
5数值解的稳定性
以上对连续系统的模拟用到了差分方程,把初始值代入递推公式进行迭代运算,如果原系统是稳定的,由数值积分法求得的解也应该是稳定的.但是,在计算机逐次计算时,初始数据的误差及计算过程中的舍入误差对后面的计算结果会产生影响.而且计算步长选择得不合理时,有可能使模拟结果不稳定.以下我们简要地讨论这一问题.
差分方程的解与微分方程的解类似,均可分为特解与通解两部分.与稳定性有关的为通解部分,它取决于该差分方程的特征根是否满足稳定性要求。
例如,为了考虑欧拉法的稳定性，可用检验方程y=λy。
其中λ为方程的特征根。
对此方程有
要使该微分方程稳定，须使下式成立
|1+λh|<1
有此可得
|λh|<2或h<2T
要保证用欧拉法进行的计算是稳定的,积分步长h必须小于系统时间常数的两倍.所以积分步长的选择就要慎重,不能随意取任何值。
2．4龙格--贝塔（Runge--Kutta）法
仍以（2.3）方程为例.已知y（t0）=y0,假设我t0开始以h增长,t1=t0+h，t1时刻为
在
附近展开成台老级数，保留
项，则有
（2-11）
上式括号内下标为0，表示括号中的函数用t=t0,y=y0代入，以下均同。
假设上式的解可写成如下形式
其中,K1=f（t0,y0）
对K2式右端得函数在t=t0,y=y0处展开成台劳级数，保留h项,可得到:
将K1与K2代入（2-12）式得:
将上式与（2-12）式比较，可得系数a1,b1,c2,b2方程如下：
以上三个方程,四个未知数,因此有无穷个解,若限定b1=b2则可得一个解:
将它们代入（2-12）式可得出一组计算公式
其中K1=f（t0,y0）,k2=f（t0+h,y0+k1h）,写成一般的递推形式如下：
式中：
四阶龙格—库塔法公式精度较高，故其应用较为普遍。
下面再给出几组系数值
读者可以写出其相应的四阶龙格—库塔法公式。
龙格—库塔法的特点如下：
1.龙格－库塔法有多组a1,c2b1b2值，故可有多种龙格－库塔公式，常使用的有：
所有龙格－库塔公式都有以下几个特点：
（1）在计算yn＋1时只用到yn，而不直接用yn－1，yn－2等项，也就是计算后一步时
只用到前一步的计算结果，故称为单步法，单步法可以自启动，且使用的存储量也小。
（2）在计算过程中，步长h可以改变，这要由计算精度来定，但是在一步中，为计算
若干个系数Ki（称成为龙格－库塔系数），必须使用同一个步长h。
（3）.不论是几阶龙格－库塔公式，计算式中均为两部分组成。
第一部分为上一步的结果yn，第二部分为h=tn--tn－1中对f（t,y）的积分。
它是步长h乘上各点斜率的加权平均。
例如上面计算的四阶公式中，取四点的斜率：
K1、K2、K3和K4,然后对K2和K3取两份，K1和K4各取一份，进行加权平均。
3.龙格－库塔法的精度取决于步长h的大小及求解的方法。
实践证明：
为达到同样的精度，四阶方法的步长h可以比二阶方法的h大10倍，而四阶方法的每步计算量仅比二阶方法大一倍，所以总的计算量仍比二阶方法小，使用的CPU机时也少。
一般工程中四阶龙格－库塔公式已能达到要求的精度，故不再使用更高阶的公式。
下表示出了f的计算次数与精度阶数的关系。
每一步计算f的次数
2
3
4
5
6
7
n=>8
精度阶数
2
3
4
4
5
6
n-2
可见精度的阶数与计算函数值f的次数之间的关系并非等量增加。
其精度最低。
如果将h取得很小，能否用欧拉公式得到很高的精度呢？
从理论上讲是可以的，但是实际上由于计算机字长有限，在计算中有舍入误差。
步长h越小，计算的步数越多，舍入误差的积累会越来越大，故用欧拉法时不能用很小的h。
2.5四阶龙格—库塔法模拟程序及应用
目前已有多种四阶龙格—库塔法模拟（仿真）软件包推出，虽然子程序稍有不同，但总的结构还是相同的。
对于连续系统的龙格—库塔法的计算机程序，其大致结构如下图所示：
主函数
输入及运行输出
初始化
微分方程数值积分显示打印绘图
右函数
图2.3龙格—库塔法程序简要框图
图2.3中每一个程序模块的功能为：
1.主函数：
主函数实现对整个模拟系统的控制,这是通过调用各个子程序实现的。
2.输入及初始化函数：
主要对系统参数输入初值，例如模拟时间，积分步长，方程阶数等。
这可以从键盘输入，也可从数据文件输入。
3．运行模块:
在运行子程序中,将反复调用数值积分和方程右函数模块,进行迭代计算，可以每计算一步便显示一次结果，也可以计算数步显示一次结果，屏幕的显示使用户可以随时监视系统运行的进程，以便人工干预。
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（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）+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"math.h"
#defineNORDER3
#defineNPARAM2
floaty[NORDER],d[NORDER],p[NPARAM],t;
voidoutput（void）
{intj;
printf（"%9s","time"）;
for（j=0;jprintf（"y[%d]",j）;
putchar（'\n'）;
}
voiddifq（void）
{d[0]=-p[0]*y[0]*y[1];
d[1]=p[0]*y[0]*y[1]-p[1]*y[1];
d[2]=p[1]*y[1];
}
voidrun（floatdt）
{inti;
floatrk[3][NORDER],ys[NORDER];
difq（）;
for（i=0;i{rk[0][i]=d[i]*0.5*dt;
ys[i]=y[i];
y[i]=ys[i]+rk[0][i];
}
t+=0.5*dt;
difq（）;
for（i=0;i{rk[1][i]=d[i]*0.5*dt;
y[i]=ys[i]+rk[1][i];
}
difq（）;
for（i=0;i{rk[2][i]=d[i]*dt;
y[i]=ys[i]+rk[2][i];
}
t+=0.5*dt;
difq（）;
for（i=0;iy[i]=ys[i]+1.0/6*（2*（rk[0][i]+rk[2][i]）+4*rk[1][i]+d[i]*dt）;
}
voidmain（void）
{charc;
do
{floatdt,tmax,co,tnext,tol;
intj;
printf（"Thisprogramwillfindtherootof:
\n"）;
printf（"{y[0]'（t）=-p[0]*y[0]*y[1]\n"）;
printf（"&&{y[1]'（t）=p[0]*y[0]*y[1]-p[1]*y[1]\n"）;
printf（"&&{y[2]'（t）=p[1]*y[1]\n"）;
for（j=0;j{printf（"Now,pleaseinputthefirsty[%d]:
",j）;
scan
升级会员 