物理学中常微分方程初值问题的数值解法

1,第四章物理学中常微分方程初值问题的数值解法,2,4.1物理学中的常微分方程,一、力学中的例子,落体运动：作用力重力、阻力,牛顿方程：,阻尼振动：作用力弹性力、阻尼力,振动方程：,一阶常微分方程,二阶常微分方程,3,对比阻尼振动的标准形式：,受迫振动：,二阶常微分方程,二阶常微分方程,4,简谐振动实例：,弹簧振子：忽略各种阻力和弹簧质量的理想模型平衡位置：弹簧原长，选为原点；回复力：,单摆：忽略阻力和摆线质量，摆锤可视为质点，摆角小于5度,平衡位置：竖直位置；回复力矩：,由转动定理：,令,令,二阶齐次线性常微分方程,二阶齐次线性常微分方程,5,由转动定理：,不计阻力和弹簧质量，竖直弹簧振子的运动也是简谐振动,在平衡位置,取为原点o,回复力：,与水平弹簧振子一样也是简谐振动,令,扭摆：忽略各种阻力，忽略弹性杆的质量回复力矩：,动力学方程：,二阶齐次线性常微分方程,二阶齐次线性常微分方程,6,二、电学中的例子,放电电路,上电压：,上电压：,一阶常微分方程,1、放电电路,7,2、电磁振荡电路,与阻尼振动方程相当!问题：,二阶常微分方程,8,三、常微分方程数值解法的原理,1高阶化为一阶方程,只需研究一阶方程的解法。,9,2.泰勒级数,在级数中取若干项，得到近似方法：,二级欧拉法：取3项,龙格库塔法：取5项,10,4.2常微分方程初值问题的一级、二级欧拉近似法,一、一级欧拉近似法,1基本公式：,递推公式：,11,电路方程：,求,12,笔算此题：,13,求：,数值解：,14,其中,15,16,曲线,曲线,17,二、二级欧拉近似法,(1),将(2)代入(1):,18,(5),(4),19,(6),递推关系：,(6)可进一步写为：（由(6)式中的第一式解出,并代入(6)式中的第二式）,递推公式：,20,例3：RC电路：,一级欧拉近似结果,二级近似结果：,21,22,求：,数值解：,23,计算程序,real*8I(0:200),q(0:200),R,C,L,dt,I1(0:200)write(*,*)inputR=?,C=?,L=?read(*,*)R,C,Lopen(1,file=RCL2.dat)q(0)=1.0I(0)=0.0t0=0.0write(1,*)t0,I(0),q(0)dt=0.1,do10K=0,199f=-(R*I(K)/L+q(K)/(L*C)I1(K+1)=I(K)+f*dtf1=(1.)*(R*I1(K+1)/L+q(K)/(L*C)I(K+1)=0.5*(I(K)+I1(K+1)+f1*dt)q(K+1)=q(K)+I(K)*dtt=float(K+1)*dt10write(1,*)t,I(K+1),q(K+1)end,一级欧拉近似结果,二级欧拉近似结果,24,例5：RL暂态电路,一阶RL暂态电路如图所示，其中U0=311V，=314rad，R=10，L=500mH。求开关K合上后电流i(t)在区间0，0.01上的数值解。取t=0.001。,解：,根据电路理论可得初值问题为：,即：,25,用一级欧拉近似法公式如下：,用二级欧拉近似法公式如下：,26,此初值问题的精确解为：,计算结果,27,同例4类似，将二阶微分方程化成：,28,曲线,曲线,29,例7：导弹追踪问题,设位于坐标原点的甲舰向位于x轴上点A(1,0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程，又乙舰行驶多远时，导弹将它击中？,解法一（解析法）,假设导弹在t时刻的位置为P(x(t),y(t)，乙舰位于Q(1,v0t),由于导弹头始终对准乙舰，故此时直线PQ,即有：,就是导弹的轨迹曲线弧OP在点P处的切线,30,又根据题意，弧OP的长度为AQ的5倍,联系(1)、(2)式消去t可得：,初值条件为：,解即为导弹的运行轨迹,当时，即当乙舰航行到点处时被导弹击中。,被击中时间为:。,若v0=1，则在t=0.21时被击中。,31,解法二（数值方法）,令，将方程(3)化为一阶微分方程。,初值条件为：,由一阶欧拉近似公式可得：,32,programmainreal*8y(0:200),z(0:200),x0,dx,xopen(1,file=DDoula.dat)y(0)=0.0z(0)=0.0 x0=0.0write(1,*)x0,y(0),z(0)dx=0.01do10k=0,100y(k+1)=y(k)+z(k)*dxz(k+1)=z(k)+sqrt(1+z(k)*2)*dx/(5-5*x)x=float(k+1)*dx10write(1,*)x,y(k+1)end,一级欧拉近似计算程序如下：,33,模拟导弹追踪轨迹：,试用二级欧拉近似求解导弹追踪问题？,34,作业,EX4-1：假设高楼顶上有一物体掉下来，在下落过程中，受到空气的阻力与下落速度的平方成正比，正比系数为，物体质量，重力加速度为9.81m/s2，设楼顶处为原点，初速度为0，试用一级欧拉近似法、二级欧拉近似法编程求解下落位移、速度随时间的变化关系（）。,35,若是受迫振动，振动方程为试取，求解微分方程并给出曲线。,36,受迫振动示意图,37,EX4-3：如右图所示的一单摆，摆长1.016m，重力加速度9.81m/s2，摆角所满足的微分方程可表示为求用二级欧拉法近似求摆角与时间的曲线关系。,38,4.3龙格-库塔法,一、龙格库塔公式,类似于欧拉法公式的推导，利用二元函数的Taylor,可以导出一类具有四阶精度的龙格-库塔公式。,不进行推导,展开式，使计算公式的局部截断误差取为,39,常用格式：,40,解：精确解：,方程初值问题,令,计算得：,而精确解为：,41,计算程序,精确解：,42,subroutinerk4(y,dydx,x,h)real*8y,dydx,yt,dyt,dym,dyk,x,xh,h6,hhhh=h*0.5!h6=h/6!xh=x+hh!,subroutinederivs1(x,y,dydx)real*8x,y(1),dydx(1)dydx(1)=-y(1)+x*2.+1returnend,求导数,43,y=y+h6*(dydx+2*dyt+2*dym+dyk)returnend,计算y1,计算K4存在dyk,计算K3存在dym,计算K2存在dyt,yt=y+hh*dydxcallderivs1(xh,yt,dyt)yt=y+hh*dytcallderivs1(xh,yt,dym)yt=y+h*dymcallderivs1(x+h,yt,dyk),44,二、常微分方程组的求解,含有两个未知函数的一阶常微分方程组初值问题,可以写作如下形式：,为了用龙格库塔公式求解以上初值问题，,方程组写成向量形式，记,通常将,45,则一阶常微分方程组初值问题式可以写为,46,解：此初值问题的精确解为：,常微分方程组初值问题,47,48,因此有：,以上数值解与精确解的误差为,同理可计算出：,49,例7:解法三（参数方程）,设时刻t乙舰的坐标为(X(t)，Y(t)，导弹的坐标为(x(t)，y(t)。,3、因乙舰以速度v0沿直线x=1运动，设v0=1，则w=5，X=1，Y=t,1、设导弹速度恒为w，则,2、由于弹头始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量，,即,由(1),(2)消去可得：,50,因此导弹运动轨迹的参数方程为：,解此参数方程即可得导弹追踪轨迹。,51,三、高阶常微分方程的求解,上式可以转化为一阶微分方程组的初值问题。,m阶微分方程初值问题,52,例如，考虑二阶微分方程的初值问题：,也可采用降阶方法求解！,53,例10：对二阶微分方程初值问题,应用龙格库塔公式，取,计算可得：,54,该初值问题的精确解为,55,作业,56,57,慢跑者与狗,一个慢跑者在平面上沿椭圆以恒定的速率v=1跑步,设椭圆方程为:x=10+20cost,y=20+5sint.突然有一只狗攻击他.这只狗从原点出发