0

我有兴趣编写一个函数，该函数将要使用的模块名称作为其输入之一。例如，我编写了一个用于求解 ODE 系统的 Runge Kutta 4 阶积分器。我更愿意编写实际的 Runge Kutta 积分函数，以便为定义 ODE 系统的函数提供输入。

这是一个示例模块，其中包含定义 ODE 系统的函数

MODULE DIFF_EQ
! Description: This module contains the function that defines the
!              system of ODEs to be solved.

CONTAINS

    FUNCTION YDOT(t, y, PP) RESULT(yd)
    ! Description: This function defines the system of ODEs to be
    !              solved.
    !
    ! Inputs:   t - current independent variable value (real, scalar)
    !           y - current dependent variable value (real, array)
    !          PP - passed parameters/constants (real, array)
    !
    ! Outputs: yd - current value of ODE (real, array)

    IMPLICIT NONE

    REAL*8, INTENT(IN)               :: t
    REAL*8, DIMENSION(2), INTENT(IN) :: y
    REAL*8, DIMENSION(3), INTENT(IN) :: PP
    REAL*8, DIMENSION(2)             :: yd

    yd = [y(2), &
          -PP(1)*SIN(y(1)) + SIN(PP(2)*t) + PP(3)]

    END FUNCTION YDOT

END MODULE DIFF_EQ

这是 Runge Kutta 集成模块

MODULE RUNGE_KUTTA
! Description: This module contains the function that implements the Runge Kutta 4th 
!              order integration scheme.

CONTAINS

    FUNCTION RK4(Neq, tspan, y0, dt, DEQ, PP) RESULT(t_y)
    ! Description: This function implements the Runge Kutta 4th order integration scheme.
    !
    ! Inputs: Neq   - number of equations in system of ODEs (integer, scalar)
    !         tspan - [t0, tF] where t0 is start time and tF is end time (real, array - 2)
    !         y0    - [y1, y2, ..., yNeq] @ t0 (real, array - Neq)
    !         dt    - time step size such that t1 = t0 + dt (real, scalar)
    !         PP    - passed parameters/constants (real, array - variable)
    !
    ! Outputs: t_y  - time and solution (real, matrix - n x Neq + 1)

    USE DIFF_EQ

    IMPLICIT NONE

    INTEGER, INTENT(IN)                  :: Neq
    REAL*8, DIMENSION(2), INTENT(IN)     :: tspan
    REAL*8, DIMENSION(Neq), INTENT(IN)   :: y0
    REAL*8                               :: dt
    REAL*8, EXTERNAL                     :: DEQ
    REAL*8, DIMENSION(:), INTENT(IN)     :: PP
    INTEGER                              :: n, i
    REAL*8, DIMENSION(:,:), ALLOCATABLE  :: t_y
    REAL*8, DIMENSION(Neq)               :: k1, k2, k3, k4

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)
    ALLOCATE(t_y(n, Neq+1))

    IF (MOD(tspan(2) - tspan(1), dt) .LT. 0.000000000000001D0) THEN
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n)]
    ELSE
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n-1), tspan(2)]
    ENDIF

    t_y(1, 2:Neq+1) = y0

    PRINT *, t_y(1, 1:Neq+1)

    DO i = 2, n
        dt = t_y(i, 1) - t_y(i-1, 1)
        k1 = DEQ(t_y(i-1, 1),            t_y(i-1, 2:Neq+1),               PP)
        k2 = DEQ(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k1, PP)
        k3 = DEQ(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k2, PP)
        k4 = DEQ(t_y(i-1, 1) + dt,       t_y(i-1, 2:Neq+1) + dt*k3,       PP)
        t_y(i, 2:Neq+1) = t_y(i-1, 2:Neq+1) + 0.16666666666666667D0*dt*(k1 + 2.0D0*k2 + 2.0D0*k3 + k4)

        PRINT *, t_y(i, 1:Neq+1)
    ENDDO

    END FUNCTION RK4

END MODULE RUNGE_KUTTA

以及接受输入并调用函数的主程序

PROGRAM RK4_TEST

    USE DIFF_EQ
    USE RUNGE_KUTTA

    IMPLICIT NONE

    INTEGER                             :: Neq
    REAL*8, DIMENSION(2)                :: tspan
    REAL*8, DIMENSION(2)                :: y0
    REAL*8                              :: dt
    REAL*8, DIMENSION(3)                :: PP
    INTEGER                             :: n
    REAL*8, DIMENSION(:,:), ALLOCATABLE :: t_y

    Neq   = 2
    tspan = [0.0D0, 1.0D0]
    y0    = [1.0D0, 0.0D0]
    dt    = 0.1D0
    PP    = [1.0D0, 5.0D0, 0.0D0]

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)

    ALLOCATE(t_y(n, Neq+1))

    t_y = RK4(Neq, tspan, y0, dt, YDOT, PP)

END PROGRAM RK4_TEST

使用 gfortran，这是由

gfortran DIFF_EQ.f90 RUNGE_KUTTA.f90 RK4_TEST.f90 -o MAIN.exe

Segmentation fault (core dumped)它编译没有错误，但是当我尝试运行它时收到一条消息。

感谢您提供的任何帮助。

4

1 回答 1

4

您可以将函数传递给求解器。这应该足以使您的求解器具有通用性，能够处理各种功能。让您的求解器将函数作为参数，然后使用您希望求解的特定函数的实际参数调用它。或者您可以使用函数指针调用它，每次调用都有指向您希望解决的函数的指针。在您的示例中，您不传递类似于 的模块DIFF_EQ，而是传递类似于YDOT. 这是函数指针的示例：Fortran 中的函数指针数组

编辑：好的，这是一个代码示例。我已经展示了如何将第二个函数传递给求解器。通过使用接口语句在求解器中声明函数，模块不必在求解器中使用。这些函数可以很容易地位于不同的模块中，并且不需要修改求解器的源代码，只需修改调用程序的源代码。

MODULE DIFF_EQ

   use ISO_FORTRAN_ENV

! Description: This module contains the function that defines the
!              system of ODEs to be solved.

CONTAINS

    FUNCTION YDOT(t, y, PP) RESULT(yd)
    ! Description: This function defines the system of ODEs to be
    !              solved.
    !
    ! Inputs:   t - current independent variable value (real, scalar)
    !           y - current dependent variable value (real, array)
    !          PP - passed parameters/constants (real, array)
    !
    ! Outputs: yd - current value of ODE (real, array)

    IMPLICIT NONE

    real (real64), INTENT(IN) :: t
    real (real64), INTENT(IN) :: y(2)
    real (real64), INTENT(IN) :: PP(3)
    real (real64)             :: yd(2)

    yd = [y(2), -PP(1)*SIN(y(1)) + SIN(PP(2)*t) + PP(3)]

    END FUNCTION YDOT

    FUNCTION YDOT2 (t, y, PP) RESULT(yd)
    ! Description: This function defines the system of ODEs to be
    !              solved.
    !
    ! Inputs:   t - current independent variable value (real, scalar)
    !           y - current dependent variable value (real, array)
    !          PP - passed parameters/constants (real, array)
    !
    ! Outputs: yd - current value of ODE (real, array)

    IMPLICIT NONE

    real (real64), INTENT(IN) :: t
    real (real64), INTENT(IN) :: y(2)
    real (real64), INTENT(IN) :: PP(3)
    real (real64)             :: yd(2)

    yd = [y(2), -PP(1)*cos(y(1)) + cos(PP(2)*t) + PP(3)]

    END FUNCTION YDOT2

END MODULE DIFF_EQ

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

MODULE RUNGE_KUTTA

   use ISO_FORTRAN_ENV
! Description: This module contains the function that implements the Runge Kutta 4th
!              order integration scheme.

CONTAINS

    FUNCTION RK4(DEQ_func, Neq, tspan, y0, dt, PP) RESULT(t_y)
    ! Description: This function implements the Runge Kutta 4th order integration scheme.
    !
    ! Inputs: Neq   - number of equations in system of ODEs (integer, scalar)
    !         tspan - [t0, tF] where t0 is start time and tF is end time (real, array - 2)
    !         y0    - [y1, y2, ..., yNeq] @ t0 (real, array - Neq)
    !         dt    - time step size such that t1 = t0 + dt (real, scalar)
    !         PP    - passed parameters/constants (real, array - variable)
    !
    ! Outputs: t_y  - time and solution (real, matrix - n x Neq + 1)


    IMPLICIT NONE

   interface

      function  DEQ_func ( t, y, pp )  result (yd)

         use ISO_FORTRAN_ENV
         real (real64), INTENT(IN) :: t
         real (real64), dimension (2), INTENT(IN) :: y
         real (real64), dimension (3), INTENT(IN) :: PP
         real (real64), dimension (2)             :: yd

      end function DEQ_func

   end interface

    INTEGER, INTENT(IN)                  :: Neq
    real (real64), DIMENSION(2), INTENT(IN)     :: tspan
    real (real64), INTENT(IN)                   :: y0(Neq)
    real (real64)                               :: dt
    real (real64), INTENT(IN)                   :: PP(:)
    INTEGER                              :: n, i
    real (real64), DIMENSION(:,:), ALLOCATABLE  :: t_y
    real (real64)                               :: k1(Neq), k2(Neq), k3(Neq), k4(Neq)

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)
    ALLOCATE(t_y(n, Neq+1))

    IF (MOD(tspan(2) - tspan(1), dt) .LT. 0.000000000000001D0) THEN
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n)]
    ELSE
        t_y(1:n, 1) = [(tspan(1) + dt*(i-1), i = 1, n-1), tspan(2)]
    ENDIF

    t_y(1, 2:Neq+1) = y0

    PRINT *, t_y(1, 1:Neq+1)
    DO i = 2, n
        dt = t_y(i, 1) - t_y(i-1, 1)
        k1 = DEQ_func(t_y(i-1, 1),              t_y(i-1, 2:Neq+1),                 PP)
        k2 = DEQ_func(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k1, PP)
        k3 = DEQ_func(t_y(i-1, 1) + 0.5D0*dt, t_y(i-1, 2:Neq+1) + 0.5D0*dt*k2, PP)
        k4 = DEQ_func(t_y(i-1, 1) + dt,         t_y(i-1, 2:Neq+1) + dt*k3,         PP)
        t_y(i, 2:Neq+1) = t_y(i-1, 2:Neq+1) + 0.16666666666666667D0*dt*(k1 + 2.0D0*k2 + 2.0D0*k3 + k4)

        PRINT *, t_y(i, 1:Neq+1)
    ENDDO

    END FUNCTION RK4

END MODULE RUNGE_KUTTA

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

PROGRAM RK4_TEST

    USE RUNGE_KUTTA
    use DIFF_EQ

    IMPLICIT NONE

    INTEGER                             :: Neq
    real (real64), DIMENSION(2)                :: tspan
    real (real64), DIMENSION(2)                :: y0
    real (real64)                              :: dt
    real (real64), DIMENSION(3)                :: PP
    INTEGER                             :: n
    real (real64), DIMENSION(:,:), ALLOCATABLE :: t_y

    Neq   = 2
    tspan = [0.0D0, 1.0D0]
    y0    = [1.0D0, 0.0D0]
    dt     = 0.1D0
    PP    = [1.0D0, 5.0D0, 0.0D0]

    n = CEILING((tspan(2) - tspan(1))/dt + 1.0D0)

    ALLOCATE(t_y(n, Neq+1))

    write (*, '( // "calling to solve ydot:")' )
    t_y = RK4( ydot, Neq, tspan, y0, dt, PP)

    write (*, '( // "calling to solve ydot2:")' )
    t_y = RK4( ydot2, Neq, tspan, y0, dt, PP)

END PROGRAM RK4_TEST
于 2013-07-04T04:59:27.630 回答