GNU GSL Numerical Integration 🔗ℹ

by Petter Pripp <petter.pripp@yahoo.com>

Interface to GNU GSL Numerical Integration.

Library hides memory allocation and other low level C stuff.

GNU GSL has to be installed separately. Development version of gsl is preferred. In Ubuntu the package is: libgsl-dev.

Naming of functions and keywords follow GNU GSL documentation https://www.gnu.org/software/gsl/doc/html/integration.html

The source code is distributed under the GNU General Public License.

1 Example🔗ℹ

(require gsl-integration)

(define (f x)
  ( / (log x) (sqrt x)))

(define expected (- 4.0))

(define res  (qags f 0 1 #:epsrel 1e-7))

(if (= 0 (first res))
    (begin
      (displayln (string-append "result          = " (~a (second res))))
      (displayln (string-append "exact result    = " (~a expected)))
      (displayln (string-append "result          = " (~a (second res))))
      (displayln (string-append "estimated error = " (~a (third res))))
      (displayln (string-append "actual error    = " (~a (- (second res) expected)))))
    (error "status = " (~a (first res))))

; result          = -4.000000000000085
; exact result    = -4.0
; result          = -4.000000000000085
; estimated error = 1.354472090042691e-13
; actual error    = -8.526512829121202e-14

For more examples look at test.rkt source file.

2 Reference🔗ℹ

The functions will always return a list.

First element is status code. Success when code = 0, otherwise error.

Success list: 0, result. Thereafter one or both (see GNU GSL documentation): abserr , neveal.

Error list: codenr, gsl-symbol, message.

procedure
(qng f a b [#:epsabs epsabs #:epsrel epsrel])
→
(or/c (list/c integer? real? real? integer?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8

The QNG algorithm is a non-adaptive procedure which uses fixed Gauss-Kronrod-Patterson abscissae to sample the integrand at a maximum of 87 points. It is provided for fast integration of smooth functions.

When success, returns:

(0 result abserr neval)

When error, returns:

(codenr gsl-symbol message)

procedure
(qng-r f a b [#:epsabs epsabs #:epsrel epsrel])
→ (list/c real? real? integer?)
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8

Same as qng, but raises an exception when error.

When success, returns:

(result abserr neval)

procedure
(qag f
a
b
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit
#:key key])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000
  key : exact-positive-integer? = 2

The QAG algorithm is a simple adaptive integration procedure. The integration region is divided into subintervals, and on each iteration the subinterval with the largest estimated error is bisected. This reduces the overall error rapidly, as the subintervals become concentrated around local difficulties in the integrand.

When success, returns:

(0 result abserr)

When error, returns:

(codenr gsl-symbol message)

procedure
(qag-r f
a
b
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit
#:key key]) → (list/c real? real?)
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000
  key : exact-positive-integer? = 2

Same as qag, but raises an exception when error.

When success, returns:

(result abserr)

procedure
(qags f
a
b
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

QAGS adaptive integration with singularities

The presence of an integrable singularity in the integration region causes an adaptive routine to concentrate new subintervals around the singularity. As the subintervals decrease in size the successive approximations to the integral converge in a limiting fashion. This approach to the limit can be accelerated using an extrapolation procedure. The QAGS algorithm combines adaptive bisection with the Wynn epsilon-algorithm to speed up the integration of many types of integrable singularities.

procedure
(qags-r f
a
b
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit]) → (list/c real? real?)
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

Same as qags, but raises an exception when error.

When success, returns:

(result abserr)

procedure
(qagp f
pts
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  pts : (listof real?)
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

QAGP adaptive integration with known singular points

Dublicates points will removed and the points will be sorted.

procedure
(qagp-r f
pts
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit]) → (list/c real? real?)
  f : (-> flonum? flonum?)
  pts : (listof real?)
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

Same as qagp, but raises an exception when error.

When success, returns:

(result abserr)

procedure
(qagi f
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

QAGI adaptive integration on infinite interval (-\infty,+\infty)

procedure
(qagi-r f
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit]) → (list/c real? real?)
  f : (-> flonum? flonum?)
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

Same as qagi, but raises an exception when error.

When success, returns:

(result abserr)

procedure
(qagiu f
a
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

QAGIU adaptive integration on semi-infinite interval (a,+\infty)

procedure
(qagil f
b
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

QAGIL adaptive integration on semi-infinite interval (-\infty,b)

procedure
(qawc f
a
b
c
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  c : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

QAWC adaptive integration for Cauchy principal values

procedure
(cquad f
a
b
[ #:epsabs epsabs
#:epsrel epsrel
#:limit limit])
→
(or/c (list/c integer? real? real? integer?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  limit : exact-positive-integer? = 1000

CQUAD is a doubly-adaptive general-purpose quadrature routine which can handle most types of singularities, non-numerical function values such as Inf or NaN, as well as some divergent integrals. It generally requires more function evaluations than the integration routines in QUADPACK, yet fails less often for difficult integrands.

procedure
(romberg f
a
b
[ #:epsabs epsabs
#:epsrel epsrel
#:n n])
→
(or/c (list/c integer? real? integer?)
      (list/c integer? symbol? string?))
  f : (-> flonum? flonum?)
  a : real?
  b : real?
  epsabs : real? = 0
  epsrel : real? = 1e-8
  n : exact-positive-integer? = 20

Romberg integration

3 Troubleshooting🔗ℹ

3.1 Romberg🔗ℹ

If you get an error that it can not find Romberg, but not error on the other functions: You have and older version of GNU GSL on our system. Romberg was added at version 2.5

3.2 ffi-lib: couldn’t open ...🔗ℹ

If you get error: ffi-lib: couldn’t open "libgslcblas.so" (libgslcblas.so: cannot open shared object file: No such file or directory):

Solution 1: Install development version of gsl. In Ubuntu the package is: libgsl-dev.

Solution 2: Modify source code in wrap.rkt with version number. This apply if you have installed package with version number. Example: libgsl27

(define-ffi-definer gslcblas (ffi-lib "libgslcblas" #:global? #t))
(define-ffi-definer gsl (ffi-lib "libgsl" #:global? #t))

Modify this to

(define-ffi-definer gslcblas (ffi-lib "libgslcblas" '("0" #f) #:global? #t))
(define-ffi-definer gsl (ffi-lib "libgsl" '("27" #f) #:global? #t))

Solution 3: Specify absolute path in wrap.rkt without .so extension (package location may differ).

(define-ffi-definer gsl (ffi-lib "/usr/lib/x86_64-linux-gnu/libgslcblas" #:global? #t))
(define-ffi-definer gslcblas (ffi-lib "/usr/lib/x86_64-linux-gnu/libgsl" #:global? #t))

3.3 DrRacket freeze/crash🔗ℹ

C/C++ exception’s do not mix well with Racket exception’s. Try run at command line to see what exception is raised.