第二章-常微分方程的初值问题

第二章常微分方程的初值问题,本章要研究的物理问题：经典动力学中的有序和混沌,本章内容,4,简单方法,1,2,3,多步法和隐式法,龙格库塔法,稳定性问题,4,动力学中的有序和混沌,5,2.0引子,常微分方程是物理学中经常碰到，以一维运动粒子为例,一般形式为M个藕合的一阶方程,实际问题往往涉及不止一类问题，例如对偏微分方程分离变量时，同时得到分离开来的初值问题与本征问题。,本章只讨论初值问题,常微分方程定解问题的分类,给定待求函数在某个初始点上的值,在自变量的两个端点上对待求函数施加约束,含有待定参数，只有在参数取特定值时，方程才有非零解,2.1简单方法,1).离散化：将区间0,1分为N个等间隔的子区间，每个区间宽度为h=1/N,待求问题,求x=1处y的值,策略,2).寻找一个递推关系，把yn同yn-1,yn-2,联系起来。,局部误差为O(h2),在xn点，将微分方程左边的微分用向前差分公式替代,欧拉法,得到递推关系,全局误差为NO(h2)O(h),精度太低！,例子：牛顿动力学方程,一维运动的粒子方程为,用欧拉法写出具体的递推关系,欧拉法的Matlab实现,标量形式,矢量形式,forn=1:Nx(n+1)=x(n)+h*p(n)/mp(n+1)=p(n)+h*g(x(n),t(n)end,f=(x,t)x(2)/m,g(x(1),t)forn=1:Nz(n+1,:)=z(n,:)+h*f(z(n,:),t(n)end,这个公式在已知f的解析形式时很好用。阶数增大时，变得复杂。,泰勒级数法,将yn+1在yn附近做泰勒展开,而又已知,从而有,练习：推导三阶泰勒级数法,2.2多步法和隐式法,在子区间xn,xn+1将微分方程写成积分形式,关键在于对积分号下的f(x,y)取合理的近似,欧拉法实际上将f(x,y)近似为f(xn,yn),利用f在xn-1和xn处的值通过线型插值得到f在积分区域的值,多步法,为了获得更高的精度，可以将yn+1不仅仅是同yn，而且还同更早的点，如yn-1,yn-2等点相联系,例如:,代入积分公式得Adams-Bashforth二步积分法,类似的，Adams-Bashforth四步积分法为,多步法有个小麻烦：最初几个格点的启动值需要通过其它的积分法如欧拉法、泰勒级数法来获得,隐式法,在区间xn,xn+1上过fn，fn+1两点对f插值，得,带入积分公式,得,隐式法意味着在每一个积分步都必须解一个方程，会非常耗时。,一个特别简单的情况是，如果f对于y是线性的，即f(x,y)=g(x)y，则上面的方程可以为,得到其显式解,利用fn-1,fn,fn+1三个值在区间xn,xn+1对f进行二次多项式插值，得到隐式递推关系,Adams-Moulton方法多步隐式法,带入积分公式,得,先通过显式法“预报”一个yn+1的值，然后利用隐式法来校正得到一个更精确的值。这样的算法有个优点，它可以持续监控积分的精度。,用三次多项式插值，得到相应的三步法为,隐式法很少直接使用，它通常用在预报校正算法中,一个常常使用的具有局域误差O(h5)的预报校正算法：,显式的Adams-Bashforth四步法,Adams-Moulton三步法,练习：简谐振子方程,2.3Runge-Kutta法,思想：从欧拉法说起,在子区间xn,xn+1将微分方程写成积分形式,中值定理,平均斜率,欧拉法就是用xn点的斜率近似xn,xn+1区间的平均斜率Kave,改进的欧拉法,K1、K2分别是xnxn+1点处的斜率,但是由于yn+1待定，因此需要做“预报”,用xnxn+1两点斜率的平均值来近似xn,xn+1平均斜率,思想：为了提高精度，多取几点的斜率值作为加权平均当作平均斜率,其中i(i=1,2,.,m)和ij(i=2,3,.,m且ji)是待定参数,Runge-Kutta法,将上式展开到O(hm)，得到m个方程，而有m+m(m1)/2个待定参数i,ij.所以还有灵活选择的空间,以m=2为例,Runge-Kutta法公式为,其中,将K2做泰勒展开到O(h2)项，得,利用,将yn+1在xn点附近做泰勒展开至O(h2)项，得,两式相比，有,可选解为,或,若取,得二阶Runge-Kutta法,二阶Runge-Kutta法与二阶泰勒级数法比较,可以看出，相比二阶泰勒级数法而言，二阶Runge-Kutta法适用性更广，使用也更为方便,三阶Runge-Kutta法,最为常用的是四阶Runge-Kutta法,通常认为四阶Runge-Kutta法在效率和精度间达到了最好的平衡,练习：四阶Runge-Kutta法处理一阶常微分方程,f=(x,y)(-x*y);y(1)=1;fori=1:nK1=f(x(i),y(i);K2=f(x(i)+h/2,y(i)+h*K1/2);K3=f(x(i)+h/2,y(i)+h*K2/2);K4(i)=f(x(i)+h,y(i)+h*K3);y(i+1)=y(i)+1/6*h*(K1+2*K2+2*K3+K4);end,练习：二阶Runge-Kutta法处理二阶常微分方程,2.4算法的稳定性,以欧拉法的一个简单扩展为例,对微分方程进行积分时，一个首要的考虑是所用算法的数值稳定性，也就是说，舍入误差或数值计算中的其它误差能被放大的程度。,为了启动上面的递推关系，还需要y1的值，这可以由泰勒展开获得,当把这一方法应用到下述问题,这个问题的解析解是y=e-x,设解为指数形式y=Arn，代入递推关系得,正根略小于1，它对应于我们要找的按指数减小的解，但是负根略小于-1，因而它对应于一个虚假的解,上述方程的解为,其大小随n增大并且在格点上逐点发生振荡。,线性差分方程的通解正是这两个指数解的一个线性组合。虽然可以精心安排初值y0和y1，使得当x值小时只呈现指数衰减的解，但是在递推过程中的舍入误差将会掺入一个小的“坏”解，它最后将增长到在解中占压倒地位。,一个好的经验规则是，每当积分一个随着迭代过程急剧衰减的解时，就应当小心不稳定性和舍入误差。,例1.强迫钟摆,一根长度为l钟摆被限制在一个垂直的平面内，在强迫外力fd和阻力fr的作用下振荡运动。钟摆的运动可以通过Newton方程来描述,其中fg=-mgsin是重力在运动方向的分力，a=ld2/dt2是沿切线方向的加速度，是杆和垂线的夹角。,2.5动力学中的有序和混沌,那么Newton运动方程就可以写成如下形式,假设强迫外力为,阻力为,并且取(l/g)1/2为时间单位。,其中,该运动方程就可以化为一阶方程组,设,q=0.5,b=0.9,0=2/3,t=100。在这种情况下,钟摆运动是一个有序的周期运动,杆和垂线的夹角、角速度随时间演化过程,角速度与夹角的轨迹,q=0.5,b=1.15,0=2/3,t=1000，在这种情况下钟摆运动是一个无序运动,呈现出分形的特征,杆和垂线的夹角、角速度随时间演化过程,角速度与夹角的轨迹,例2.自激振动范德波尔方程,范德波尔描述非线性有阻尼的自激动振动系统,其中是一个正的小量,自激系统能将非振动的能源通过系统本身的反馈调节吸收进来，以补充被损耗的能量。自激系统的例子包括心脏等。,VDP方程不能准确的表示心脏振荡图形的细节，为了更准确的表示心脏的波动，FitzHugh对VDP方程进行了修改，提出了所谓的BVP方程，形式为,例3.化学振荡BZ反应,CA,CB，CC，CD分别为四种化合物浓度。,例4.二维粒子的运动,考虑一个单位质量的粒子,它在一个位势V中作二维运动。,设在t时刻,粒子的平面坐标为(x,y),它的共轭动量为(px,py),则Hamilton量的形式为,粒子的轨迹就由坐标和动量随时间的演化由如下Hamilton方程规定,该约束条件把粒子运动的轨道限制在四维相空间的三维流形上.,它是四个耦合的一阶微分方程组.这些方程使能量E守恒,即满足约束条件,在可分离变量的系统中,位势V具有(x,y)的可分离形式,即V(x,y)=Vx(x)+Vy(y)，此时Hamilton方程可写为,可分离变量的系统,和,由此可见,x方向和y方向的运动是相互不耦合的，每一个Hamilton量单独都是一个运动常数，其中,作为一个例子,考虑两个坐标方向的运动都是简谐运动的情况,即位势,对应的Hamilton方程为,和,该问题是精确可解的。我们可以通过直接积分,求得满足初始条件(x0,y0,px0,py0)的解为,(x,px)平面和(y,py)平面的轨迹图,(x,y)平面的运动轨迹,在Henon-Heiles位势系统中，位势的形式为,Henon-Heile