模拟物理 第二章

1、模拟物理 第二章第1页，共81页，2022年，5月20日，6点35分，星期四主要内容简单方法：Euler方法多步法和隐式法Runge-Kutta方法例子：二维运动中的有序和混沌第2页，共81页，2022年，5月20日，6点35分，星期四最常遇到的任务许多物理定律通过微分方程表述。比如求微分方程的数值解是模拟物理系统时最常遇到的一项任务。第3页，共81页，2022年，5月20日，6点35分，星期四方程的表达形式常微分方程最一般的形式是一组M个耦合的一阶方程X是自变量，y是M个因变量，f有M个分量。高阶微分方程通过引入若干辅助函数可以写成这种一阶形式。第4页，共81页，2022年，5月20日，6点

2、35分，星期四牛顿方程-哈密顿方程例如，一个质量为m的粒子在力场F(z)作用下的一维运动由二阶方程描述若定义动量方程可以改为第5页，共81页，2022年，5月20日，6点35分，星期四只要详细讨论求解一阶方程组的方法就够了。只有一个因变量的情况，很容易推广到有多个因变量的情况。因此只讨论单个因变量的情况。本章重点讨论初值问题。即给定y(x=0)=y0，求y(x)第6页，共81页，2022年，5月20日，6点35分，星期四2.1 简单方法设我们想求解带有初始条件y(x=0)=y0的微分方程更具体地说，我们通常感兴趣的是某一特定x值（比如说x=1）上的y值。第7页，共81页，2022年，5月20日

3、，6点35分，星期四解法概述总的策略是把区间0,1分成N个等间隔的子区间，每个子区间的宽度h=1/N。然后找出一个递推公式，把yn同yn-1,yn-2,联系起来。其中yn是对y(xn=nh)的近似。这样一个递推关系将允许对这个微分方程进行从x=0到x=1的逐步积分。第8页，共81页，2022年，5月20日，6点35分，星期四Euler方法一个最简单的算法是Euler法。考虑xn点上的情况，并且把微商换成前向差分近似。递推关系第9页，共81页，2022年，5月20日，6点35分，星期四误差局部误差O(h2)整体误差：经过N步迭代，y(1)上的误差为NO(h2)O(h)这一误差线性地随h减小。为使

4、最后结果的不精确度减半，需要h减半，因而步数增加一倍。第10页，共81页，2022年，5月20日，6点35分，星期四例子考虑微分方程和初条件它的解析解为下面我们用Euler方法从x=0积分到x=3，步长由外界输入。第11页，共81页，2022年，5月20日，6点35分，星期四开始输入步长：h计算总步数 N=3/h初始条件 x=0,Y=1.0迭代一步: y=y+h*fnf(x,y)X=x+h循环计数器i=i+1输出误差iNN结束Yfnf(x, y)=-x*y第12页，共81页，2022年，5月20日，6点35分，星期四精度一般来说Euler法精度太低。这使我们不能通过采用较大的h来减少步数，从而

5、减少计算时间。在上面的例子中，当我们试图积分到更大的x值时，这一缺点就更明显。当x1/h是，y=0通常我们采用更高阶精度的算法第13页，共81页，2022年，5月20日，6点35分，星期四简单的高阶方法：Taylor级数方法一类简单的高阶方法可以由Taylor级数展开式导出对其中的导数可以进一步处理第14页，共81页，2022年，5月20日，6点35分，星期四局部误差O(h3)整体误差O(h2)，比Euler方法的精度高一阶。缺点：只有当f的解析形式已知，并且足够简单可以求微商时，才可以使用。第15页，共81页，2022年，5月20日，6点35分，星期四2.2 多步法和隐式法达到更高精度的另一

6、种方法是：使yn+1不仅同yn相联系而且同更早的点比如yn-1, yn-2,相联系，构造一个包含多步的递推关系。第16页，共81页，2022年，5月20日，6点35分，星期四多步法推导多步法公式：对微分方程做一步积分f(x,y)未知。取一个可以解析积分的近似。用xn和xn-1点上的f来做一个线性近似。带入积分二步法第17页，共81页，2022年，5月20日，6点35分，星期四更高阶的方法可以通过用更高次的多项式外插而得出。例如，若f用一个与fn,fn-1,fn-2,fn-3拟合的多项式来外插，就得到“四步法”第18页，共81页，2022年，5月20日，6点35分，星期四由于多步法的递推关系式包

7、含前面的好几步，单单关于y0的信息不能使它启动。因此必须通过别的方法，比如Euler方法、Taylor级数方法，或者下面讨论的Runge-Kutta方法，先得出前几个格点上的y值。第19页，共81页，2022年，5月20日，6点35分，星期四上述方法都是“显式的”。它意味着：yn+1是用已经知道的yn直接给出的（迭代）。隐式法：求解一个方程来决定yn+1。它提供了达到更高的精度的另一个手段。第20页，共81页，2022年，5月20日，6点35分，星期四我们设两个格点的中点xn+1/2=(n+1/2)h。考虑方程使用对称差分近似得到递推公式但是yn+1出现在两边。必须解方程。第21页，共81页，

8、2022年，5月20日，6点35分，星期四解方程可能很花时间。回顾：求根公式，搜索法，牛顿法，弦割法。如果f对于y是线性的，比方说f(x,y)=g(x)y，那么方程可以解出：第22页，共81页，2022年，5月20日，6点35分，星期四多步隐式法二步法：对方程 ，使用通过fn-1,fn,fn+1的二次多项式拟合f。解析积分得到隐式递推公式第23页，共81页，2022年，5月20日，6点35分，星期四用三次多项式内插可以推出对应的三步公式第24页，共81页，2022年，5月20日，6点35分，星期四隐式法真的要解方程吗？很少以通过解隐式方程的方式使用。用于“预估-校正”算法的基础。先以显式法得出

9、yn+1的预估值，再通过隐式法对它加以校正，给出更好的近似值。第25页，共81页，2022年，5月20日，6点35分，星期四2.3 Runge-Kutta方法对微分方程进行积分的算法是有些自由的。实际上，的确存在许多算法。每种算法都有其特点和优点。一种非常方便和广泛使用的方式是Runge-Kutta算法。它有不同阶的精度。第26页，共81页，2022年，5月20日，6点35分，星期四推导二阶公式对于 ，f用它在积分区间中点附近的Taylor级数展开式逼近。做解析积分其中hf来自零次项，一次项积分为0，误差来自二次项。第27页，共81页，2022年，5月20日，6点35分，星期四怎样处理yn+1

10、/2呢？用Euler方法产生它。误差为O(h2)二阶RK算法：第28页，共81页，2022年，5月20日，6点35分，星期四特点：它体现了把y的近似值带入隐式表达式右边的想法。优点：它同Taylor级数法或隐式法同样精确，但是并不对f加特殊的约束。不要求容易求微商或者f关于y是线性的。仅使用y在前面的一个格点上的值。第29页，共81页，2022年，5月20日，6点35分，星期四高阶算法用高阶算法计算xn到xn+1的积分，用高阶算法估计积分区间中格点上的y值可以推出高阶算法。第30页，共81页，2022年，5月20日，6点35分，星期四三阶算法第31页，共81页，2022年，5月20日，6点35

11、分，星期四四阶算法第32页，共81页，2022年，5月20日，6点35分，星期四经验发现， 四阶算法在精度和计算量之间给出最佳的折中。第33页，共81页，2022年，5月20日，6点35分，星期四使用RK方法的例子用4阶Runge-Kutta方法积分方程：求x=3时y的值。第34页，共81页，2022年，5月20日，6点35分，星期四第35页，共81页，2022年，5月20日，6点35分，星期四第36页，共81页，2022年，5月20日，6点35分，星期四2.5 二维运动中的有序和混沌在物理学中使用计算机带来的一个基本好处是：能够处理不能解析求解的系统。通常，数值结果同我们通过研究可解模型所发

12、展出的直观图象在定性上一致。而定量的数值结果又实际兴趣。在少数情况下，计算机结果否定了我们的直观。这时数值结果对正确理解这种现象有根本的重要性。第37页，共81页，2022年，5月20日，6点35分，星期四2.5 二维运动中的有序和混沌基本好处： 能够处理不能解析求解的系统。通常在少数情况下数值结果直观图象定性实际兴趣定量结果数值结果直观图象第38页，共81页，2022年，5月20日，6点35分，星期四否定经典图象令人惊奇的是，数值结果也否定了一些经典系统动力学的直观图象。牛顿以后，一切运转正常相对论指出高速运动下经典的直观图象不可行量子力学指出微观下经典的直观图象不可行计算机的应用指出经典情

13、形下，一些经典的直观图象不可行第39页，共81页，2022年，5月20日，6点35分，星期四做二维运动的粒子本例题中，我们将通过用数值积分计算一个做二维运动的粒子的轨道以此为例研究这种令人惊奇的行为。考虑一个单位质量的粒子。它在一个位势V中作二维运动。假定如果粒子的能量足够低，它将永远被约束在V中。第40页，共81页，2022年，5月20日，6点35分，星期四若与两个坐标(x, y)共轭的动量是(px, py)，则Hamilton量的形式为粒子的轨道由坐标和动量随时间的演化所确定。注意质量第41页，共81页，2022年，5月20日，6点35分，星期四Hamilton方程坐标和动量的随时间的演化

14、由四个耦合的一阶微分方程规定，即Hamilton方程第42页，共81页，2022年，5月20日，6点35分，星期四约束对于任何V，这些方程都能使能量E守恒，约束条件：因此，约束条件把轨道限制在嵌在四维相空间中的一个三维流形上。除此之外，关于系统的演化很难作出别的什么一般性判断。第43页，共81页，2022年，5月20日，6点35分，星期四可积系统能够对轨道作出进一步判断的一类重要的二维Hamilton量是可积的Hamilton量。对于这种情况，除了能量以外，还有第二个坐标和动量的函数是运动常数。于是轨道被约束在相空间的一个二维流形上。两种熟悉的可积系统是可分离变量的系统和中心位势。第44页，共

15、81页，2022年，5月20日，6点35分，星期四可积系统1：可分离变量可分离变量的情况下其中Vx, Vy是两个独立的函数。因此Hamilton量分离为两部分。两个方向的运动没有耦合。每一个方向上的Hamilton单独是一个运动常数。（等价于，Hx-Hx和E=Hx+Hy）H=Hx+HyHx=px2/2+Vx , Hy=py2/2+Vy第45页，共81页，2022年，5月20日，6点35分，星期四可积系统2：中心位势在中心位势的场合向心力导致角动量是第二个运动常数Hamilton量可以写为其中pr是r的共轭动量。第46页，共81页，2022年，5月20日，6点35分，星期四可积性的困难可积系统的

16、动力学是简单的。但是要找出这种简单性常常很不容易。没有一个普遍的解析方法可以判断，在任意一个位势中是否存在第二个运动常数，以及如果有的话如何求出它。第47页，共81页，2022年，5月20日，6点35分，星期四数值计算的固有局限数值计算看起来也不会改善情况。因为数值计算只给出给定初始条件下的轨迹。而这条轨迹即使在我们熟悉的情况下也可以是相当复杂的。比如Lissajous图形。（示例）第48页，共81页，2022年，5月20日，6点35分，星期四第49页，共81页，2022年，5月20日，6点35分，星期四第50页，共81页，2022年，5月20日，6点35分，星期四分析相空间通过对相空间的分析

17、，可以得出一个根据轨迹来检测可积性的办法。考虑一个可分离变量的位势的情况。因为两个坐标上的运动是独立的，在(x, px)平面和(y, py)平面上画出的轨迹看起来可能是第51页，共81页，2022年，5月20日，6点35分，星期四可积性的标志一个粒子在一个可分离变量的二维位势中的轨迹在(x, px)和(y, py)平面上的样子。这些闭合回路的存在是系统可积性的标志。第52页，共81页，2022年，5月20日，6点35分，星期四数值做法：截面根据轨迹可以得到一幅低维空间上的图。假设每当我们观察到一个坐标例如x通过零时，我们就在(y, py)平面上画出粒子的位置。如果x运动和y运动的周期是不可通约

18、的（即它们的比值是一个无理数），那么随着轨迹前进，这些观测值将逐渐描绘出完整的(y, py)回路。如果它们的周期是可通约的（即比值是有理数）那么将得到沿回路的一系列分立点。第53页，共81页，2022年，5月20日，6点35分，星期四利用截面方法，可以只从轨迹出发，来研究同任何给定的Hamilton量相联系的相空间的拓扑性质。第54页，共81页，2022年，5月20日，6点35分，星期四在一个可积的Hamilton量的情形下，相空间的一般的拓扑性质，可以通过考虑中心位势中的运动来说明。对于固定的能量值和角动量值，径向运动被限制在两个转折点rin和rout之间。这两个r值是方程的两个解： pr=

19、0 例如：卫星轨道第55页，共81页，2022年，5月20日，6点35分，星期四这两个半径在(x, y)平面中定义了一个圆环面。对于给定的r值，能量守恒允许的径向运动的动量只能是这两个动量在(x, y, pr)空间中确定了包含全部轨道的二维流形第56页，共81页，2022年，5月20日，6点35分，星期四第57页，共81页，2022年，5月20日，6点35分，星期四如果通过考虑x=0平面来画出一幅(y, py)截面图，我们将会得到两个闭合回路。画截面比画流形简单的多。x=0第58页，共81页，2022年，5月20日，6点35分，星期四可以证明，中心位势的相空间的环面拓扑结构，是一切可积系统所共

20、有的。对于给定的运动常数，轨迹所在的流形叫做一个“不变环面”。第59页，共81页，2022年，5月20日，6点35分，星期四环面的一般的截面图的样子第60页，共81页，2022年，5月20日，6点35分，星期四对可积性的扰动问题：在一个破坏掉系统的可积性的扰动下，一个可积系统的环面会发生什么变化？对于小扰动，椭圆型不动点周围的大部分环面发生轻微的畸变，但仍保持它们的拓扑结构（KAM定理）。但是相空间中邻接区域则变成混沌的，给出的截面图看起来像是随机撒的一些点子。在这些混沌区域中又嵌入别的椭圆型不动点和别的混沌区域，构成层次结构。第61页，共81页，2022年，5月20日，6点35分，星期四He

21、non-Heiles位势对可积性的重大偏离必须用数值方法研究。进行这一研究的一个方便的例子是这个位势原来是Henon和Heiles在研究恒星穿过星系的轨道时引入的。恒星的运动被认为限制在二维平面上。第62页，共81页，2022年，5月20日，6点35分，星期四第63页，共81页，2022年，5月20日，6点35分，星期四Hamiltion运动方程第64页，共81页，2022年，5月20日，6点35分，星期四有了运动方程，给定一组初始条件（x, y, px, py）就可以计算出一条轨道。这里使用4阶Runge-Kutta方法。第65页，共81页，2022年，5月20日，6点35分，星期四选择初条件考虑束缚轨道总能量 E1/6为了方便满足这个条件，我们直接选定E。第66页，共81页