Consider the following: I would like to solve a differential equation in Fortran (with a Runge-Kutta routine), e.g. x'=-a*x
, where a
is a constant. I then want to solve this for different a
. I now would like to write the code such that it can be applied to other problems. But in general I do not know what additional input my differential equation may need (and what size and type they have). So my question is if it is possible in Fortran to pass all additional parameters to my Runge-Kutta method (at best with a single property) so that it passes it on to my differential equation so that it can be used there? With a known set of parameters this is easy but what about an unknown set?
My approach so far is this:
PROGRAM main
USE RKModule
IMPLICIT NONE
INTEGER :: l
REAL*8, PARAMETER :: num_t = 1e3
REAL*8, PARAMETER :: min_t = 0.0d0
REAL*8, PARAMETER :: max_t = 10.0d0
REAL*8, DIMENSION(1), PARAMETER :: a = 1.8d0
REAL*8, DIMENSION(num_t) :: dGrid_t, grid_t
COMPLEX*16, DIMENSION(num_t) :: sol
COMPLEX*16, DIMENSION(1) :: x_in, x_out
DO l = 1, num_t
dGrid_t(l) = (max_t - min_t) / DBLE(num_t)
grid_t(l) = min_t + DBLE(l-0.5d0)*dGrid_t(l)
END DO
x_in(1) = dcmplx(1.0d0, 0.0d0)
x_out(1) = dcmplx(0.0d0, 0.0d0)
DO l = 1, num_t
sol(l) = x_in(1)
CALL rk45(x_out, DGL, x_in, grid_t(l), dGrid_t(l), a)
x_in = x_out
END DO
OPEN(UNIT = 1000, FILE = 'sol.dat', ACTION = 'WRITE')
WRITE(1000, FMT = '(e14.7e3, 1x, e14.7e3)', ADVANCE = "NO") sol
CLOSE(1000)
CONTAINS
SUBROUTINE DGL(t,f,df,vararg_in)
REAL*8, INTENT(IN) :: t
COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in
REAL*8, DIMENSION(1) :: b
b = 1.0d0
IF (PRESENT(vararg_in)) THEN
b = vararg_in
END IF
df = -b*f
END SUBROUTINE DGL
END PROGRAM main
And the Runge Kutta Routine:
MODULE RKModule
IMPLICIT NONE
CONTAINS
SUBROUTINE rk45(x_out, RHS, x_in, t, dt, vararg_in)
INTERFACE
SUBROUTINE RHS(t,f,df,a)
REAL*8, INTENT(IN) :: t
COMPLEX*16, DIMENSION(:), INTENT(IN) :: f
COMPLEX*16, DIMENSION(:), INTENT(OUT) :: df
REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: a
END SUBROUTINE
END INTERFACE
REAL*8, INTENT(IN) :: t, dt
COMPLEX*16, DIMENSION(:), INTENT(IN) :: x_in
COMPLEX*16, DIMENSION(:), INTENT(OUT) :: x_out
REAL*8, DIMENSION(:), INTENT(IN), OPTIONAL :: vararg_in
REAL*8 :: tHlp
COMPLEX*16, DIMENSION(SIZE(x_in)) :: k1 ,k2, k3, k4, dx, x_hlp
REAL*8, DIMENSION(:), ALLOCATABLE :: vararg_in_copy
IF (PRESENT(vararg_in)) THEN
ALLOCATE(vararg_in_copy(SIZE(vararg_in)))
vararg_in_copy = vararg_in
tHlp = t
x_hlp = x_in
CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
k1 = dt * dx
tHlp = t + 0.5d0 * dt
x_hlp = x_in + 0.5d0 * k1
CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
k2 = dt * dx
x_hlp = x_in + 0.5d0 * k2
CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
k3 = dt * dx
tHlp = t + dt
x_hlp = x_in + k3
CALL RHS(tHlp, x_hlp, dx, vararg_in_copy)
k4 = dt * dx
x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0
DEALLOCATE(vararg_in_copy)
ELSE
tHlp = t
x_hlp = x_in
CALL RHS(tHlp, x_hlp, dx)
k1 = dt * dx
tHlp = t + 0.5d0 * dt
x_hlp = x_in + 0.5d0 * k1
CALL RHS(tHlp, x_hlp, dx)
k2 = dt * dx
x_hlp = x_in + 0.5d0 * k2
CALL RHS(tHlp, x_hlp, dx)
k3 = dt * dx
tHlp = t + dt
x_hlp = x_in + k3
CALL RHS(tHlp, x_hlp, dx)
k4 = dt * dx
x_out = x_in + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0
END IF
END SUBROUTINE rk45
END MODULE RKModule
ISO_C_BINDING
or just move to another language like C/C++ that supports variadic arguments if you have to do it this way. $\endgroup$ – Doug Lipinski May 24 '14 at 13:06