Nyström Integration Methods
Runge-Kutta Nyström methods can be used to integrate second-order ordinary differential equations of the form:
- (a) \(\ddot{\mathbf{x}}=f(t,\mathbf{x})\)
There are also versions for the more general form:
- (b) \(\ddot{\mathbf{x}}=f(t,\mathbf{x},\dot{\mathbf{x}})\)
For second-order systems, Nyström methods are generally more efficient than the more familiar Runge-Kutta methods. The following algorithm is a 4th order Nyström method from Reference [1], encapsulated into a simple Fortran module. It integrates equations of type (b). Several others are also given in the same paper. A high-accuracy variable step size Nyström method for systems of type (a) can also be found in Reference [2].
module nystrom
use, intrinsic :: iso_fortran_env, wp=>real64
implicit none
private
! derivative function prototype
abstract interface
function deriv(n,t,x,xd) result(xdd)
import :: wp
integer,intent(in) :: n
real(wp),intent(in) :: t
real(wp),dimension(n),intent(in) :: x
real(wp),dimension(n),intent(in) :: xd
real(wp),dimension(n) :: xdd
end function deriv
end interface
public :: nystrom_lear_44
contains
subroutine nystrom_lear_44(n,df,t0,x0,xd0,dt,xf,xdf)
! 4th order Nystrom integration routine
! Take a step from t0 to t0+dt (x0,xd0 -> xf,xdf)
implicit none
integer,intent(in) :: n !number of equations
procedure(deriv) :: df !derivative routine
real(wp),intent(in) :: t0 !initial time
real(wp),dimension(n),intent(in) :: x0,xd0 !initial state
real(wp),intent(in) :: dt !time step
real(wp),dimension(n),intent(out) :: xf,xdf !final state
real(wp) :: t
real(wp),dimension(n) :: x,xd,k1,k2,k3,k4
real(wp),parameter :: s5 = sqrt(5.0_wp)
real(wp),parameter :: delta2 = (5.0_wp-s5)/10.0_wp
real(wp),parameter :: delta3 = (5.0_wp+s5)/10.0_wp
real(wp),parameter :: delta4 = 1.0_wp
real(wp),parameter :: a1 = (3.0_wp-s5)/20.0_wp
real(wp),parameter :: b1 = 0.0_wp
real(wp),parameter :: b2 = (3.0_wp+s5)/20.0_wp
real(wp),parameter :: c1 = (-1.0_wp+s5)/4.0_wp
real(wp),parameter :: c2 = 0.0_wp
real(wp),parameter :: c3 = (3.0_wp-s5)/4.0_wp
real(wp),parameter :: alpha1 = 1.0_wp/12.0_wp
real(wp),parameter :: alpha2 = (5.0_wp+s5)/24.0_wp
real(wp),parameter :: alpha3 = (5.0_wp-s5)/24.0_wp
real(wp),parameter :: alpha4 = 0.0_wp
real(wp),parameter :: ad1 = (5.0_wp-s5)/10.0_wp
real(wp),parameter :: bd1 = -(5.0_wp+3.0_wp*s5)/20.0_wp
real(wp),parameter :: bd2 = (3.0_wp+s5)/4.0_wp
real(wp),parameter :: cd1 = (-1.0_wp+5.0_wp*s5)/4.0_wp
real(wp),parameter :: cd2 = -(5.0_wp+3.0_wp*s5)/4.0_wp
real(wp),parameter :: cd3 = (5.0_wp-s5)/2.0_wp
real(wp),parameter :: beta1 = 1.0_wp/12.0_wp
real(wp),parameter :: beta2 = 5.0_wp/12.0_wp
real(wp),parameter :: beta3 = 5.0_wp/12.0_wp
real(wp),parameter :: beta4 = 1.0_wp/12.0_wp
if (dt==0.0_wp) then
xf = x0
xdf = xd0
else
t = t0
x = x0
xd = xd0
k1 = dt*df(n,t,x,xd)
t = t0+delta2*dt
x = x0+delta2*dt*xd0+a1*dt*k1
xd = xd0+ad1*k1
k2 = dt*df(n,t,x,xd)
t = t0+delta3*dt
x = x0+delta3*dt*xd0+dt*(b1*k1+b2*k2)
xd = xd0+bd1*k1+bd2*k2
k3 = dt*df(n,t,x,xd)
t = t0+dt
x = x0+dt*xd0+dt*(c1*k1+c3*k3)
xd = xd0+cd1*k1+cd2*k2+cd3*k3
k4 = dt*df(n,t,x,xd)
xf = x0+dt*xd0+dt*(alpha1*k1+alpha2*k2+alpha3*k3)
xdf = xd0+beta1*k1+beta2*k2+beta3*k3+beta4*k4
end if
end subroutine nystrom_lear_44
end module nystrom
References
- W. M. Lear, "Some New Nyström Integrators", 78-FM-11, JSC-13863, Mission Planning and Analysis Division, NASA Johnson Space Center, Feb. 1978.
- R. W. Brankin, I. Gladwell, J. R. Dormand, P. J. Prince, and W. L. Seward, "Algorithm 670: a Runge-Kutta-Nyström code", ACM Transactions on Mathematical Software, Volume 15 Issue 1, March 1989. [code]