nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

立体角積分 - Lebedev quadrature

数値計算で立体角積分をしたいとき、安直に二重積分するよりも、Lebedev求積法を使った方が効率が良い。
Lebedev quadrature - Wikipedia
http://people.maths.ox.ac.uk/beentjes/Essays/QuadratureSphere.pdf

効率が良いという意味は、メッシュ数が少なくても精度が出るという意味で、言い換えれば同じ精度を出すのに二重積分よりもずっと速く計算が出来るということである。

求積法は、つまるところ、積分を「上手く精度が出るようにメッシュ間隔等を調整しかつ各点に重み付けをした和」に変換したものである。
Lebedev法の場合には、

 \displaystyle
\int sin \theta \, d\theta \, \int d\phi \, f(\theta,\phi) \approx 4 \pi \sum^N_{i=1} \omega_i \, f(\theta_i, \phi_i )

というように 4 \pi は後で掛けることを想定している。

Lebedev求積法のFortranのフリーのソースコードとしては、以下のものが知られている。
http://ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F

subroutineの構成としては、

get_oh( ... )
LD0006( X, Y, Z, W, N )
LD0014( X, Y, Z, W, N )
...
LD5810( X, Y, Z, W, N )

という感じで、使うのはLDXXXX(...)の方である。get_ohは外から使うことは無い。

使い方のポイントを箇条書きすると、

  • 引数の"X, Y, Z, W, N"は全てoutputであるため、こちらから何か値を用意する必要は無い。
  • ただし、"X, Y, Z, W"は配列サイズを事前に指定する必要がある
  • 例えばLD0006(...)の"0006"はメッシュ数を表していて、これを使いたければ少なくとも6より大きい配列サイズの"X, Y, Z, W"を用意しておく。
  • "N"はメッシュ数であるためLD0006(...)ではN=6が返って来る。
  • "X Y Z"がLebedevメッシュであり、 ( \theta, \phi )に直したければ、"acos(z)"や"atan2(y,x)"で自分で変換する。
  • "W"が重み付けで、"sum( W )"はほぼ1になるようになっている。
  • どのLDXXXX(...)を使うべきかは、計算時間と計算精度との相談になる。

gfortranでコンパイルしてWの和が1になるかチェックする簡単なテストプログラム。

program check_Lebedev
implicit none
double precision, allocatable :: x(:), y(:), z(:), w(:)
! double precision, allocatable :: theta(:), phi(:)
integer :: n, n_lbdv

write( *, * ) 'number of mesh:'
read( *, * ) n
allocate( x(n), y(n), z(n), w(n) )
select case( n )
    case(    6 )
        call LD0006( x, y, z, w, n_lbdv )
    ...
    case( 5810 )
        call LD5810( x, y, z, w, n_lbdv )
    case default
        write( *, * ) ' n /= n_lbdv'
        stop
end select

write( *, * ) sum( w )

! allocate( theta( n_lbdv ), phi( n_lbdv ) )
! theta = acos( z )
! phi = atan2( y, x )

end program check_Lebedev

include "Lebedev-Laikov_CommentOutBy!.F"

リンク先のLebedev-Laikov.Fはコメントアウトが行頭の"c"で行われているため、sedか何かを使って"!"でコメントアウトするようにする必要がある(適当なコンパイルオプションを使えば、そんなことしなくても良いのかもしれないが)。

これで球面調和関数とかを掛けた計算とかが比較的簡単に出来るようになると思われる。