虚拟程序中的接口不匹配：形状不匹配

Question

我正在尝试为微分方程系统实现 runge-kutta。 我的子程序有点问题。我想让它成为一个通用代码，以传递任何方程组来求解。

程序的主要代码是：

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

myfun 是一个示例函数，我有解析解，所以我可以将我的实现与系统的正确响应进行比较。这个myfun 函数得到一个变量t 和一个由两个元素组成的一维数组y 并返回它的导数，所以我可以对它进行数值积分。

我在一个单独的模块中实现的 runge-kutta 算法：

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

rk子程序的思想是t是时间向量（一维数组），y是响应向量（二维数组，每行对应一个变量，每列对应对应时间变量的值）。

我找不到让它工作的方法，因为我在编译代码时遇到错误，所有这些都与在主代码中调用子程序rk 时的函数myfun 有关。

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

我试图寻找这个问题的答案，并且已经阅读了之前提出的问题：

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

其中一些已经帮助我解决了一些错误。但是，我想要做的与第一个问题几乎相同，但它根本不起作用。

我正在使用 GNU Fortran 编译器，并在 Code::Blocks 16.01 中编写代码。 找不到解决办法。

以下是我尝试过的一些没有奏效的方法：

• 改变界面中的变量：

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

• 不在接口中声明函数的结果

interface
    pure function f(t,y)
        double precision, intent(in) :: t, y(:)
    end function
end interface

它会产生错误：

    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))

• 虽然我可以使用过程语句，但我认为要做到这一点，我应该在子例程的同一模块中声明我的函数 myfun ，对吗？但是有一个通用代码的想法是我完成后不需要打开这个模块，所以我什么时候想使用我的rk4积分器，我只是调用它，不需要在同一个模块中指定函数。李>

Answer 1

您传递的函数myfun 和子程序rk4 预期的函数f 不是同一类型：

• f 是 pure，而 myfun 不是。
• 在f 中，y 是延迟形式，而在myfun 中是显式形式（对dy 产生影响）。

这实际上是编译器所抱怨的 :) 一旦你解决了这些问题，代码就可以正常编译。