数値計算・虎の穴#3 数値微分(2)～超実数を使った実装
この”数値計算・虎の穴”シリーズの記事を書いていたことをすっかり忘れていました。『数値微分(1)』を書いてから、ずいぶんと時間が経過してしまいました。もう少し短い間隔で記事を書くはずだったのですが･･･。今回はおちゃらけた部分がない超真面目な記事ですので、覚悟して下さい。
数値計算で使う微分は、簡単な数式なら元の式の微分から直接求めることができますが、複雑な式の場合は数値微分を使って計算します。数値微分といっても難しいことは無く、微分の定義式から極限の部分を省いたいわゆる差分を使います。
しかし、この差分で近似する方法も、うまく行く場合とそうでない場合があります。単調増加や単調減少などの大きな変化のない関数なら比較的うまく行きますが、増加と減少を繰り返すような振動を伴うような関数では、差分を行なう微小量の数値によっては、大きな誤差が生じます。
大きな誤差が生じない数値微分の方法には何種類かありますが、ここでは”超実数(Hyperreal numbers)”を使った実装法を紹介します。超実数というのは、通常の実数に微小部分を付加した拡張実数です。虚数計算の実数版を考えれば何となくわかると思います。ただし、虚数では虚数単位の二乗が-1となるのですが、超実数の微小単位の二乗は0となるのが大きな特徴です。
Fortranプログラムでは、超実数同志の計算をモジュール化し、微分演算を演算子”//”に多重定義しました。超実数モジュールでは、通常の四則演算に加え、標準の三角関数や対数関数、指数関数を実装しています。また、Fortran標準ではありませんが、ベッセル関数(J0とJ1)と誤差関数(erf)も実装しています。
実際のHyperrealモジュールはかなり長いのでここでは載せませんが、下記のサンプルプログラムを読めば、中身が想像できると思います。サンプルプログラムでは、実際の使用法を中心にテストをしています。test1からtest7までが1変数の微分で、test8とtest9が2変数の偏微分のテストプログラムになっています。テストプログラムでは、数値微分と実際の微分の比較を行なっていますが、数値微分は真値と一致します。
program main
call test1
call test2
call test3
call test4
call test5
call test6
call test7
call test8
call test9
end program main
!--------------------------------------------
! Test No.1
!
! f(x) = (c*x+d)/(a*x+b)
! df(x)/dx = (b*c-a*d)/(a*x+b)**2
!--------------------------------------------
subroutine test1
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.1 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)
x=independent(1,z)
y=(3.0*x+4.0)/(x+2.0)
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = 2.0/(xx+2.0)**2.0
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
enddo
end subroutine test1
!--------------------------------------------
! Test No.2
!
! f(x) = x+sqrt(x*x+1.0)
! df(x)/dx = (x+sqrt(x*x+1.0))/sqrt(x*x+1.0)
!--------------------------------------------
subroutine test2
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.2 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)
x=independent(1,z)
y=x+sqrt(x**2.0+1.0)
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = (xx+sqrt(xx**2.0+1.0))/sqrt(xx**2.0+1.0)
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
enddo
end subroutine test2
!--------------------------------------------
! Test No.3
!
! f(x) = (exp(x)*sin(x))/x
! df(x)/dx = exp(x)*(x*cos(x)+(x-1.0)*sin(x))/x**2
!--------------------------------------------
subroutine test3
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.3 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)
x=independent(1,z)
y=exp(x)*sin(x)/x
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = exp(xx)*(xx*cos(xx)+(xx-1.0)*sin(xx))/xx**2
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
enddo
end subroutine test3
!--------------------------------------------
! Test No.4
!
! f(x) = exp(sqrt(x))
! df(x)/dx = exp(sqrt(x))/(2.0*sqrt(x))
!--------------------------------------------
subroutine test4
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.4 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)
x=independent(1,z)
y=exp(sqrt(x))
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = exp(sqrt(xx))/(2.0*sqrt(xx))
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
enddo
end subroutine test4
!--------------------------------------------
! Test No.5
!
! f(x) = exp(x)*log(x)
! df(x)/dx = exp(x)*(1.0/x+log(x))
!--------------------------------------------
subroutine test5
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.5 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)
x=independent(1,z)
y=exp(x)*log(x)
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = exp(xx)*(1.0/xx+log(xx))
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
enddo
end subroutine test5
!--------------------------------------------
! Test No.6
!
! f(x) = asin(x)/sqrt(1.0-x**2)
! df(x)/dx = (sqrt(1.0-x**2)+x*asin(x))/(1.0-x**2)**(3/2)
!--------------------------------------------
subroutine test6
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.6 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)/20.0
x=independent(1,z)
y=asin(x)/sqrt(1.0-x**2.0)
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = (sqrt(1.0-xx**2.0)+xx*asin(xx))/(1.0-xx**2.0)**1.5
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
enddo
end subroutine test6
!--------------------------------------------
! Test No.7
!
! f(x) = exp(atan(x))
! df(x)/dx = exp(atan(x))/(1.0+x**2)
!--------------------------------------------
subroutine test7
use HyperReal
implicit none
type(Hreal) :: x, y
real(dpr) :: xx, z, diff1, diff2
integer :: i
write(6,'(/,a)') ' Test No.7 '
write(6,'(/,a)') ' No. f(x) df/dx.Cal df/dx.Theo '
do i=1,10
z=dble(i)
x=independent(1,z)
y=exp(atan(x))
! numerical solution
diff1 = y // x
! analytical solution
xx=real(x)
diff2 = exp(atan(xx))/(1.0+xx**2.0)
write(6,'(i4,3e24.16)') i,real(y),diff1,diff2
! write(6,'(e24.16)') (diff1-diff2)/
enddo
end subroutine test7
!--------------------------------------------
! Test No.8
!
! f(x,y) = (3*x*x+2*x+1)*(4*y**3+3*y*y+2*y+1)
! Df(x,y)/Dx = (6*x+2)*(4*y**3+3*y*y+2*y+1)
! Df(x,y)/Dy = (3*x*x+2*x+1)*(12*y**2+6*y+2)
!--------------------------------------------
subroutine test8
use HyperReal
implicit none
type(Hreal) :: x, y, z
real(dpr) :: xx, yy, w, pdiff1, pdiff2
integer :: i
write(6,'(/,a)') ' Test No.8 '
write(6,'(/,a)') ' No. f(x,y) Df/Dx.Cal Df/Dx.Theo '
do i=1,10
w=dble(i)/10.0
x=independent(1,w)
y=independent(2,w)
z=(3.0*x*x+2.0*x+1.0)*(4.0*y**3.0+3.0*y*y+2.0*y+1.0)
! numerical solution
pdiff1 = z // x ! Df(x,y)/Dx
! analytical solution
xx=real(x)
yy=real(y)
pdiff2 = (6.0*xx+2.0)*(4.0*yy**3.0+3.0*yy*yy+2.0*yy+1.0)
write(6,'(i4,3e24.16)') i,real(y),pdiff1,pdiff2
enddo
write(6,'(/,a)') ' No. f(x,y) Df/Dy.Cal Df/Dy.Theo '
do i=1,10
w=dble(i)/10.0
x=independent(1,w)
y=independent(2,w)
z=(3.0*x*x+2.0*x+1.0)*(4.0*y**3.0+3.0*y*y+2.0*y+1.0)
! numerical solution
pdiff1 = z // y ! Df(x,y)/Dy
! analytical solution
xx=real(x)
yy=real(y)
pdiff2 = (3.0*xx*xx+2.0*xx+1.0)*(12.0*yy**2.0+6.0*yy+2.0)
write(6,'(i4,3e24.16)') i,real(y),pdiff1,pdiff2
enddo
end subroutine test8
!--------------------------------------------
! Test No.9
!
! f(x,y) = exp(sin(x*y))
! Df(x,y)/Dx = exp(sin(x*y))*y*cos(x*y)
! Df(x,y)/Dy = exp(sin(x*y))*x*cos(x*y)
!--------------------------------------------
subroutine test9
use HyperReal
implicit none
type(Hreal) :: x, y, z
real(dpr) :: xx, yy, w, pdiff1, pdiff2
integer :: i
write(6,'(/,a)') ' Test No.9 '
write(6,'(/,a)') ' No. f(x,y) Df/Dx.Cal Df/Dx.Theo '
do i=1,10
w=dble(i)/10.0
x=independent(1,w)
y=independent(2,w)
z=exp(sin(x*y))
! numerical solution
pdiff1 = z // x ! Df(x,y)/Dx
! analytical solution
xx=st(x) ! standard part of HyperReal
yy=st(y) ! standard part of HyperReal
pdiff2 = exp(sin(xx*yy))*yy*cos(xx*yy)
write(6,'(i4,3e24.16)') i,real(y),pdiff1,pdiff2
enddo
write(6,'(/,a)') ' No. f(x,y) Df/Dy.Cal Df/Dy.Theo '
do i=11,20
w=dble(i)/10.0
x=independent(1,w)
y=independent(2,w)
z=exp(sin(x*y))
! numerical solution
pdiff1 = z // y ! Df(x,y)/Dy
! analytical solution
xx=real(x) ! function overload = st(x)
yy=real(y) ! function overload = st(y)
pdiff2 = exp(sin(xx*yy))*xx*cos(xx*yy)
write(6,'(i4,3e24.16)') i,real(y),pdiff1,pdiff2
enddo
end subroutine test9この記事が気に入ったらサポートをしてみませんか？