Я хочу решить ODE типа
в современном фортране. Я хочу написать, что интегратор (например, Рунге-Кутта 4-го порядка) должен быть достаточно общим, чтобы один и тот же интегратор можно было использовать для разных правых частей уравнения, включая, по крайней мере, следующие разные случаи:
- f (x, t) - функция, где и x, и возвращаемое значение являются скалярами
- f (x, t) - функция, где x и возвращаемое значение являются массивами ( такой же формы)
- f (x, t) - это процедура с привязкой к типу (или аналогичная), где целью производного типа является интерполяция некоторых базовых данных в положение и время, заданные (x , t)
На основании ответа на аналогичный вопрос Я реализовал приведенный ниже код, который, кажется, работает как ожидалось.
Моя проблема что я хотел бы сохранить rk4(x, t, h, f) общий, такой, что x в f(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