interface - 实现一个 ODE 积分器，它可以将一般 f(x,t) 作为 fortran 中的输入

问题描述

我想解决类型的 ODE

在现代 fortran 中。我想把积分器（例如四阶龙格库塔）写得非常笼统，这样同一个积分器可以用于等式的不同右侧，至少包括这些不同的情况：

f(x, t) 是一个函数，其中 x 和返回值都是标量
f(x, t) 是一个函数，其中 x 和返回值都是数组（形状相同）
f(x, t) 是一个类型绑定过程（或类似过程），其中派生类型的目的是将一些基础数据插入到 (x, t) 给出的位置和时间

根据对类似问题的回答，我已经实现了下面包含的代码，它似乎按预期工作。

我的问题是我想保持rk4(x, t, h, f)一般性，这样xinf(x,t)可以是一个假定形状的数组（或者甚至是一个标量，理想情况下），但同时我希望能够指定它x，例如dimension(2)当我实际实现一个函数来插入一些 2D 数据。但是，如果我尝试修改evaluate下面示例中x的函数dimension(2)，那么当我尝试编译时，我会收到一个错误，抱怨接口不匹配。有没有办法解决这个问题？

module interpolator_module
    implicit none
    integer, parameter :: WP = kind(1.0D0)

    interface
        ! This is the general form of the right-hand side of an ODE
        function rhs(x, t) result( val )
            import :: WP
            real(WP), dimension(:),      intent(in) :: x
            real(WP),                    intent(in) :: t
            real(WP), dimension(size(x))            :: val
        end function
    end interface

    type interpolator_type
        ! This type would in practice store arrays,
        ! of discrete data to be interpolated.
        real(WP) :: stored_data
        procedure(rhs), nopass, pointer :: eval
    contains
        procedure :: init
    endtype

    class(interpolator_type), pointer :: interpolator

contains

    subroutine init( this, stored_data )
        implicit none
        class(interpolator_type), target :: this
        real(WP) :: stored_data
        this % stored_data = stored_data
        this % eval => evaluate
        interpolator => this
    end subroutine

    function evaluate(x, t) result( val )
        implicit none
        real(WP), dimension(:),      intent(in) :: x
        real(WP),                    intent(in) :: t
        real(WP), dimension(size(x))            :: val
        ! This is where interpolation would happen
        val = interpolator % stored_data * x
    end function
end module


program main
    use interpolator_module, only : interpolator_type
    implicit none
    integer, parameter :: WP = kind(1.0D0)
    type(interpolator_type) :: interp
    real(WP), dimension(2) :: x
    real(WP) :: t, h

    ! initialise interpolator with some data
    call interp % init(-0.1_WP)

    x = (/ 2.0_WP, 1.0_WP /)
    t = 0.0_WP
    h = 1.0_WP

    ! Example of calling rk1 with the "type-bound procedure"
    ! which evaluates an interpolator
    call rk4(x, t, h, interp % eval )
    print *, x

    ! Example of calling rk1 with analytical function
    call rk4(X, t, h, f )
    print *, x

    contains

    subroutine rk4(x, t, h, f)
        ! Makes one step with 4th-order Runge-Kutta.
        ! Calculates next position using timestep h.
        implicit none
        real(WP), intent(inout), dimension(:)   :: x
        real(WP), intent(inout)                 :: t
        real(WP), intent(in)                    :: h
        interface
            function f(x, t) result(val)
                import WP
                real(WP), dimension(:),      intent(in) :: x
                real(WP),                    intent(in) :: t
                real(WP), dimension(size(x))            :: val
            end function
        end interface
        ! Local variables
        real(WP), dimension(size(x))            :: k1, k2, k3, k4

        ! Evaluations of f(x, t)
        k1 = f(x,          t)
        k2 = f(x + k1*h/2, t + h/2)
        k3 = f(x + k2*h/2, t + h/2)
        k4 = f(x + k3,     t + h)

        ! Next position
        x  = x + h*(k1 + 2*k2 + 2*k3 + k4)/6
        t  = t + h
    end subroutine

    pure function f(x, t) result(val)
        implicit none
        real(WP), dimension(:),      intent(in) :: x
        real(WP),                    intent(in) :: t
        real(WP), dimension(size(x))            :: val
        val = -0.1_WP*x
    end function
end program

标签： interfacefortrannumericode

interface - 实现一个 ODE 积分器，它可以将一般 f(x,t) 作为 fortran 中的输入

问题描述

解决方案

推荐阅读