program main !==================================================================== ! Initial value problem for a single 2nd-order ODE ! Call Methods: Runge-Kutta 4th order ! comment: can be used for a system of two first order ODEs !==================================================================== implicit none double precision ti,xi,vi,tf,xf,vf double precision d1x,d2x double precision ei, dt double precision tmax, k, m integer i external d1x, d2x ! open (unit=6,file="result12.dat") !* initial conditions ti = 0.0 xi = 0.0 vi = 1.0 !* time step and max time dt = 0.1 tmax = 2.00 ei = vi**2/2 + xi**2/2 ! energy for a harmonic oscillator !* print the header and initial conditions write (6,*) ' 2-nd order single ODE ' write (6,*) ' Method: Runge-Kutta 4th order' write (6,*) ' t x(t) v(t) energy' write (6,100) ti, xi, vi, ei do while(ti.le.tmax) tf = ti + dt call RK4D12(d1x,d2x,ti,xi,vi,tf,xf,vf) ei = vf**2/2 + xf**2/2 write (6,100) tf, xf, vf, ei ti = tf xi = xf vi = vf end do 100 format(4f10.5) end program main Function d1x(t,x,v) !-------------------------------------- ! function dx/dt !-------------------------------------- implicit none double precision d1x, t, x, v d1x = v end Function d1x Function d2x(t,x,v) !-------------------------------------- ! function d2x/dt2 !-------------------------------------- implicit none double precision d2x, t, x, v ! simple harmonic oscillator double precision, parameter:: k=1.0, m=1.0 d2x = (-1.0*k/m)*x end Function d2x Subroutine RK4D12(d1x,d2x,ti,xi,vi,tf,xf,vf) !=========================================================== ! Solution of a single second-order ODE d2x/dt2=f(t,x,dx/dt) ! Method: Runge-Kutta 4th-order ! Alex G. February 2010 !----------------------------------------------------------- ! input ... ! d1x(t,x,v)- function dx/dt (supplied by a user) ! d2x(t,x,v)- function d2x/dt2 (supplied by a user) ! ti - initial time ! xi - initial position ! vi - initial velocity (first order derivative) ! tf - final time (find solution at this time) ! output ... ! xf - position at point tf ! vf - velocity at point tf !=========================================================== implicit none double precision d1x,d2x,ti,xi,vi,tf,xf,vf double precision h,t,k1x,k2x,k3x,k4x,k1v,k2v,k3v,k4v h = tf-ti t = ti k1x = h*d1x(t,xi,vi) k1v = h*d2x(t,xi,vi) k2x = h*d1x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0) k2v = h*d2x(t+h/2.0,xi+k1x/2.0,vi+k1v/2.0) k3x = h*d1x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0) k3v = h*d2x(t+h/2.0,xi+k2x/2.0,vi+k2v/2.0) k4x = h*d1x(t+h,xi+k3x,vi+k3v) k4v = h*d2x(t+h,xi+k3x,vi+k3v) xf = xi + (k1x + 2.0*(k2x+k3x) + k4x)/6.0 vf = vi + (k1v + 2.0*(k2v+k3v) + k4v)/6.0 end subroutine RK4D12