Subroutine rk4_d22(d1x,d2x,d1y,d2y,ti,tf,xi,xf,yi,yf) c========================================================== c rk4_d22.f: Solution of the second-order 2D ODE c method: Runge-Kutta 4th-order c written by: Alex Godunov c last revision: 21 October 2006 c---------------------------------------------------------- c input ... c d1x(t,x,y)- function dx/dt (supplied by a user) c d2x(t,x,y)- function d2x/dt2 (supplied by a user) c d1y(t,x,y)- function dy/dt (supplied by a user) c d2y(t,x,y)- function d2y/dt2 (supplied by a user) c where x(2) and y(2) (x(1)-position, x(2)-speed, etc.) c ti - initial time c tf - time for a solution c xi(2) - initial position and speed for x component c yi(2) - initial position and speed for y component c c output ... c xf(2) - solutions (x position and speed) at point tf c yf(2) - solutions (y position and speed) at point tf c c the program written as an educational program c therefore, the text is clear but far from concise c the program can be easily rewritten in a compact form c========================================================== implicit none Real*8 d1x, d2x, d1y, d2y, ti, tf Real*8 xi(2), xf(2), yi(2), yf(2) Real*8 h,t, x(2), y(2) Real*8 k1x(2),k2x(2),k3x(2),k4x(2),k1y(2),k2y(2),k3y(2),k4y(2) h = tf-ti t = ti c k1 k1x(1) = h*d1x(t,xi,yi) k1x(2) = h*d2x(t,xi,yi) k1y(1) = h*d1y(t,xi,yi) k1y(2) = h*d2y(t,xi,yi) c k2 x(1) = xi(1)+k1x(1)/2.0 x(2) = xi(2)+k1x(2)/2.0 y(1) = yi(1)+k1y(1)/2.0 y(2) = yi(2)+k1y(2)/2.0 k2x(1) = h*d1x(t+h/2.0,x,y) k2x(2) = h*d2x(t+h/2.0,x,y) k2y(1) = h*d1y(t+h/2.0,x,y) k2y(2) = h*d2y(t+h/2.0,x,y) c k3 x(1) = xi(1)+k2x(1)/2.0 x(2) = xi(2)+k2x(2)/2.0 y(1) = yi(1)+k2y(1)/2.0 y(2) = yi(2)+k2y(2)/2.0 k3x(1) = h*d1x(t+h/2.0,x,y) k3x(2) = h*d2x(t+h/2.0,x,y) k3y(1) = h*d1y(t+h/2.0,x,y) k3y(2) = h*d2y(t+h/2.0,x,y) c k4 x(1) = xi(1)+k3x(1) x(2) = xi(2)+k3x(2) y(1) = yi(1)+k3y(1) y(2) = yi(2)+k3y(2) k4x(1) = h*d1x(t+h,x,y) k4x(2) = h*d2x(t+h,x,y) k4y(1) = h*d1y(t+h,x,y) k4y(2) = h*d2y(t+h,x,y) c result xf(1) = xi(1)+k1x(1)/6.0+k2x(1)/3.0+k3x(1)/3.0+k4x(1)/6.0 xf(2) = xi(2)+k1x(2)/6.0+k2x(2)/3.0+k3x(2)/3.0+k4x(2)/6.0 yf(1) = yi(1)+k1y(1)/6.0+k2y(1)/3.0+k3y(1)/3.0+k4y(1)/6.0 yf(2) = yi(2)+k1y(2)/6.0+k2y(2)/3.0+k3y(2)/3.0+k4y(2)/6.0 Return End