Fortran | Passing functions as arguments in other functions

Question

I have a simple system of 1st-order differential equations to be solved in Matlab. It looks like the following:

function passing
y0 = [0; 0];
t = 0:0.05:1;
y = myprocedure(@myfunc,0.05,t,y0);
myans = y'
end

function f = myfunc(t,y)
f = [-y(1) + y(2);
     -0.8*t + y(2)];
end

function y=myprocedure(f,h,t,y0)
n = length(t);
y = zeros(length(y0),n);
y(:,1) = y0;
for i=1:n-1
    k1 = f(t(i),y(:,i));
    k2 = f(t(i)+h/2, y(:,i)+h/2*k1);
    k3 = f(t(i)+h/2, y(:,i)+h/2*k2);
    k4 = f(t(i)+h, y(:,i)+h*k3);
    y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4);
end
end

So, after running the passing.m file, I have the result:

>> passing

myans =

     0         0
-0.0000   -0.0010
-0.0001   -0.0041
-0.0005   -0.0095
-0.0011   -0.0171
-0.0021   -0.0272
-0.0036   -0.0399
-0.0058   -0.0553
-0.0086   -0.0735
-0.0123   -0.0946
-0.0169   -0.1190
-0.0225   -0.1466
-0.0293   -0.1777
-0.0374   -0.2124
-0.0469   -0.2510
-0.0579   -0.2936
-0.0705   -0.3404
-0.0849   -0.3917
-0.1012   -0.4477
-0.1196   -0.5086
-0.1402   -0.5746

Now, I am trying converting the passing.m code to Fortran 90. I did try lots of things, but I am still lost -- particularly in doing array initializations and passing functions to other functions.

Any help is appreciated. Thanks!

== EDIT:

I succeeded replicating the above Matlab code in Fortran 90. I still don't know how to construct pointers, but the "interface" statement seems to work. I would appreciate it if you look at the code and provide some input. Thanks.

module odemodule
implicit none
contains
function odef(a,b)
 implicit none
 real::a
 real::b(2)
 real::odef(2)
 odef(1) = -b(1) + b(2)
 odef(2) = -0.8*a + b(2)
end function odef
subroutine rk4(f,h,t,y0,y)
 implicit none
 real,dimension(2,21)::y
 real,dimension(2)::k1,k2,k3,k4
 real::h
 real::t(21)
 real::y0(2)
 integer::i
 interface
  function f(a,b)
   real::a
   real::b(2), f(2)
  end function f
 end interface
 y(1,1)=y0(1)
 y(2,1)=y0(2)
 do i=1,20     
 k1 = f(t(i),y(:,i))
 k2 = f(t(i)+h/2, y(:,i)+h/2*k1)
 k3 = f(t(i)+h/2, y(:,i)+h/2*k2)
 k4 = f(t(i)+h, y(:,i)+h*k3)
 y(:,i+1) = y(:,i)+h/6*(k1+2*k2+2*k3+k4)
 end do
end subroutine rk4
end module odemodule

program passing
  use odemodule
implicit none
real,parameter::dt=0.05
real::yinit(2)
real::y(2,21)
integer::i
real::t(21)
t(1)=0.0
do i=2,21
t(i)=t(i-1)+dt
end do
yinit(1)=0
yinit(2)=0
call rk4(odef,dt,t,yinit,y)
do i=1,21
  write(*,100) y(1,i), y(2,i)
enddo
100 format(f7.4,3x,f7.4)
end program passing

Answer 1

The best way in Fortran to select functions to pass to other procedures (subroutine or functions) is to use procedure pointers. That way you can write the procedure in a general manner and invoke it with a specific function. Here are some examples: How to alias a function name in Fortran and Function pointer arrays in Fortran

There are many ways to initialize arrrays. For example, you can initialize all elements to the same value, e.g.,

array = 0.0

See http://en.wikipedia.org/wiki/Fortran_95_language_features

Answer 2

There is another way, I want to make it in a single file.

program main
    implicit none
    integer      , parameter       :: wp = 8
    integer      , parameter       ::  n = 2       ! number of equation
    real(kind=wp), parameter       ::  h = 0.05_wp ! time step
    integer      , parameter       ::  m = ceiling(1.0/h) + 1 ! total number of time step
    real(kind=wp), dimension(m)    ::  t
    real(kind=wp), dimension(n)    ::  y0
    real(kind=wp), dimension(n, m) ::  y
    integer :: i
    
    t(1:m) = [ (i*h ,i = 0,m-1) ] ! if h=0.05 m=21 ( t = 0:0.05:1 )
    y0(:)  = [  0.0_wp,  0.0_wp ]
    y = ode_rk4(f, h, t, y0)

    write(*,'(a,10x,a4,20x,a,30x,a)' ) '#', 'time', 'x', 'v'
    do i = 1, m
       write(*,*) t(i), y(:,i)
    end do

contains
    function ode_rk4(f, h, t, y0) result(y)
        implicit none
        real(kind=wp), intent(in), dimension(:)     :: y0
        real(kind=wp), intent(in), dimension(:)     :: t
        real(kind=wp), intent(in)                   :: h
        real(kind=wp), dimension(size(y0), size(t)) :: y

        interface
            pure function f(t, x) result(y)
                implicit none
                real(kind=8)              , intent(in) :: t
                real(kind=8), dimension(2), intent(in) :: x
                real(kind=8), dimension(2)             :: y
            end function f
        end interface

        real(kind=wp), dimension( size(y0) ) :: k1, k2, k3, k4

        y(:, :) = 0.0_wp
        y(:, 1) = y0

        do i = 1, size(t) - 1
            k1(:) = f( t(i), y(:,i) )
            k2(:) = f( t(i) + h*0.5, y(:,i) + h*k1(:)*0.5 )
            k3(:) = f( t(i) + h*0.5, y(:,i) + h*k2(:)*0.5 )
            k4(:) = f( t(i) + h, y(:,i) + h*k2(:) )

            y(:, i+1) = y(:, i) + h/6.0 * ( k1 + 2.0*k2 + 2.0*k3 + k4 )

        end do
    end function ode_rk4

    
    pure function f(t, x) result(y)
        implicit none
        real(kind=wp)              , intent(in) :: t
        real(kind=wp), dimension(2), intent(in) :: x
        real(kind=wp), dimension(2)             :: y

        y(1:2) = [ -x(1) + x(2), &
                  -0.8*t + x(2) ]

    end function f

end program main

But the problem with the block interface of the function f(t,x).

interface
    pure function f(t, x) result(y)
        implicit none
        real(kind=8)              , intent(in) :: t
        real(kind=8), dimension(2), intent(in) :: x
        real(kind=8), dimension(2)             :: y
    end function f
end interface

I want it all declaration variable depending on the precision variable wp not on hard coded values, I could use module, for exemple create file kind_mod.f90

module kind
    implicit none
    integer, private, parameter :: sp = 4
    integer, private, parameter :: dp = 8
    integer, private, parameter :: qp = 16
    integer, public , parameter :: wp = dp
end module kind

and then

interface
    pure function f(t, x) result(y)
        use kind, only: wp
        implicit none
        real(kind=wp)              , intent(in) :: t
        real(kind=wp), dimension(2), intent(in) :: x
        real(kind=wp), dimension(2)             :: y
    end function f
end interface

The code could be improve more, but i hope that was enough to demonstrate how you make function as an argument of another function by using interface, you can also execute it with any precision you want.

#          time                    x                              v
   0.0000000000000000        0.0000000000000000        0.0000000000000000     
   5.0000000000000003E-002  -1.6666666915019357E-005  -1.0166666818161806E-003
  0.10000000000000001       -1.3337500198744241E-004  -4.1362920755244432E-003
  0.15000000000000002       -4.5043763692761095E-004  -9.4666966267490885E-003
  0.20000000000000001       -1.0686681413439925E-003  -1.7121228825211946E-002
  0.25000000000000000       -2.0896331089569026E-003  -2.7219048632159955E-002
  0.30000000000000004       -3.6159061266336280E-003  -3.9885425439353160E-002
  0.35000000000000003       -5.7513242612989464E-003  -5.5252051304658302E-002
  0.40000000000000002       -8.6012477060707013E-003  -7.3457370247492326E-002
  0.45000000000000001       -1.2272823234869468E-002  -9.4646924427517695E-002
  0.50000000000000000       -1.6875252124274188E-002 -0.11897371807220791     
  0.55000000000000004       -2.2520063212565732E-002 -0.14659860006328257     
  0.60000000000000009       -2.9321391778745664E-002 -0.17769066613866774     
  0.65000000000000002       -3.7396264938870084E-002 -0.21242768171568613     
  0.70000000000000007       -4.6864894273334769E-002 -0.25099652639274433     
  0.75000000000000000       -5.7850976416828570E-002 -0.29359366124099234     
  0.80000000000000004       -7.0482002362582521E-002 -0.34042562005441557     
  0.85000000000000009       -8.4889576254331925E-002 -0.39170952578672846     
  0.90000000000000002      -0.10120974446313261      -0.44767363346641825     
  0.95000000000000007      -0.11958333577188944      -0.50855790094749531     
   1.0000000000000000      -0.14015631351823021      -0.57461458892310990