GNU GSL Numerical Integration

7.7

top ← prev up next →

GNU GSL Numerical Integration

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

(require gsl-integration)

package: gsl-integration

Interface to GNU GSL Numerical Integration.

Library hides memory allocation and other low level C stuff.

GNU GSL has to be installed separately. Tested for version 2.5.

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-07))

(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-1

For more examples look at test.rkt source file.

Tip: When working with math formulas it is recommended to use a infix library, for better readability. For example: https://pkgs.racket-lang.org/package/infix

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-08

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
(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-08
  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.

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-08
  limit : exact-positive-integer? = 1000

QAGS adaptive integration with singularities

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-08
  limit : exact-positive-integer? = 1000

QAGP adaptive integration with known singular points

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-08
  limit : exact-positive-integer? = 1000

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

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-08
  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-08
  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-08
  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-08
  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-08
  n : exact-positive-integer? = 20

Romberg integration

3 Troubleshooting

Some linux systems have precompiled package for GNU GSL library. Howeever this package can be of an older version. It is recommended to compile and install GNU GSL library from source. Beware that installation directory from source can be different that from precompiled package.

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

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

Solution: Edit ~/.bashrc

LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/local/lib

export LD_LIBRARY_PATH

top ← prev up next →