Array-Zuweisung löscht vorherige Werte auf Array

Question

Ich führe den folgenden Code, das ist die Implementierung einer Runge-Kutta-Methode, um ein System von Differentialgleichungen zu lösen. Der Hauptcode ruft nur die Unterroutine rk auf, die die Implementierung selbst ist, und myfun ist nur ein Beispiel, um den Code zu testen.

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 

und das Unterprogramm ist in einem Modul:

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 

Das Problem hierbei ist, dass die Zuordnung in der rk Subroutine

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

berechnet die vorherigen Werte zu löschen. In der i-ten Iteration der do-Schleife löscht er die vorherigen Werte des Arrays y und weist nur die i-te Spalte des Arrays y zu. Wenn das Unterprogramm endet, hat y nur den letzten gespeicherten Wert. Da Fortran elementare Operationen und Zuweisungen zu Arrays implementiert hat, denke ich, dass der Code leichter zu lesen ist und wahrscheinlich schneller läuft als Zuweisungen für jedes Element in einer Schleife. Also, warum funktioniert es nicht? Was fehlt mir in der Aufgabe hier? Sollte es nicht nur die Werte in der i-ten Zeile ändern, anstatt auch den Rest des Arrays zu löschen?

Answer 1

Dies ist ein typischer Fall des Zugriffs auf ein Array außerhalb seiner Grenzen. Sie können diese Fehler leicht mit den entsprechenden Compilerflags finden. Mit gfortran wäre dies -fbounds-check.

solche Kontrollen mit der Fehler bestimmt wird eine fehlerhafte Größe des Funktionsergebnisses im Interface-Block sein - dy sollte die gleiche Länge wie y0 haben (die eindimensionale Scheinargument von f) und nicht y:

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

Zusätzlich, wenn auch nicht auf Ihre Fehler im Zusammenhang begannen Sie die Indizierung von t und die zweite Dimension von y mit Null. Sie müssen also die Schleife im Hauptprogrammlauf nur auf size(t)-1 einstellen oder ubound(t) verwenden. Sonst überschreiten Sie wiederum die Grenzen der Arrays.