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

第二章常微分方程的初值问题本章要研究的物理问题：经典动力学中的有序和混沌本章内容4简单方法123多步法和隐式法龙格库塔法稳定性问题4动力学中的有序和混沌52.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))endf=@(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(m−1)/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-Heiles系统这个位势原来是由Henon和Heiles在研究恒星穿过星系的轨道时引入的。它具有三重对称性，并且在原点处位势为零,当坐标值大时位势无限增大。对于小于1/6的能量，轨道被限制在一个等边三角形之内Henon-Heiles位势所确定的Hamilton方程为该问题不能精确求解，所以必须用数值方法研究.二维粒子的运动轨迹(x,px)平面和(y,py)平面的轨迹图