2013-07-04 41 views
0

我有兴趣编写一个函数,该函数将一个模块名称作为其输入之一来使用。例如,我写了一个用于解决ODE系统的Runge Kutta四阶积分器。我更愿意编写实际的Runge Kutta集成函数,以便为定义ODE系统的函数提供输入。在Fortran中传递模块名称作为函数输入

下面是一个包含限定常微分方程

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 

这里是龙格库塔积分模块

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

您可以将函数传递给求解器。这应该足以让你的求解器一般,能够处理各种功能。让你的求解器把一个函数当作一个参数,然后用你想要解决的特定函数的实际参数来调用它。或者你可以用一个函数指针来调用它,每次调用的指针都指向你想要解析的函数。在你的例子中,你没有传递类似DIFF_EQ的模块,你传递类似于YDOT的函数。这里是一个功能指针的例子:Function pointer arrays in 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 
+0

所以,你会建议我包含所有我的函数定义常微分方程的系统,如'DIFF_EQ'模块内'YDOT',然后通过常微分方程我想在那个时候解决的功能?例如,'DIFF_EQ'将包含'YDOT_1','YDOT_2',....,'YDOT_N'。这样,我的Runge Kutta求解器总是会使用DIFF_EQ,但是我传递给求解器的函数会有所不同。 – Stephen

+0

是的。一个或多个包含功能的模块。直接或通过函数指针传递你想要解决的任何函数到龙格库塔求解器。 Runge Kutta求解器可以使用模块,也可以有一个声明函数解决的接口声明。第二种方式更一般,如果以后你想解决不同模块中的功能。 –

+0

非常感谢代码示例。我已经阅读过有关接口块的内容,但并不完全了解它如何适合我的需求。你的代码示例应该为我清除很多东西。 – Stephen

相关问题