Fortranで関数を引数にとる方法

Fortranで関数を引数にとる方法についての備忘録です．

Fortranで数値解析のコードを書いているとき，引数に関数をとりたいときってありますよね．
たとえば，関数$f(x)$のある点における微分$\frac{\partial f(x)}{\partial x}|_{x=a}$を計算したいときに，

dfdx = d_dx(func=f, x=a)

のように書けると嬉しい場合があります．
そこで今回は，$f(x) = x^2 + 1$を実装した以下の関数

function myfunc(x) result(ret)
    real(real64),parameter :: a=1.0d0
    real(real64),parameter :: b=1.0d0
    real(real64),intent(in) :: x
    real(real64) :: ret
    ret = a*x*x + b
end function

の数値微分を例に，実装例を紹介します．

呼び出し型のmainプログラム

この記事では，以下のmain.f90が動くように，DerivativeClassモジュールを作成します．

main.f90

program main
    use DerivativeClass
    implicit none
    
    ! 数値微分 
    ! 第一引数は関数，第二引数は微分点のxの値
    print *, d_dx(myfunc      ,x=1.0d0)
    
contains

! 微分したい関数
function myfunc(x) result(ret)
    real(real64),parameter :: a=1.0d0
    real(real64),parameter :: b=1.0d0
    real(real64),intent(in) :: x
    real(real64) :: ret
    ret = a*x*x + b
end function

end program

手順1．普通の関数と同様に実装を行う．

中心差分の式

$
\begin{equation}
\frac{\partial f(x)}{\partial x} = \frac{f(x+\epsilon) -f(x-\epsilon)}{2 \epsilon}
\end{equation}
$

により，数値微分を行う関数を書きます．ここで，funcは関数型の引数のつもりですが，今は未定義のまま，他を実装します．

DerivativeClass.f90

module DerivativeClass
    use iso_fortran_env
    implicit none
contains


! ########################################################
real(real64) function derivative_scalar(func,x,eps)
    ! >>> Define func()
    ! funcについてはまだ未定義
    ! <<<

    ! >>> arg
    real(real64),intent(in) :: x
    real(real64),optional,intent(in) :: eps !小さい値
    ! <<< 

    real(real64)  :: eps_val =dble(1.0e-4)
    if(present(eps) )then
        eps_val = eps
    endif
    
    ! >>> operation
    ! numerical derivative 
    ! 中心差分による数値微分
    derivative_scalar = (func(x+eps_val) - func(x-eps_val) )/(2.0d0*eps_val)
    ! <<<

end function
! ########################################################

end module DerivativeClass

手順2．Interface文により，第一引数を関数とする．

Interface文を使い，funcが，real(real64)型の戻り値・引数を持つことを定義します．
この定義により，real(real64)型の戻り値・引数をもつ任意の関数を，引数funcにとれるようになります．

DerivativeClass.f90

! ########################################################
real(real64) function derivative_scalar(func,x,eps)
    ! >>> Define func()
    interface 
        real(real64) function func(x)
            use iso_fortran_env
            real(real64),intent(in) :: x
        end function
    end interface 

    ! <<<

    ! >>> arg
    real(real64),intent(in) :: x
    real(real64),optional,intent(in) :: eps
    ! <<< 

    real(real64)  :: eps_val =dble(1.0e-4)
    if(present(eps) )then
        eps_val = eps
    endif
    
    ! >>> operation
    ! numerical derivative 
    
    derivative_scalar = (func(x+eps_val) - func(x-eps_val) )/(2.0d0*eps_val)
    ! <<<

end function
! ########################################################

end module DerivativeClass

あとは，コンパイルして動かすと，

   2.0000000000000000 

と数値微分が行われます．

おわりに

汎関数がつくれたり，最適化問題が実装できたりと素敵ですね！

参考

NAG様の記事(https://www.nag-j.co.jp/fortran/tips/tips_HowToPassArgumentsSafely.html)