複数の積分を直接評価することは、計算上困難です。おそらく重要度サンプリングを使用して、モンテカルロを使用する方が良いかもしれません。ただし、メソッドの検証には、力ずくの直接統合が重要な場合があります。
私が使用している統合ルーチンは、Luke Mo によって 1970 年頃に書かれた "QuadMo" です。これを再帰的にしてモジュールに入れました。QuadMo は、要求された統合精度を得るために必要なメッシュを改良しました。これは、QuadMo を使用して n 次元積分を行うプログラムです。
これは、G95 コンパイルを使用して、6 までの nDim のすべての次元で SD 0.1 で 0.5 を中心とするガウス分布を使用したプログラムの検証です。数秒で実行されます。
nDim ans expected nlvl
1 0.249 0.251 2
2 6.185E-02 6.283E-02 2 2
3 1.538E-02 1.575E-02 2 2 2
4 3.826E-03 3.948E-03 2 2 2 2
5 9.514E-04 9.896E-04 2 2 2 2 2
6 2.366E-04 2.481E-04 2 2 2 2 2 2
コードは次のとおりです。
!=======================================================================
module QuadMo_MOD
implicit none
integer::QuadMo_MinLvl=6,QuadMo_MaxLvl=24
integer,dimension(:),allocatable::QuadMo_nlvlk
real*8::QuadMo_Tol=1d-5
real*8,save,dimension(:),allocatable::thet
integer,save::nDim
abstract interface
function QuadMoFunct_interface(thet,k)
real*8::QuadMoFunct_interface
real*8,intent(in)::thet
integer,intent(in),optional::k
end function
end interface
abstract interface
function MultIntFunc_interface(thet)
real*8::MultIntFunc_interface
real*8,dimension(:),intent(in)::thet
end function
end interface
procedure(MultIntFunc_interface),pointer :: stored_func => null()
contains
!----------------------------------------------------------------------
recursive function quadMoMult(funct,lower,upper,k) result(ans)
! very powerful integration routine written by Luke Mo
! then at the Stanford Linear Accelerator Center circa 1970
! QuadMo_Eps is error tolerance
! QuadMo_MinLvl determines initial grid of 2**(MinLvl+1) + 1 points
! to avoid missing a narrow peak, this may need to be increased.
! QuadMo_Nlvl returns number of subinterval refinements required beyond
! QuadMo_MaxLvl
! Modified by making recursive and adding argument k
! for multiple integrals (GuthrieMiller@gmail.com)
real*8::ans
procedure(QuadMoFunct_interface) :: funct
real*8,intent(in)::lower,upper
integer,intent(in),optional::k
real*8::Middle,Left,Right,eps,est,fLeft,fMiddle,fRight
& ,fml,fmr,rombrg,coef,estl,estr,estint,area,abarea
real*8::valint(50,2), Middlex(50), Rightx(50), fmx(50), frx(50)
& ,fmrx(50), estrx(50), epsx(50)
integer retrn(50),i,level
level = 0
QuadMo_nlvlk(k) = 0
abarea = 0
Left = lower
Right = upper
if(present(k))then
fLeft = funct(Left,k)
fMiddle = funct((Left+Right)/2,k)
fRight = funct(Right,k)
else
fLeft = funct(Left)
fMiddle = funct((Left+Right)/2)
fRight = funct(Right)
endif
est = 0
eps = QuadMo_Tol
100 level = level+1
Middle = (Left+Right)/2
coef = Right-Left
if(coef.ne.0) go to 150
rombrg = est
go to 300
150 continue
if(present(k))then
fml = funct((Left+Middle)/2,k)
fmr = funct((Middle+Right)/2,k)
else
fml = funct((Left+Middle)/2)
fmr = funct((Middle+Right)/2)
endif
estl = (fLeft+4*fml+fMiddle)*coef
estr = (fMiddle+4*fmr+fRight)*coef
estint = estl+estr
area= abs(estl)+ abs(estr)
abarea=area+abarea- abs(est)
if(level.ne.QuadMo_MaxLvl) go to 200
QuadMo_nlvlk(k) = QuadMo_nlvlk(k)+1
rombrg = estint
go to 300
200 if(( abs(est-estint).gt.(eps*abarea)).or.
1(level.lt.QuadMo_MinLvl)) go to 400
rombrg = (16*estint-est)/15
300 level = level-1
i = retrn(level)
valint(level, i) = rombrg
go to (500, 600), i
400 retrn(level) = 1
Middlex(level) = Middle
Rightx(level) = Right
fmx(level) = fMiddle
fmrx(level) = fmr
frx(level) = fRight
estrx(level) = estr
epsx(level) = eps
eps = eps/1.4d0
Right = Middle
fRight = fMiddle
fMiddle = fml
est = estl
go to 100
500 retrn(level) = 2
Left = Middlex(level)
Right = Rightx(level)
fLeft = fmx(level)
fMiddle = fmrx(level)
fRight = frx(level)
est = estrx(level)
eps = epsx(level)
go to 100
600 rombrg = valint(level,1)+valint(level,2)
if(level.gt.1) go to 300
ans = rombrg /12
end function quadMoMult
!-----------------------------------------------------------------------
recursive function MultInt(k,func) result(ans)
! MultInt(nDim,func) returns multi-dimensional integral from 0 to 1
! in all dimensions of function func
! variable QuadMo_Mod: nDim needs to be set initially to number of dimensions
procedure(MultIntFunc_interface) :: func
real*8::ans
integer,intent(in)::k
stored_func => func
if(k.eq.nDim)then
if(allocated(thet))deallocate(thet)
allocate (thet(nDim))
if(allocated(QuadMo_nlvlk))deallocate(QuadMo_nlvlk)
allocate(QuadMo_nlvlk(nDim))
endif
if(k.eq.0)then
ans=func(thet)
return
else
ans=QuadMoMult(MultIntegrand,0d0,1d0,k)
endif
end function MultInt
!-----------------------------------------------------------------------
recursive function MultIntegrand(thetARG,k) result(ans)
real*8::ans
real*8,intent(in)::thetARG
integer,optional,intent(in)::k
if(present(k))then
thet(k)=thetARG
else
write(*,*)'MultIntegrand: not expected, k not present!'
stop
endif
ans=MultInt(k-1,stored_func)
end function MultIntegrand
!-----------------------------------------------------------------------
end module QuadMo_MOD
!=======================================================================
module test_MOD
use QuadMo_MOD
implicit none
contains
!-----------------------------------------------------------------------
real*8 function func(thet) ! multidimensional function
! this is the function defined in nDim dimensions
! in this case a Gaussian centered at 0.5 with SD 0.1
real*8,intent(in),dimension(:)::thet
func=exp(-sum(((thet-5d-1)/1d-1)
& *((thet-5d-1)/1d-1))/2)
end function func
!-----------------------------------------------------------------------
end module test_MOD
!=======================================================================
! test program to evaluate multiple integrals
use test_MOD
implicit none
real*8::ans
! these values are set for speed, not accuracy
QuadMo_MinLvl=2
QuadMo_MaxLvl=3
QuadMo_Tol=1d-1
write(*,*)' nDim ans expected nlvl'
do nDim=1,6
! expected answer is (0.1 sqrt(2pi))**nDim
ans=MultInt(nDim,func)
write(*,'(i10,2(1pg10.3),999(i3))')nDim,ans,(0.250663)**nDim
& ,QuadMo_nlvlk
enddo
end
!-----------------------------------------------------------------------
