Fortran "global" variable vs Matlab function handle

Question

Given a linear system of differential equations, represented in state space form,

dy = A*y + b

In Matlab, I can define the function

function [dy] = fun(t,y,A,b)
    dy = A*y+b;
end

And if I want to do a time integration, I can write in the main program:

n = size(M,1);
A = [zeros(n,n) eye(n,n) ; -M\K -M\C];
b = [zeros(n,1) M\F];
odefun = @(t,y) fun(t,y,A,b)

[TOUT,YOUT] = ode15s(odefun,tspan,y0);

provided I have previously defined the M, C and K matrices and the force vector F, the time span tspan and initial conditions y0.

Now I want this in Fortran. To do that, I found an ode solver online. This is Shampine and Gordon ODE solver (presumably the same solver utilized in MATLAB), which integrates differential equations of the form

y' = f(t,y)

Exactly like MATLAB! If I define a subroutine like (the online example)

subroutine fun( t, y, yp )
    implicit none
    double precision :: t
    double precision :: y(2)
    double precision :: yp(2)

    yp(1) = y(2)
    yp(2) = -y(1)

    return
end subroutine fun

I can integrate this system of equations.

subroutine test01()
    implicit none

    integer :: neqn = 2
    double precision :: abserr
    external fun
    integer :: i
    integer :: iflag
    integer :: iwork(5)
    double precision :: pi = 3.141592653589793D+00
    double precision :: relerr
    integer, parameter :: step_num = 12
    double precision :: t
    double precision :: tout
    double precision :: work(100+21*neqn)
    double precision :: y(neqn)

    abserr = 0.00001D+00
    relerr = 0.00001D+00

    iflag = 1

    t = 0.0D+00
    y(1) = 1.0D+00
    y(2) = 0.0D+00

    write ( *, '(2x,f8.4,2x,2g14.6)' ) t, y(1:neqn)

    do i = 1, step_num

    tout = dble(i)*2.0d0 * pi / dble(step_num)
    call ode ( f01, neqn, y, t, tout, relerr, abserr, iflag, work, iwork )

    if ( iflag /= 2 ) then
      write ( *, '(a)' ) 'TEST01 - Fatal error!'
      write ( *, '(a,i8)' ) '  ODE returned IFLAG = ', iflag
      exit
    end if

    write ( *, '(2x,f8.4,2x,2g14.6)' ) t, y(1:neqn)

    end do

    return
end subroutine test01

This works fine and I can test the ode solver on the main program:

program main
    call test01()
end program main

Here is the problem. Now I want to do a time integration of the equations in state space form. The problem: how to pass the variables A and b to the function fun to be integrated?

subroutine fun( t, y, yp )
    implicit none
    double precision :: t
    double precision :: y(:)
    double precision :: yp(:)

    yp = matmul(A,y) + b

    return
end subroutine fun

I cannot change the number of arguments in the function definition, because the ode integrator requires a fun(t,y,yp) and nothing else.

Is there something in Fortran like MATLAB, where I can define a function fun(t,y,yp,A,b) and a second function fun2 = @(t,x) fun(t,y,yp,A,b) and integrate the second function?

Or tell Fortran A and b are sort of global variables? (I though using the common blocks, but it does not accept allocatables)

Fortran "global" variable vs Matlab function handle

0 Answers0