I want to solve ODEs of the type
in modern fortran. I want to write the integrator (for example 4th-order Runge-Kutta) to be quite general, such that the same integrator can be used for different right-hand sides of the equation, including at least these different cases:
- f(x, t) is a function, where both x and the return value are scalars
- f(x, t) is a function, where both x and the return value are arrays (of the same shape)
- f(x, t) is a type-bound procedure (or similar), where the purpose of the derived type is to interpolate some underlying data to the position and time given by (x, t)
Based on an answer to a similar question I have implemented the code included below, which seems to work as expected.
My problem is that I would like to keep rk4(x, t, h, f) general, such that x in f(x,t) can be an assumed-shape array (or even a scalar, ideally), but at the same time I would like to be able to specify that x has for example dimension(2) when I actually implement a function to interpolate some 2D data. However, if I try to modify the function evaluate in the example below such that x has dimension(2), then I get an error complaining about interface mismatch when I try to compile. Is there a way around this problem?
module interpolator_module
implicit none
integer, parameter :: WP = kind(1.0D0)
interface
! This is the general form of the right-hand side of an ODE
function rhs(x, t) result( val )
import :: WP
real(WP), dimension(:), intent(in) :: x
real(WP), intent(in) :: t
real(WP), dimension(size(x)) :: val
end function
end interface
type interpolator_type
! This type would in practice store arrays,
! of discrete data to be interpolated.
real(WP) :: stored_data
procedure(rhs), nopass, pointer :: eval
contains
procedure :: init
endtype
class(interpolator_type), pointer :: interpolator
contains
subroutine init( this, stored_data )
implicit none
class(interpolator_type), target :: this
real(WP) :: stored_data
this % stored_data = stored_data
this % eval => evaluate
interpolator => this
end subroutine
function evaluate(x, t) result( val )
implicit none
real(WP), dimension(:), intent(in) :: x
real(WP), intent(in) :: t
real(WP), dimension(size(x)) :: val
! This is where interpolation would happen
val = interpolator % stored_data * x
end function
end module
program main
use interpolator_module, only : interpolator_type
implicit none
integer, parameter :: WP = kind(1.0D0)
type(interpolator_type) :: interp
real(WP), dimension(2) :: x
real(WP) :: t, h
! initialise interpolator with some data
call interp % init(-0.1_WP)
x = (/ 2.0_WP, 1.0_WP /)
t = 0.0_WP
h = 1.0_WP
! Example of calling rk1 with the "type-bound procedure"
! which evaluates an interpolator
call rk4(x, t, h, interp % eval )
print *, x
! Example of calling rk1 with analytical function
call rk4(X, t, h, f )
print *, x
contains
subroutine rk4(x, t, h, f)
! Makes one step with 4th-order Runge-Kutta.
! Calculates next position using timestep h.
implicit none
real(WP), intent(inout), dimension(:) :: x
real(WP), intent(inout) :: t
real(WP), intent(in) :: h
interface
function f(x, t) result(val)
import WP
real(WP), dimension(:), intent(in) :: x
real(WP), intent(in) :: t
real(WP), dimension(size(x)) :: val
end function
end interface
! Local variables
real(WP), dimension(size(x)) :: k1, k2, k3, k4
! Evaluations of f(x, t)
k1 = f(x, t)
k2 = f(x + k1*h/2, t + h/2)
k3 = f(x + k2*h/2, t + h/2)
k4 = f(x + k3, t + h)
! Next position
x = x + h*(k1 + 2*k2 + 2*k3 + k4)/6
t = t + h
end subroutine
pure function f(x, t) result(val)
implicit none
real(WP), dimension(:), intent(in) :: x
real(WP), intent(in) :: t
real(WP), dimension(size(x)) :: val
val = -0.1_WP*x
end function
end program
