模拟物理 第二章精选文档.ppt

《模拟物理 第二章精选文档.ppt》由会员分享，可在线阅读，更多相关《模拟物理 第二章精选文档.ppt（81页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、模拟物理 第二章本讲稿第一页，共八十一页主要内容1.简单方法：Euler方法2.多步法和隐式法3.Runge-Kutta方法4.例子：二维运动中的有序和混沌本讲稿第二页，共八十一页最常遇到的任务许多物理定律通过微分方程表述。比如求微分方程的数值解求微分方程的数值解是模拟物理系统时最常遇到的一项任务。本讲稿第三页，共八十一页方程的表达形式常微分方程最一般的形式是一组M个耦合的一阶方程X是自变量，y是M个因变量，f有M个分量。高阶微分方程高阶微分方程通过引入若干辅助函数可以写成这种一阶形式。本讲稿第四页，共八十一页牛顿方程-哈密顿方程例如，一个质量为m的粒子在力场F(z)作用下的一维运动由二阶方程

2、描述若定义动量方程可以改为本讲稿第五页，共八十一页1.只要详细讨论求解一阶方程组的方法就够了。2.只有一个因变量的情况，很容易推广到有多个因变量的情况。因此只讨论单个因变量的情况。3.本章重点讨论初值问题。即给定y(x=0)=y0，求y(x)本讲稿第六页，共八十一页2.1 简单方法设我们想求解带有初始条件y(x=0)=y0的微分方程更具体地说，我们通常感兴趣的是某一特定x值（比如说x=1）上的y值。本讲稿第七页，共八十一页解法概述总的策略是把区间0,1分成N个等间隔的子区间，每个子区间的宽度h=1/N。然后找出一个递推公式，把yn同yn-1,yn-2,联系起来。其中yn是对y(xn=nh)的近

3、似。这样一个递推关系将允许对这个微分方程进行从x=0到x=1的逐步积分。本讲稿第八页，共八十一页Euler方法方法一个最简单的算法是Euler法。考虑xn点上的情况，并且把微商换成前向差分把微商换成前向差分近似近似。递推关系本讲稿第九页，共八十一页误差局部误差O(h2)整体误差：经过N步迭代，y(1)上的误差为NO(h2)O(h)这一误差线性地随h减小。为使最后结果的不精确度减半，需要h减半，因而步数增加一倍。本讲稿第十页，共八十一页例子考虑微分方程和初条件它的解析解为下面我们用Euler方法从x=0积分到x=3，步长由外界输入。本讲稿第十一页，共八十一页开始输入步长：h计算总步数 N=3/h

4、初始条件 x=0,Y=1.0迭代一步:y=y+h*fnf(x,y)X=x+h循环计数器i=i+1输出误差iNN结束Yfnf(x,y)=-x*y本讲稿第十二页，共八十一页精度一般来说Euler法精度太低。这使我们不能通过采用较大的h来减少步数，从而减少计算时间。在上面的例子中，当我们试图积分到更大的x值时，这一缺点就更明显。当x1/h是，y=0通常我们采用更高阶精度的算法本讲稿第十三页，共八十一页简单的高阶方法：Taylor级数方法一类简单的高阶方法可以由Taylor级数展开式导出对其中的导数可以进一步处理本讲稿第十四页，共八十一页局部误差O(h3)整体误差O(h2)，比Euler方法的精度高一

5、阶。缺点：只有当f的解析形式已知，并且足够简单可以求微商时，才可以使用。本讲稿第十五页，共八十一页2.2 多步法和隐式法达到更高精度的另一种方法是：使yn+1不仅同yn相联系而且同更早的点比如yn-1,yn-2,相联系，构造一个包含多步的递推关系。本讲稿第十六页，共八十一页多步法推导多步法公式：对微分方程做一步积分f(x,y)未知。取一个可以解析积分的近似。用xn和xn-1点上的f来做一个线性近似。带入积分二步法本讲稿第十七页，共八十一页更高阶的方法可以通过用更高次的多项式外插而得出。例如，若f用一个与fn,fn-1,fn-2,fn-3拟合的多项式来外插，就得到“四步法”本讲稿第十八页，共八十

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

7、：求根公式，搜索法，牛顿法，弦割法。如果f对于y是线性的，比方说f(x,y)=g(x)y，那么方程可以解出：本讲稿第二十二页，共八十一页多步隐式法二步法：对方程 ，使用通过fn-1,fn,fn+1的二次多项式拟合f。解析积分得到隐式递推公式本讲稿第二十三页，共八十一页用三次多项式内插可以推出对应的三步公式本讲稿第二十四页，共八十一页隐式法真的要解方程吗？很少以通过解隐式方程的方式使用。用于“预估-校正”算法的基础。先以显式法得出yn+1的预估值，再通过隐式法对它加以校正，给出更好的近似值。本讲稿第二十五页，共八十一页2.3 Runge-Kutta方法对微分方程进行积分的算法是有些自由的。实际上

8、，的确存在许多算法。每种算法都有其特点和优点。一种非常方便和广泛使用的方式是一种非常方便和广泛使用的方式是Runge-Kutta算法。算法。它有不同阶的精度。本讲稿第二十六页，共八十一页推导二阶公式对于 ，f用它在积分区间中点附近的Taylor级数展开式逼近。做解析积分其中hf来自零次项，一次项积分为0，误差来自二次项。本讲稿第二十七页，共八十一页怎样处理yn+1/2呢？用Euler方法产生它。误差为O(h2)二阶RK算法：本讲稿第二十八页，共八十一页特点：它体现了把y的近似值带入隐式表达式右边的想法。优点：它同Taylor级数法或隐式法同样精确，但是并不对f加特殊的约束。不要求容易求微商或者

9、f关于y是线性的。仅使用y在前面的一个格点上的值。本讲稿第二十九页，共八十一页高阶算法用高阶算法计算xn到xn+1的积分，用高阶算法估计积分区间中格点上的y值可以推出高阶算法。本讲稿第三十页，共八十一页三阶算法本讲稿第三十一页，共八十一页四阶算法本讲稿第三十二页，共八十一页经验发现，四阶算法在精度和计算量之间给出最佳的折中。本讲稿第三十三页，共八十一页使用RK方法的例子用4阶Runge-Kutta方法积分方程：求x=3时y的值。本讲稿第三十四页，共八十一页本讲稿第三十五页，共八十一页本讲稿第三十六页，共八十一页2.5 二维运动中的有序和混沌在物理学中使用计算机带来的一个基本好处是：能够处理不能

10、解析求解的系统。能够处理不能解析求解的系统。通常，数值结果同我们通过研究可解模型所发展出的直观图象在定性上一致。而定量的数值结果又实际兴趣。在少数情况下，计算机结果否定了我们的直观。这时数值结果对正确理解这种现象有根本的重要性。本讲稿第三十七页，共八十一页2.5 二维运动中的有序和混沌基本好处：能够处理不能解析求解不能解析求解的系统。通常在少数情况下数值结果直观图象定性实际兴趣定量结果数值结果直观图象本讲稿第三十八页，共八十一页否定经典图象令人惊奇的是，数值结果也否定了一些经典系统动力学的直观图象。牛顿以后，一切运转正常相对论指出高速运动下经典的直观图象不可行量子力学指出微观下经典的直观图象不

11、可行计算机的应用指出经典情形下，一些经典的直观图象不可行本讲稿第三十九页，共八十一页做二维运动的粒子本例题中，我们将通过用数值积分计算一个做二维运动的粒子的轨道以此为例研究这种令人惊奇的行为。考虑一个单位质量单位质量的粒子。它在一个位势V中作二维运动。假定如果粒子的能量足够低，它将永远被约束在V中。本讲稿第四十页，共八十一页若与两个坐标(x,y)共轭的动量是(px,py)，则Hamilton量的形式为粒子的轨道由坐标坐标和动量动量随时间的演化所确定。注意质量本讲稿第四十一页，共八十一页Hamilton方程坐标和动量的随时间的演化由四个耦合的一阶微分方程规定，即Hamilton方程本讲稿第四十二

12、页，共八十一页约束对于任何V，这些方程都能使能量E守恒，约束条件：因此，约束条件把轨道限制在嵌在四维相空间相空间中的一个三维流形上。除此之外，关于系统的演化很难作出别的什么除此之外，关于系统的演化很难作出别的什么一般性判断。一般性判断。本讲稿第四十三页，共八十一页可积系统能够对轨道作出进一步判断的一类重要的二维Hamilton量是可积的可积的Hamilton量。对于这种情况，除了能量以外，还有第二个坐标和动量的函数是运动常数运动常数。于是轨道被约束在相空间的一个二维二维流形上。两种熟悉的可积系统是可分离变量的系统和中心位势中心位势。本讲稿第四十四页，共八十一页可积系统1：可分离变量可分离变量的

13、情况下其中Vx,Vy是两个独立的函数。因此Hamilton量分离为两部分。两个方向的运动没有耦合。每一个方向上的Hamilton单独是一个运动常数。（等价于，Hx-Hx和E=Hx+Hy）H=Hx+HyHx=px2/2+Vx,Hy=py2/2+Vy本讲稿第四十五页，共八十一页可积系统2：中心位势在中心位势的场合向心力导致角动量角动量是第二个运动常数Hamilton量可以写为其中pr是r的共轭动量。本讲稿第四十六页，共八十一页可积性的困难可积系统的动力学是简单的。但是要找出这种简单性常常很不容易。没有一个普遍的解析方法可以判断，在任意一没有一个普遍的解析方法可以判断，在任意一个位势中是否存在第二个

14、运动常数，以及如果个位势中是否存在第二个运动常数，以及如果有的话如何求出它。有的话如何求出它。本讲稿第四十七页，共八十一页数值计算的固有局限数值计算看起来也不会改善情况。因为数值计算因为数值计算只给出只给出给定初始条件下的轨迹。给定初始条件下的轨迹。而这条轨迹即使在我们熟悉的情况下也可以是相当复杂的。比如Lissajous图形。（示例）本讲稿第四十八页，共八十一页本讲稿第四十九页，共八十一页本讲稿第五十页，共八十一页分析相空间通过对相空间的分析，可以得出一个根据轨迹来通过对相空间的分析，可以得出一个根据轨迹来检测可积性的办法。检测可积性的办法。考虑一个可分离变量的位势的情况。因为两个坐标上的运

15、动是独立的，在(x,px)平面和(y,py)平面上画出的轨迹看起来可能是本讲稿第五十一页，共八十一页可积性的标志一个粒子在一个可分离变量的二维位势中的轨迹在(x,px)和(y,py)平面上的样子。这些闭合回路的存在是系统可积性的标志。本讲稿第五十二页，共八十一页数值做法：截面根据轨迹可以得到一幅低维低维空间上的图。假设每当我们观察到一个坐标例如x通过零时，我们就在(y,py)平面上画出粒子的位置。如果x运动和y运动的周期是不可通约的（即它们的比值是一个无理数），那么随着轨迹前进，这些观测值将逐渐描绘出完整的(y,py)回路。如果它们的周期是可通约的（即比值是有理数）那么将得到沿回路的一系列分立

16、点。本讲稿第五十三页，共八十一页利用截面方法，可以只从轨迹出发，来研究利用截面方法，可以只从轨迹出发，来研究同任何给定的同任何给定的Hamilton量相联系的相空间的量相联系的相空间的拓扑性质。拓扑性质。本讲稿第五十四页，共八十一页在一个可积的Hamilton量的情形下，相空间的一般的拓扑性质，可以通过考虑中心位势中的运动来说明。对于固定的能量值和角动量值，径向运动被限制在两个转折点rin和rout之间。这两个r值是方程的两个解：pr=0 例如：卫星轨道本讲稿第五十五页，共八十一页这两个半径在(x,y)平面中定义了一个圆环面。对于给定的r值，能量守恒允许的径向运动的动量只能是这两个动量在(x,

17、y,pr)空间中确定了包含全部轨道的二维流形本讲稿第五十六页，共八十一页本讲稿第五十七页，共八十一页如果通过考虑x=0平面来画出一幅(y,py)截面图，我们将会得到两个闭合回路闭合回路闭合回路闭合回路。画截面比画流形简单的多。x=0本讲稿第五十八页，共八十一页可以证明，中心位势的相空间的环面拓扑结构，是一切可积系统所共有的。对于给定的运动常数，轨迹所在的流形叫做一个“不变环面”。本讲稿第五十九页，共八十一页环面的一般的截面图的样子本讲稿第六十页，共八十一页对可积性的扰动对可积性的扰动问题：在一个破坏掉系统的可积性的扰动下，一个可积系统的环面会发生什么变化？对于小扰动，椭圆型不动点周围的大部分环

18、面发生轻微的畸变，但仍保持它们的拓扑结构（KAM定理）。但是相空间中邻接区域则变成混沌的，给出的截面截面截面截面图看起来像是随机撒的一些点子图看起来像是随机撒的一些点子图看起来像是随机撒的一些点子图看起来像是随机撒的一些点子。在这些混沌区域中又嵌入别的椭圆型不动点和别的混沌区域，构成层次结构。本讲稿第六十一页，共八十一页Henon-Heiles位势对可积性的重大偏离必须用数值方法研究。进行这一研究的一个方便的例子是这个位势原来是Henon和Heiles在研究恒星穿过星系的轨道时引入的。恒星的运动被认为限制在二维平面上。本讲稿第六十二页，共八十一页本讲稿第六十三页，共八十一页Hamiltion运

19、动方程本讲稿第六十四页，共八十一页有了运动方程，给定一组初始条件（x,y,px,py）就可以计算出一条轨道。这里使用4阶Runge-Kutta方法。本讲稿第六十五页，共八十一页选择初条件1.考虑束缚轨道总能量 E1/6为了方便满足这个条件，我们直接选定E。本讲稿第六十六页，共八十一页选择初条件2.选择 x=03.选择 y1.确定y的范围：当能量全部转化为势能Epot时，y方向上距离原点最远。2.在这个范围内选择一个y值。4.选择py1.确定py的范围：动能全部在y方向时，py最大2.在这个范围内选择一个Py值。5.计算Px本讲稿第六十七页，共八十一页计算截面当轨道穿过x=0的面时，x变号。检测到x变号以后，计算轨道与x=0面的交点的精确位置。把运动方程转换为以x为自变量。把x的值作为步长，使用RK方法反向积分一步。本讲稿第六十八页，共八十一页本讲稿第六十九页，共八十一页本讲稿第七十页，共八十一页本讲稿第七十一页，共八十一页本讲稿第七十二页，共八十一页本讲稿第七十三页，共八十一页本讲稿第七十四页，共八十一页本讲稿第七十五页，共八十一页本讲稿第七十六页，共八十一页本讲稿第七十七页，共八十一页本讲稿第七十八页，共八十一页画图脚本1只输出一个截面图本讲稿第七十九页，共八十一页画图脚本2同时输出轨道和截面本讲稿第八十页，共八十一页本讲稿第八十一页，共八十一页