Numerical Integration Library

Help
1

Could you please direct me to the standard Fortran numerical integration library? I need to solve integrals using the Gauss-Hermite quadrature rule.

2

Welcome to the forum, @MichaelRedenti. Would the routines in GitHub - fortran-lang/stdlib: Fortran Standard Library help? I did not find this specific rule there, but perhaps others can point in the right direction.

3

Thank you for welcoming to the Fortran community @Arjen !! The standard library does have Gauss-Legendre integration rule (integration over a finite interval) but does not have Gauss-Hermite (integration over an infinite interval). Perhaps this is something that can be added.

1 Like
4

Here is the routine for Gauss-Hermite quadrature from Stroud & Secrest:

I haven’t used it in a while, so I’ve forgotten the accuracy it achieves.

Another library you could take a look at is IQPACK: https://dl.acm.org/doi/abs/10.1145/35078.214351

1 Like
5

Here’s a small program demonstrating a Gauss-Hermite rule with tabulated coefficients taken from a book by V. I. Krylov (1902-1994),

program gausshermite

implicit none

integer, parameter :: dp = kind(1.0d0)

! Gauss-Hermite quadrature formula for n = 20
!
!   taken from:
!       V.I. Krylov (1962). Approximate calculation of integrals
!       (translated by Arthur Stroud). The Macmillan Company, New York.
!   originally published in:
!       H.E. Salzer, R. Zucker, and R. Capuano (1952). Table of zeros
!       and wieght factors of the first twenty Hermite polynomals.
!       J. Res. Nat. Bur. Standards, Vol. 48, pp. 111-116
!
! (only the values for the positive half-line [0,\infty) are given)
!
real(dp), parameter :: x(*) = [ &
    0.2453407083009_dp, &
    0.7374737285454_dp, &
    1.2340762153953_dp, &
    1.7385377121166_dp, &
    2.2549740020893_dp, &
    2.7888060584281_dp, &
    3.3478545673832_dp, &
    3.9447640401156_dp, &
    4.6036824495507_dp, &
    5.3874808900112_dp]

real(dp), parameter :: w(*) = [ &
    0.4622436696006_dp, &
    0.2866755053628_dp, &
    0.1090172060200_dp, &
    0.2481052088746e-1_dp, &
    0.3243773342238e-2_dp, &
    0.2283386360163e-3_dp, &
    0.7802556478532e-5_dp, &
    0.1086069370769e-6_dp, &
    0.4399340992273e-9_dp, &
    0.2229393645534e-12_dp ]

real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp)
real(dp) :: approx

!
! \int_{-\infty}^{\infty} e^{-x^2} dx = \sqrt(\pi)
!

approx = sum([w, w])

print '(A)',    "Integral 1"
print '(A,G0)', "Approx (GH): ", approx
print '(A,G0)', "True value:  ", sqrt(pi)

!
! \int_{-\infty}^{\infty} e^{-x^2} x^2 dx = \sqrt(\pi)/2
!

approx = sum([w*(-x)**2, w*x**2])

print '(/,A)',  "Integral 2"
print '(A,G0)', "Approx (GH): ", approx
print '(A,G0)', "True value:  ", sqrt(pi)/2._dp

!
! \int_{-\infty}^{\infty} e^{-x^2} J_0(x) dx = \sqrt(\pi) e^{-1/8} I_0(1/8) 
!                                            \approx 1.5703011006678
!   where, J_0 is the Bessel function of order zero, and
!          I_0 is the modified Bessel function of order zero 
!

approx = sum([w*bessel_j0(-x), w*bessel_j0(x)])

print '(/,A)',  "Integral 3"
print '(A,G0)', "Approx (GH): ", approx
print '(A,G0)', "Approx:      ", 1.5703011006678_dp

!
! \int_{-\infty}^{\infty} \frac{1}{1 + 25 x^2} dx = \frac{\pi}{5}  
!
approx = sum([w*expxsqr_runge(-x), w*expxsqr_runge(x)])

print '(/,A)',  "Integral 4"
print '(A,G0)', "Approx (GH): ", approx
print '(A,G0)', "True value:  ", pi/5._dp

contains

   !> Runge's function, 1/(1 + 25*x^2) 
   !>   We multiply by exp(x^2) because of the weight in the Gauss-Hermite rule
   elemental real(dp) function expxsqr_runge(x)
      real(dp), intent(in) :: x
      expxsqr_runge = exp(x**2)*(1._dp/(1._dp + 25._dp*x**2))
   end function

end program

Output:

$ gfortran -Wall gausshermite.f90 
$ ./a.out
Integral 1
Approx (GH): 1.7724538509053740
True value:  1.7724538509055159

Integral 2
Approx (GH): 0.88622692545262205
True value:  0.88622692545275794

Integral 3
Approx (GH): 1.5703011006676562
Approx:      1.5703011006678000

Integral 4
Approx (GH): 0.52471814908103476
True value:  0.62831853071795862

As you can see the number of matching digits (apart from Runge’s function) is roughly equal to the precision with which the tabulated values are given. For Runge’s function, the Hermite polynomials obviously do a poor job. Here’s a look at the two function:

Runge's function

image

I remembered one more resource from Daniele Funaro, Professor of Numerical Analysis at UNIMORE:

The routines needed are called zehega (for zeros) and wehega (for weights).

1 Like
6

FYI: The modernized Quadpack has various methods for numerical integration (quadrature):

1 Like
7

For use on the infinite interval quadpack has:

QAGI : Handles integration over infinite intervals. The infinite range is mapped onto a finite interval and then the same strategy as in QAGS is applied.

I’m curious to test it on Runge’s function.