There’s no silver bullet. As always it depends on what are your requirements and what are your preferences. In Modern Fortran in Practice by @Arjen Markus provides a table comparing some of the different callback options:
In a way I think that callbacks were easier in the era of punch-cards and serial programs. You would just place your callback at the top/bottom of your stack of punch cards and declare it as external in the driver program. Still, the callback interface of an early code from Hindmarsh (epsode.f) remains just as valid today as it was in 1974:
c the user must furnish the following double precision subroutines
c and declare them external in his calling program:
c diffun(n,t,y,ydot) computes the function ydot = f(y,t),
c the right hand side of the ordinary
c differential equation system, where y
c and ydot are vectors of length n.
If you don’t like this interface you can always use an internal routine and host association to create an adaptor routine with your interface of choice. This comes at a cost of one indirection (function call), which is probably irrelevant in the grand scheme of things.
The same indirection occurs in MATLAB, but due to the dynamic typing is less intrusive. The ODE solvers in MATLAB generally expect a function of the form:
function dydt = odefun(t,y)
dydt = zeros(2,1); # size declared explicitly, or assumed from input y
dydt(1) = ...;
dydt(2) = ...;
end
If you require a parameterized function, you create a function handle to an anonymous function (here’s the indirection):
# lorenz.m
function dydt = lorenz(t,y,sigma,beta,rho)
dydt = zeros(3,1);
dydt(1) = sigma*(y(2) - y(1))
dydt(2) = y(1)*(rho - y(3)) - y(2)
dydt(3) = y(1)*y(2) - beta*y(3)
end
# solve_lorenz.m
sigma = 10;
beta = 8/3;
rho = 28;
f = @(t,y) lorenz(t,y,sigma,beta,rho);
y0 = [1,2,1];
[tspan,y] = ode45(f,[0 100],y0);
figure
plot3(y(:,1),y(:,2),y(:,3))
It’s the same thing in Python and R. With Julia it’s slightly different in that their JIT compiler can perform some tricks to inline the callback directly into the solver. This can lead to time savings for small systems of ODEs. With Fortran you could do something similar “manually” with a preprocessor such as fypp (I don’t recommend doing this). For large (& small) systems of ODE’s I don’t think it’s worth the trouble. Just for fun, I tried to inline the function call with a solver I’m familiar with:
! lorenz.fypp
module lorenz
implicit none
public
integer, parameter :: wp = kind(1.0e0)
type :: lorenz_params
real(wp) :: sigma, beta, rho
end type
contains
#! N.B.: The argument names y, dydt, and params are fixed!
#:def lorenz()
dydt(1) = params%sigma*(y(2) - y(1))
dydt(2) = y(1)*(params%rho - y(3)) - y(2)
dydt(3) = y(1)*y(2) - params%beta*y(3)
#:enddef
#! The name jac is also fixed!
#:def lorenz_jac()
jac(1, 1) = -params%sigma
jac(2, 1) = params%rho - y(3)
jac(3, 1) = y(2)
jac(1, 2) = params%sigma
jac(2, 2) = -1
jac(3, 2) = y(1)
jac(1, 3) = 0
jac(2, 3) = -y(1)
jac(3, 3) = -params%beta
#:enddef
! subroutine solve_lorenz( ... )
#:include "stiff3_solver.fypp"
$:solve("lorenz",3,lorenz,lorenz_jac,"lorenz_params")
end module
! lorenz_main.f90
program main
use lorenz
implicit none
integer, parameter :: n = 3
real(wp) :: y(n), w(n)
real(wp) :: h0, eps, x0, x1
integer, parameter :: nprint = 1
type(lorenz_params) :: p
p = lorenz_params( &
sigma = 10.0_wp, &
beta = 8.0_wp/3.0_wp, &
rho = 28.0_wp)
! initial value
y = [1.0_wp, 1.0_wp, 1.0_wp]
! initial step size
h0 = 0.001_wp
! tolerance parameters
w = 1
eps = 1.e-4_wp
! time interval
x0 = 0.0_wp
x1 = 100.0_wp
call output(x0,y,0,0.0_wp)
call solve_lorenz(output,nprint,x0,x1,h0,eps,w,y,p)
contains
! Callback routine for output
subroutine output(x,y,iha,qa)
real(wp), intent(in) :: x
real(wp), intent(in) :: y(:)
integer, intent(in) :: iha
real(wp), intent(in) :: qa
write(*,'(4(E19.12,2X),I4,2X,E19.12)') &
x, y(1), y(2), y(3), iha, qa
end subroutine
end program
The final output:
In case anyone is interested here is the full example: lorenz.tar.gz - Google Drive (4.2 KiB).