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
モジュールを作成します.
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
は関数型の引数のつもりですが,今は未定義のまま,他を実装します.
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
にとれるようになります.
! ########################################################
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)