实验六 计算机数值模拟实验

实验六 计算机数值模拟实验

计算机数值模拟方法是从基本的物理定律出发，用离散化变量描述物理体系的状态， 然后利用电子计算机计算这些离散变量在基本物理定律制约下的演变，从而体现物理过 程的规律．

计算机数值模拟实验是在计算机中进行的实验．虽然它不能替代真实的物理实验，但 确实是一种极其重要的实验方法．它是通过大量“个例”来研究特定的物理过程，能够反复 进行，方便地控制和调整参数，在理论研究和实验研究之间搭起了一座“桥梁”．数值模拟 可以研究一些非常复杂的过程，而理论研究必须作出许多简化假设才能处理这些过程，简 化则意味着可能丢失许多重要的因素，这就使得数值模拟可以更全面地了解一个物理过 程，而且还可能发现新的物理现象．另一方面，数值模拟也能够为实验观测方案提供理论 的支持，对大型实验装置进行评估，对实验条件或参数进行优化选择，以避免造成极大的 经济损失和人力浪费．随着计算机性能的高速发展，数值模拟在各门学科的研究中应用将 更加广泛，起到越来越重要的作用．

本实验选择一个非线性动力学系统中的混沌吸引子作为研究实例．这类问题用常规 的理论方法和实验方法是不容易了解其运动规律的，但通过计算机模拟实验我们便可得 到具体的物理图像．实验中要求读者通过在IBMPC 或兼容机上的计算机编程和实验，掌 握数值模拟的基本方法和步骤．

一、实验原理

1．数值模拟的基本方法． 由于各门学科研究的对象都具有自身的特点，使得数值模 拟的方法在不同学科中有其不同的特征，但任何数值模拟都需要求解描述相应物理过程 的数学方程．这些数学方程的数值求解方法有其共性，因而各门学科中数值模拟方法有共 同之处，一般都涉及到如下几个步骤：

(1)建立物理模型．对任何物理过程的数值模拟都首先要建立模型，所建立的模型的 合理性在很大程度上决定了模拟结果是否可靠．建立模型包括如下步骤：

a) 找出决定所研究的物理过程的主要因素；

b) 导出适当的数学方程；

c) 给出切合物理实际的边值条件和初始条件．

(2)方程和初值、边值条件的离散化：

a) 选择合适的数值方法，常用的有差分法、有限元法和边界元法；

b) 将计算区域划分为离散网格点，网格点(时间、空间) 可以是规则的，也可以是不规则 则的，取决于计算区域的几何形状是否规则，网格点多少(决定于时间和空间步长的大的选择需要综合考虑计算时间的长短、计算机内存的大小等因素；

c) 将方程和初值、边值条件化为网格点上的代数方程(组) ．

(3)选择适当的代数方程组求解方法．

(4)在计算机上实现数值求解：

a) 设计流程图；

b) 编写计算机程序；

c) 调试程序，检查程序是否有语法错误、·数学公式的程序语言表达是否正确，根据计 算结果检查算法的计算精度，以及是否有数值的不稳定性存在．

(5)计算结果的诊断．诊断是将数值模拟结果以一定的形式(通常是图形) 表达出来， 这是调试程序和输出结果的重要手段．

2．利用数值模拟的实例． 作为一个实验的具体例子，我们研究美国气象学家罗伦兹 (E．N ．Lorenz) 于1963年在大气科学杂志上提出的第一个表现奇异吸引子的动力学系统．该

混沌系统描述了从水桶的底部加热时，桶内液体的运动情况．加热时，底部的液体越来越热，并开始逐渐上升，产生对流．当提供足够的热量并保持不变时，对流便会产生不规则的运动和湍流．

(1)建立模型．该混沌系统模型可以用下列三个微分方程描述：

dx =ay (t ) -az (t ); (6．1) dt

dy =r (t ) x (t ) -y (t ) -x (t ) z (t ); (6．2) dt

dz =x (t ) y (t ) +bz (t )) (6．3) dt

其中x 正比于对流运动的速度，y 正比于水平方向温度的变化，z 正比于竖直方向温度的变化，系数a 通常取值为10.0，b 通常取值为8／3，r 正比于水桶底部和水桶顶部之间的温度变化，是该动力学系统模型中重要的参数，在实验中可以采用常数，或采用周期瑞利数 r(t)＝r 0+r1coswt ．关于该模型的建立和参数的详细讨论可参阅有关参考文献．

该方程组的初始条件为

x(t=0)=x0 , y(t=0)=y0 , z(t=0)=z0.

此问题没有边值条件．

(2)方程的离散化，方程(6．1) 、(6．1) 、(6．3) 与初始条件一起构成一个一阶微分方程组，可以采用四阶龙格—库塔(Runge—Kutta) 法求解．龙格—库塔方法计算公式为：如果微分方程为dx ／dt ＝f(x，y ，z ，t) ，则：

xinc1=f （x ，y ，z ，t ）*dt

xinc2=f(x+xinc1/2，y +yinc1/2，z +zinc1/2)*dt ，

xinc3=f(x+xinc2/2，y +yinc2/2，z +zinc2/2)*dt ，

xinc4=f(x+xinc3，y+yinc3，z +zinc3，t+dt )*dt ，

x n+1=x n +(xinc1+2*xinc2+2*xinc3+xinc4)/6

其中x n 为第n 个迭代点，X n+1第n+1个迭代点，dt 为时间步长．根据上式，从初始值 x(t＝0) ＝x 0，y(t＝0) ＝y 0，z(t＝0) ＝z 0开始就可以计算出以后各个时间的x ，y 和z 的值．

图6．1 实验方框图

(3)实验的演示程序．附录给出了罗伦兹吸引子的C 语言演示程序．程序中利用了四 阶龙格—库塔算法，在TurboC++3．0中编译通过．

程序流程如图1所示，其主要步骤如下：

a) 要求用户输入xyz 的初值，r 。，r 1，以及DisplayAfter ．其中DisplayAfter 指定了在开始显示x--y 曲线以前，计算机所进行的迭代运算的次数．该参数的加入，是便于用户观察到程序进行了DisplayAfter 个迭代运算后的输出结果．在演示程序中，固定的参数和变量为： 参数a ，固定为a =10；

参数b ，固定为b ＝8．0／3．0；

时间步进dt ，固定为dt ＝0.001；

r(t)中的角频率ω，固定为ω＝7.62．

b) 迭代次数置0，时间置0．

c) 判断是否有键按下，如有则停止迭代，否则执行以下几步：

d) 迭代次数是否大于DisplayAfter ，如是则显示该次迭代结果，否则迭代结果不显示 出来(但仍然进行下一步迭代计算) ．迭代结果的显示可以根据需要显示x-t 图、y-t 图、z-t 图、x-y 图、x-z 图、y-z 图．在演示程序仅显示x-y 图．

e) 进行x ，y ，z ，的迭代计算。在程序中，先计算该次的xincl ，xinc2，xinc3，xinc4，

yincl ，yinc2，yinc3，yinc4，zincl ，zinc2，zinc3，zinc4后，再计算新的x ，y ，z 迭代结果．如果计算完xincl ，xinc2，xinc3，xinc4后立即迭代出新的x ，再用该新的x 和原来的y,z 去计

算yincl ，yinc2，yinc3，yinc4，再迭代出的y 会有错误．这一点在编程中需要引起重视． f) 迭代次数加1，时间增加一个步进dt ．

g) 执行第3步．

(4)罗伦兹混沌系统的实验观察．

a) 初值对混沌系统的影响(蝴蝶效应) ．混沌系统是一个非线性系统，初值对系统的敏感性是混沌运动的一个基本特征．x ，y, z的初值对混沌动力学系统有很大影响．

我们以天气预报系统为例进行说明．在气象系统模型中，天气预报所用的信息可以包 括从各个气象站获得的风速、温度、气压等测量参数，将这些信息输入计算机模型中观察天气预报．每次从气象站获得的信息不能保证都是准确的，例如在某次测量风速时可能在距离不远的公路上正好有一辆大卡车通过，或者近处有一只鸟儿的翅膀在拍打，甚至是一只蝴蝶的翅膀在振翅．

在混沌现象发现以前，人们通常认为这些测量的微小误差对天气预报的影响只是短时的，它对长时间的预报不会有影响．而图示模拟结果表明，即使是一个蝴蝶的拍打都会影响到三个月后的天气．这就是蝴蝶效应的意义．

不过，长时间的天气预报并没有因此而消失．随着计算机技术的发展，气象学家有了办法对付蝴蝶效应：先用一些初始条件进行仿真，再使这些初始条件有细小的变化，如果天气

图是完全改变的，则表明初始的气象系统是混沌的，不能用该模型进行长时间间隔的天气预测．但是，如果初始条件的变化并没有影响到天气图，则表明初始的气象系统不是混沌的，此时就可以进行预测了．

在本实验中，可以先设定初值为x ＝y ＝z ＝1，r 0＝28(r 1＝0保证了系统内部处于相当不稳定的对流状态) ，r 1＝0．观察程序刚开始迭代的几百次的输出和迭代了3000次后的输出(DisplayAfter＝1和DisplayAfter ＝30000) ．

然后，在其他条件不变的情况下，如果将z 的初值从z=1改变为z=1.001，并重复实验。尽管z 的初值变化只有一千分之一，但是观察在DisplayAfter ＝30000时，初值变化前后输出有了明显的变化．

从上述初值对系统的敏感性我们知道，如果在实际的应用模型中采用了微分方程组，又有可能出现混沌的话，那么在测量中任何初值的误差都将对系统产生很大的影响，变得不能被忽略．

b) 罗伦兹混沌吸引子．混沌是一种非周期的动力学过

程，看似无序，杂乱无章，但却隐含着丰富的内容，如混

沌吸引子、分支、窗 口等等．它是一种无序中的有序，

决不仅仅是一个无从控制的随机过程．

混沌吸引子是相空间的某部分，从它附近出发的任何

点都逐渐趋近于它．在(x,y) 平面上，我们可以看到形如

肾脏的两叶的罗伦兹混沌吸引子图(图6．2) ．其中的点时

而转到左页，时而转到右页．所有点的轨迹都螺旋趋近两

叶中心不会离开闭合的曲面．

c) 倍周期运动研究．尽管罗伦兹方程所反映

的是一个非线性的混沌系统．但是，进一步的实验表明，在某些条件下，仍然可以有周期性的轨迹出现．例如，在采用周期瑞利数r(t)＝r 0+r1coswt 时，无论系统的初始状态(x,y,z 的初值) 如何，如果r 1从0增加至5.0，在r 0=26.5时，如果我们观察系统经过较长时间迭代后的状态(DisplayAfter=50000)，可以发现此时系统出现了稳态的单周期(图6．1．3a) ，在r 0增大的过程中，可以获得稳态的双周期(图6．1．3b) 、稳态的4周期(图6．1．3c) 、稳态的64周期(图6．1．3d) 等等。

又例如，当对r 1=10时，系统对不同的x,y,z 初值的敏感性很快消失了．系统此时处于稳定对流状态．如果r 1很大，系统将处于周期状态而不是混沌状态．

二、实验装置

1．计算机硬件． IBMPC 或兼容机一台．计算机性能的好坏主要影响模拟计算的速 度．建议采用80486或速度更快的计算机．打印机一台，需要支持图形打印方式．

2．计算机软件． DOS5．0以上操作系统．高级编程语言一种，例如BASIC ，PAS — CAL ，C ，FORTRAN 等．

三、实验内容

1．自选一种编程语言，对Lorenz 吸引子进行模拟计算．

2．编程获得x-y 图的罗伦兹混沌吸引子图．

3．验证混沌系统对初值的敏感性z=1和z=0.001时获得x-y 图、z-t 图，

4．验证混沌系统的倍周期运动．

5．将你的实验结果和附录程序的结果进行比较．

四、思考和讨论

1．试说明数值模拟方法的特点，它与理论研究、实验研究有什么关系?

2．从混沌系统的基本特征出发，联系天气预报系统，说明蝴蝶效应的意义．

3．将时间步进dt 由0．001改为0．01，观察实验结果并进行比较，如果有差异试分析 其原因．

4．你如何理解混沌现象，倍周期运动和初值的敏感性。

附录：Lorenz 吸引子的C 语言演示程序