配列の割り当てによって配列の前の値が消去される

Question

次のコードを実行しています。これは、ルンゲクッタ法を実装して微分方程式を解く方法です。 メインコードは実装自体であるrkサブルーチンを呼び出し、myfunはコードをテストする単なる例に過ぎません。

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 

、サブルーチン、モジュール内部にある：

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 

ここでの問題は、rkサブルーチン

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

における割り当てが計算前の値を消去することです。 doループのi番目の反復では、配列yの以前の値が消去され、配列yのi番目の列が割り当てられるので、サブルーチンが終了すると、yには最後の値のみが保存されます。 Fortranは要素ごとの操作と配列への代入を実装しているので、コードを読みやすく、おそらくループ内の各要素に代入するよりも高速に動作すると思います。だから、なぜそれは動作していないのですか？ここで私は何が欠けていますか？配列の残りの部分を消去するのではなく、i行目の値を変更するだけでいいですか？

Answer 1

これは、境界からアレイにアクセスする典型的なケースです。これらのエラーは、適切なコンパイラフラグを使用して簡単に見つけることができます。 gfortranでは、これは-fbounds-checkになります。もしインターフェースブロックにおける関数の結果の誤ったサイズにエラーがありますこのようなチェックを

からdyがy0と同じ長さ（fの一次元仮引数）を有するべきではなく、y：

あなたの特定のエラーに関連していないが

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

また、あなたはtのインデックスとゼロとyの第2次元を開始しました。したがって、メインプログラムのループをsize(t)-1にのみ調整するか、ubound(t)を使用する必要があります。それ以外の場合は、配列の境界を超えます。