# Numerical Integration Library **URL:** \ **Category:** Help\ **Created:** [May 26, 2022, 9:34am UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557 "2022-05-26T09:34:45Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![MichaelRedenti](https://avatars.discourse-cdn.com/v4/letter/m/3ab097/32.png) [@MichaelRedenti](https://fortran-lang.discourse.group/u/MichaelRedenti)\ **Post date:** [May 26, 2022, 9:34am UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/1 "2022-05-26T09:34:45Z")
Could you please direct me to the standard Fortran numerical integration library? I need to solve integrals using the Gauss-Hermite quadrature rule. ---
**Author:** ![Arjen](https://avatars.discourse-cdn.com/v4/letter/a/b9bd4f/32.png) [@Arjen](https://fortran-lang.discourse.group/u/Arjen)\ **Post date:** [May 26, 2022, 9:43am UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/2 "2022-05-26T09:43:15Z")
Welcome to the forum, @MichaelRedenti. Would the routines in [GitHub - fortran-lang/stdlib: Fortran Standard Library](https://github.com/fortran-lang/stdlib) help? I did not find this specific rule there, but perhaps others can point in the right direction. ---
**Author:** ![MichaelRedenti](https://avatars.discourse-cdn.com/v4/letter/m/3ab097/32.png) [@MichaelRedenti](https://fortran-lang.discourse.group/u/MichaelRedenti)\ **Post date:** [May 26, 2022, 9:54am UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/3 "2022-05-26T09:54:35Z")
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. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [May 26, 2022, 10:32am UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/4 "2022-05-26T10:32:32Z")
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](https://dl.acm.org/doi/abs/10.1145/35078.214351) ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [May 26, 2022, 12:22pm UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/5 "2022-05-26T12:22:49Z")
Here’s a small program demonstrating a Gauss-Hermite rule with tabulated coefficients taken from a book by [V. I. Krylov (1902-1994)](http://www.mathnet.ru/eng/person35937), ```fortran 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: ```nohighlight $ 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](https://en.wikipedia.org/wiki/Runge%27s_phenomenon)) 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](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/9/9bf35d4730b8d4f3766dca338c08c46a07e91036.png) * * * I remembered one more resource from Daniele Funaro, Professor of Numerical Analysis at UNIMORE: - [Polynomial Approximation of Differential Equations | PDF, 6.3 MB](http://morespace.unimore.it/danielefunaro/wp-content/uploads/sites/2/2014/06/bube.pdf) - [FORTRAN Routines for Spectral methods](http://morespace.unimore.it/danielefunaro/routines/) The routines needed are called `zehega` (for zeros) and `wehega` (for weights). ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [May 26, 2022, 1:43pm UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/6 "2022-05-26T13:43:55Z")
FYI: The modernized Quadpack has various methods for numerical integration (quadrature): - [GitHub - jacobwilliams/quadpack: Modern Fortran QUADPACK Library for 1D numerical quadrature](https://github.com/jacobwilliams/quadpack) ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [May 26, 2022, 1:49pm UTC](https://fortran-lang.discourse.group/t/numerical-integration-library/3557/7 "2022-05-26T13:49:27Z")
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.