微分プログラム [fortran]

こんにちは　いつも夜に更新するのに、昼間のように書く管理人のさくらえびです。
今日は、基本中の基本。\(x\)を変数とする関数\(f(x)\)の微分プログラムを作成したいと思います。
functionの引数にfunctionを用いることで、任意の1変数関数\(f(x)\)の微分を求めます。
まずは、おさらい
微分はご存じのとおり、こんな定義式です。
\[
f'(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h}
\]
さて、これをコードに起こすのも実は簡単

function DIF_F(fun,x)
implicit none
double precision fun,x,DIF_F,dx
dx = 1.0E-8
dIF_F = (fun(x+dx)-fun(x))/dx
end function

以上！テストしてみましょう。
テストは、2次関数でやってみます。
使う式は\(f(x) = x^2 \)です。\(2x \)と\(f’(x) \)の結果を比較してみます。
テストコードは、こんな感じ

program test
implicit none
external func
double precision x,func,DIF_F
print *, "x","f'(x)","2x","error"
x = 0.d0
print *, x,DIF_F(func,x),2*x,DIF_F(func,x)-2*x
x = 1.d0
print *, x,DIF_F(func,x),2*x,DIF_F(func,x)-2*x
x = 2.d0
print *, x,DIF_F(func,x),2*x,DIF_F(func,x)-2*x
x = 3.d0
print *, x,DIF_F(func,x),2*x,DIF_F(func,x)-2*x
x = 4.d0
print *, x,DIF_F(func,x),2*x,DIF_F(func,x)-2*x
x = 5.d0
print *, x,DIF_F(func,x),2*x,DIF_F(func,x)-2*x
end program
function DIF_F(fun,x)
implicit none
double precision fun,x,DIF_F,dx
dx = 1.0E-8
dIF_F = (fun(x+dx)-fun(x))/dx
end function
function func(x)
implicit none
double precision func,x
func = x**2
end function


結果は下記

いい感じで結果が出てますね。
ちなみに、あまり\(h\)の値を小さくしすぎると、誤差が乗るので、お勧めではありません。
今回は倍精度実数を用いてますが、実数の時は、もうちょっと桁を落とした方がよさそうです。
by さくらえび

英語版はこちら：
Differential program in fortran