I'm trying to implement a runge-kutta for a system of differential equations. I'm having some trouble with the subroutine. I would like to make it a general code, to pass any system of equations to solve.
The main code of the program is:
program main
use ivp_odes
implicit none
double precision, allocatable :: t(:), y(:,:)
double precision :: t0, tf, y0(2), h
integer :: i
t0 = 0d0
tf = 0.5d0
y0 = [0d0, 0d0]
h = 0.1d0
do i=lbound(t,1),ubound(t,1)
print *, t(i), y(1,i), y(2,i)
end do
contains
function myfun(t,y) result(dy)
! input variables
double precision, intent(in) :: t, y(2)
! output variables
double precision :: dy(2)
dy(1) = -4*y(1) + 3*y(2) + 6
dy(2) = -2.4*y(1) + 1.6*y(2) + 3.6
end function myfun
end program main
The myfun is an example function, which I have the analytical solution, so I can compare my implementation with the correct response of the system. This myfun function gets a variable t and a one-dimensional array of two elements y and returns its derivatives, so I can numerically integrate this.
The runge-kutta algorithm I implement in a separate module:
module ivp_odes
implicit none
contains
subroutine rk4(t, y, f, t0, tf, y0, h)
! input variables
double precision, intent(in) :: t0, tf, y0(1:)
double precision, intent(in) :: h
interface
pure function f(t,y) result(dy)
double precision, intent(in) :: t, y(:)
double precision :: dy(size(y))
end function
end interface
! output variables
double precision, allocatable :: t(:), y(:,:)
! auxiliar variables
integer :: i, m, N
double precision :: hh
double precision, allocatable :: k1(:), k2(:), k3(:), k4(:)
N = ceiling((tf-t0)/h)
m = size(y0)
allocate(k1(m),k2(m),k3(m),k4(m))
if (.not. allocated(y)) then
allocate(y(m,0:N))
else
deallocate(y)
allocate(y(m,0:N))
end if
if (.not. allocated(t)) then
allocate(t(0:N))
else
deallocate(t)
allocate(t(0:N))
end if
t(0) = t0
y(:,0) = y0
do i=1,N
k1(:) = h * f(t(i-1) , y(:,i-1) )
k2(:) = h * f(t(i-1)+h/2 , y(:,i-1)+k1(:)/2)
k3(:) = h * f(t(i-1)+h/2 , y(:,i-1)+k2(:)/2)
k4(:) = h * f(t(i-1)+h , y(:,i-1)+k3(:) )
y(:,i) = y(:,i-1) + (k1(:) + 2*k2(:) + 2*k3(:) + k4(:))/6
t(i) = t(i-1) + h
end do
deallocate(k1,k2,k3,k4)
end subroutine rk4
end module ivp_odes
The idea of the rk subroutine is that t is the time vector (one dimensional array) and y is the response vector (two dimensional array, each row corresponds to a variable and each column corresponds to the value of the variable in the correspondent time).
I can't find a way to make it work, since I get errors when compiling the code, all related to the function myfun when calling the subroutine rk in the main code.
call rk4(t, y, myfun, t0, tf, y0, h)
1
Error: Interface mismatch in dummy procedure 'f' at (1): Shape mismatch in dimension 1 of function result
I tried to look for an answer to this problem, and already read this previously asked questions:
Interface mismatch - higher order functions
Fortran array cannot be returned in function: not a DUMMY variable
Function in fortran, passing array in, receiving array out
some of which already helped me solve some errors. However, what I'm trying to do is pretty much the same thing as the first question, but it simply doesn't work.
I'm using the GNU Fortran Compiler, and writing the code in Code::Blocks 16.01. Can't find a way to solve it.
Here a few things I tried that hasn't worked:
- change the variable in the interface:
interface
pure function f(t,yy) result(dy)
double precision, intent(in) :: t, yy(:)
double precision :: dy(size(y))
end function
end interface
- not declare the result of the funciont in interface
interface
pure function f(t,y)
double precision, intent(in) :: t, y(:)
end function
end interface
it produces the error:
call rk4(t, y, myfun, t0, tf, y0, h)
1
Error: Interface mismatch in dummy procedure 'f' at (1): Type mismatch in function result (REAL(4)/REAL(8))
- I though I could use the procedure statement, but I think that to do that, I should declare my function myfun in the same module of the subroutine, right? But the idea to have a general code is that I don't need to open this module after finished, so anytime I want to use my rk4 integrator, I just call it, no noeed to specify the function in the same module.
