quadpack で数値積分

インストール

前準備として fpm をインストールしておく (ココを参照).

Netlibから FORTRAN77 のソースをダウンロードしてもいいが、modern fortranに移植された quadpack2 のが便利。

git clone https://github.com/jacobwilliams/quadpack.git
cd quadpack/
fpm build
fpm test
fpm install

必要に応じて .a と .mod を /usr/local/lib /usr/local/include などに移しておく。

計算

ガウス関数の積分 : $\int e^{-x^2} dx = \sqrt{\pi} \fallingdotseq 1.77245385091$ を計算してみた。

program test
use quadpack!, only: dqag
implicit none

real(8), parameter :: a=-1d4, b=1d4, epsabs=0.d0, epsrel=1d-6
integer, parameter :: key = 6, limit=1000, lenw = limit*4
real(8) :: abserr, result, work(lenw)
integer :: ier, iwork(limit), last, neval
real(8) :: pi = 4.d0*atan(1.d0)

call dqag(f, a, b, epsabs, epsrel, key, result, &
        abserr, neval, ier, limit, lenw, last, iwork, work)

if(ier /= 0) stop "error !"

print*, result, abs(result-sqrt(pi)), abserr

contains

real(8) function f(x)
  implicit none
  real(8), intent(in) :: x
  f = exp(-x**2)
end function f
end program test

• 変数 a,b は積分区間 -10^4 - +10^4
• epsabs epsrel は絶対誤差と相対誤差
• key = 6 は推奨設定 (keyを選ぶと積分方式を変えられる)
• limit と lenw は配列サイズの上限
• abserr に積分後の推定誤差、result に積分値が入る
• ier は計算が正常に終了したかの info で ier = 0 なら正常終了
• neval は f(x) が呼び出された回数

コンパイル・結果

gfortran test.f90 -I/usr/local/include/ -L/usr/local/lib -lquadpack
./a.out
1.7724538435571988        7.3483170659471853E-009   6.8051839196672293E-010

左から積分値、真値との差、推定誤差で、 -9乗のオーダーで正確に計算できていることがわかる。