Fortran 2003 was a significant update to Fortran 95 (probably something on the order of the update from C to C++). It brought Fortran into the modern computing world with the addition of object oriented programming features such as type extension and inheritance, polymorphism, dynamic type allocation, and type-bound procedures.
The following example shows an implementation of an object oriented class to perform numerical integration of a user-defined function using a Runge-Kutta method (in this case RK4). First, an abstract rk_class is defined, which is then extended to create the RK4 class by associating the deferred step method to rk4. Other RK classes could also be created in a similar manner. In order to use the class, it is only necessary to set n (the number of variables in the state vector) and f (the user-supplied derivative function), and then call the integrate method.
module rk_module
use, intrinsic :: iso_fortran_env, wp => real64 !double precision reals
implicit none
real(wp),parameter :: zero = 0.0_wp
!main integration class:
type,abstract,public :: rk_class
private
!user specified number of variables:
integer :: n = 0
!user-specified derivative function:
procedure(deriv_func),pointer :: f => null()
contains
procedure,non_overridable,public :: integrate !main integration routine
procedure(step_func),deferred :: step !the step routine for the selected method
end type rk_class
!extend the abstract class to create an RK4 method:
! [all we need to do is set the step function]
type,extends(rk_class),public :: rk4_class
contains
procedure :: step => rk4
end type rk4_class
interface
subroutine deriv_func(me,t,x,xdot) !derivative function
import :: rk_class,wp
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),dimension(me%n),intent(out) :: xdot
end subroutine deriv_func
subroutine step_func(me,t,x,h,xf) !rk step function
import :: rk_class,wp
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),intent(in) :: h
real(wp),dimension(me%n),intent(out) :: xf
end subroutine step_func
end interface
contains
subroutine integrate(me,t0,x0,h,tf,xf)
! main integration routine
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t0
real(wp),dimension(me%n),intent(in) :: x0
real(wp),intent(in) :: h
real(wp),intent(in) :: tf
real(wp),dimension(me%n),intent(out) :: xf
real(wp) :: t,dt,t2
real(wp),dimension(me%n) :: x
logical :: last
if (h==zero) then
xf = x0
else
t = t0
x = x0
dt = h
do
t2 = t + dt
last = ((dt>=zero .and. t2>=tf) .or. & !adjust last time step
(dt<zero .and. t2<=tf)) !
if (last) dt = tf-t !
call me%step(t,x,dt,xf)
if (last) exit
x = xf
t = t2
end do
end if
end subroutine integrate
subroutine rk4(me,t,x,h,xf)
! Take one Runge Kutta 4 integration step: t -> t+h (x -> xf)
implicit none
class(rk4_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),intent(in) :: h
real(wp),dimension(me%n),intent(out) :: xf
!local variables:
real(wp),dimension(me%n) :: f1,f2,f3,f4
real(wp) :: h2
!parameters:
real(wp),parameter :: half = 0.5_wp
real(wp),parameter :: six = 6.0_wp
h2 = half*h
call me%f(t,x,f1)
call me%f(t+h2,x+h2*f1,f2)
call me%f(t+h2,x+h2*f2,f3)
call me%f(t+h,x+h*f3,f4)
xf = x + h*(f1+f2+f2+f3+f3+f4)/six
end subroutine rk4
end module rk_module
An example use of this module follows. Here, we are integrating the equations of motion of a spacecraft in Earth orbit, using two-body assumptions: \(\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}\)
The rk4_class is extended to create a new spacecraft class, which contains the data required by the derivative routine (in this case, simply the central body gravitational parameter mu). A more complicated version could also include spherical harmonic gravity, solar radiation pressure, atmospheric drag, engine thrust, etc. Note that the select type statement is necessary in the derivative routine in order to access the data from the extended type.
program main
use rk_module
!spacecraft propagation type:
! extend the rk class to include data used in the deriv routine
type,extends(rk4_class) :: spacecraft
real(wp) :: mu = zero ! central body gravitational
! parameter (km^3/s^2)
end type spacecraft
integer,parameter :: n=6 !number of state variables
type(spacecraft) :: s
real(wp) :: t0,tf,x0(n),dt,xf(n)
!constructor:
s = spacecraft(n=n,f=twobody,mu=398600.436233_wp) !main body is Earth
!initial conditions:
x0 = [10000.0_wp,10000.0_wp,10000.0_wp,& !initial state [x,y,z] (km)
1.0_wp,2.0_wp,3.0_wp] ! [vx,vy,yz] (km/s)
t0 = zero !initial time (sec)
dt = 10.0_wp !time step (sec)
tf = 1000.0_wp !final time (sec)
call s%integrate(t0,x0,dt,tf,xf)
write(*,'(A/,*(F15.6/))') 'Final state:',xf
contains
subroutine twobody(me,t,x,xdot)
! derivative routine for two-body orbit propagation
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),dimension(me%n),intent(out) :: xdot
real(wp),dimension(3) :: r,v,a_grav
real(wp) :: rmag
select type (me)
class is (spacecraft)
r = x(1:3)
v = x(4:6)
rmag = norm2(r)
a_grav = -me%mu/rmag**3 * r !accel. due to gravity
xdot(1:3) = v
xdot(4:6) = a_grav
end select
end subroutine twobody
end program main
The result is:
Final state: 10667.963305 11658.055962 12648.148619 0.377639 1.350074 2.322509
See also
- Object-Oriented Programming in Fortran 2003 Part 1: Code Reusability, PGI Insider, Feb. 2011.
- Object-Oriented Programming in Fortran 2003 Part 2: Data Polymorphism, PGI Insider, June 2011.
- Fortran Astrodynamics Toolkit (Github)
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]
