CODEWITHSUNDEEP
arraysfortran
Thales
Thales

Reputation: 1316

Array assignment erases previous values on array

I'm running the following code, that is the implementation of a Runge-Kutta method to solve a system of differential equations. The main code just calls the rk subroutine, which is the implementation itself, and myfun is just an example to test the code.

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

    call rk4(t, y, myfun, t0, tf, y0, h)

    do i=0,size(t)
        print *, t(i), y(:,i)
    end do

contains

pure function myfun(t,y) result(dy)
    ! input variables
    double precision, intent(in) :: t, y(:)

    ! output variables
    double precision :: dy(size(y))

    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

and the subroutine is inside a 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,y0) result(dy)
            double precision, intent(in) :: t, y0(:)
            double precision :: dy(size(y))
        end function
    end interface

    ! output variables
    double precision, allocatable :: t(:), y(:,:)

    ! Variáveis auxiliares
    integer :: i, m, NN
    double precision, allocatable :: k1(:), k2(:), k3(:), k4(:)

    m = size(y0)
    allocate(k1(m),k2(m),k3(m),k4(m))

    NN = ceiling((tf-t0)/h)
    if (.not. allocated(y)) then
        allocate(y(m,0:NN))
    else
        deallocate(y)
        allocate(y(m,0:NN))
    end if

    if (.not. allocated(t)) then
        allocate(t(0:NN))
    else
        deallocate(t)
        allocate(t(0:NN))
    end if

    t(0) = t0
    y(:,0) = y0


    do i=1,NN
        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)
    return

end subroutine rk4

end module ivp_odes

The problem here is that the assignment in the rk subroutine

y(:,i) = y(:,i-1) + (k1(:) + 2*k2(:) + 2*k3(:) + k4(:))/6

is erasing the previous values calculated. In the i-th iteration of the do-loop, it erases the previous values of the array y and assigns just the i-th column of the array y, so when the subroutine ends, y has only the last value saved. Since Fortran has implemented element-wise operations and assignments to arrays, I think the code this is easier to read and probably runs faster than doing assignments to each element in a loop. So, why is it not working? What am I missing in the assignment here? Shouldn't it just change the values in the i-th row, instead of also erasing the rest of the array?

Upvotes: 1

Views: 84

Answers (1)

Alexander Vogt
Alexander Vogt

Reputation: 18098

This is a typical case of accessing an array out of its bounds. You can find these errors easily using the appropriate compiler flags. With gfortran, this would be -fbounds-check.

With such checks you will find the error to be an erroneous size of the function result in the interface block - dy should have the same length as y0 (the one-dimensional dummy argument of f), and not y: 

    interface
        pure function f(t,y0) result(dy)
            double precision, intent(in) :: t, y0(:)
            double precision :: dy(size(y0))
        end function
    end interface

Additionally, although not related to your particular error, you started indexing of t and the second dimension of y with zero. So you need to adjust the loop in the main program run to size(t)-1 only, or use ubound(t). Otherwise you will, again, exceed the boundaries of the arrays.

Upvotes: 4

Related Questions