阵列分配会擦除阵列上的先前值

我正在运行以下代码，即Runge-Kutta方法的实现来解决微分方程组。主代码只是调用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行的值，而不是擦除数组的其余部分？

来源

2016-05-22 Thales

您如何确定'y'中除最后一列以外的每一列都已被擦除？你确定你通过了2D数组吗？ –

这是一个访问数组超出范围的典型情况。您可以使用适当的编译器标志轻松找到这些错误。用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零的第二个维度。因此，您需要将主程序中的循环调整为size(t)-1，或使用ubound(t)。否则，你会再次超出数组的边界。

来源

2016-05-22 16:00:51

阵列分配会擦除阵列上的先前值

回答

相关问题