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

⒈ 弹簧振子：忽略各种阻力和弹簧质量的理想模型 平衡位置：弹簧原长，选为原点 ；回复力： 模拟导弹追踪轨迹： 解法二（数值方法） 令 ，将方程(3)化为一阶微分方程。 初值条件为： 由一阶欧拉近似公式可得： program main real*8 y(0:200), z(0:200), x0, dx, x open(1,file=DDoula.dat) y(0)=0.0 z(0)=0.0 x0=0.0 write(1,*) x0,y(0),z(0) dx=0.01 do 10 k=0,100 y(k+1)=y(k)+z(k)*dx z(k+1)=z(k)+sqrt(1+z(k)**2)*dx/(5-5*x) x=float(k+1)*dx 10 write(1,*) x, y(k+1) end 一级欧拉近似计算程序如下： 试用二级欧拉近似求解导弹追踪问题？ 作业 EX4-1：假设高楼顶上有一物体掉下来，在下落过程中，受到空气的阻力与下落速度的平方成正比，正比系数为 ，物体质量 ，重力加速度为9.81m/s2，设楼顶处为原点，初速度为0，试用一级欧拉近似法、二级欧拉近似法编程求解下落位移、速度随时间的变化关系（ ）。 EX4-2：已知弹簧振子阻尼振动的方程为： 时 取0.5秒，试分别取 用二级欧拉法近似画出 曲线 若是受迫振动，振动方程为 试取 ，求解微分方程并给出 曲线。 受迫振动示意图 EX4-3：如右图所示的一单摆，摆长 1.016m，重力加速度9.81m/s2，摆角 所满足的微分方程可表示为 求用二级欧拉法近似求摆角 与时间 的曲线关系。 §4.3 龙格-库塔法 求： 已知: 一、龙格－库塔公式 类似于欧拉法公式的推导，利用二元函数的Taylor 可以导出一类具有四阶精度的龙格-库塔公式。 不进行推导 展开式，使计算公式的局部截断误差取为 常用格式： 解：精确解： 例8：取步长 ，用龙格-库塔公式求解常微分 方程初值问题 输入初值 ， , 令 计算得： 而精确解为： 计算结果比欧拉方法 精度要高！(见程序） program main real*8 y,dydx,x x=0.0d0 y=1.0d0 dydx=0.0d0 h=0.1d0 do j=1,10 call rk4(y,dydx,x,h) x=x+h yy=-2*exp(-x)+x**2.-2*x+3 write(*,*) ‘real=’,x , y,yy call derivs1(x,y,dydx) end do end 求下一个点导数 精确解 龙格库塔法求解存在y中 计算程序 精确解： subroutine rk4(y,dydx,x,h) real*8 y,dydx,yt,dyt,dym,dyk,x,xh,h6,hh hh=h*0.5 ! h6=h/6 ! xh=x+hh ! subroutine derivs1(x,y,dydx) real*8 x,y(1),dydx(1) dydx(1)=-y(1)+x**2.+1 return end 求导数 y=y+h6*(dydx+2*dyt+2*dym+dyk) return end 计算y1 计算K4存在dyk 计算K3存在dym 计算K2存在dyt yt=y+hh*dydx call derivs1(xh,yt,dyt) yt=y+hh*dyt call derivs1(xh,yt,dym) yt=y+h*dym call derivs1(x+h,yt,dyk) * * 第四章 物理学中常微分方程初值问题的数值解法 §4.1 物理学中的常微分方程 §4.2 常微分方程初值问题的一级、二级欧拉近似法 §4.3 龙格-库塔法 §4.1 物理学中的常微分方程 一、力学中的例子 落体运动：作用力—重力、阻