(require gsl-integration)

(define (fx)
(/ (log x)(sqrt x)))

(define expected(- 4.0))

(define res(qagsf01#:epsrel1e-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.

2Reference🔗 i

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
( qngfab[#:epsabsepsabs#:epsrelepsrel])

→
(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:

(0resultabserrneval)

When error, returns:

(codenrgsl-symbolmessage)

procedure
( qng-rfab[#:epsabsepsabs#:epsrelepsrel])
→(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:

(resultabserrneval)

procedure
( qag f
a
b
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit
#:keykey])

→
(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:

(0resultabserr)

When error, returns:

(codenrgsl-symbolmessage)

procedure
( qag-r f
a
b
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit
#:keykey]) → (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:

(resultabserr)

procedure
( qags f
a
b
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit]) → (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:

(resultabserr)

procedure
( qagp f
pts
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit]) → (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:

(resultabserr)

procedure
( qagi f
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit]) → (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:

(resultabserr)

procedure
( qagiu f
a
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:limitlimit])

→
(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
[ #:epsabsepsabs
#:epsrelepsrel
#:nn])

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

3Troubleshooting🔗 i

3.1Romberg🔗 i

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.2ffi-lib: couldn’t open ...🔗 i

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-definergslcblas(ffi-lib"libgslcblas"#:global?#t))
(define-ffi-definergsl(ffi-lib"libgsl"#:global?#t))

Modify this to

(define-ffi-definergslcblas(ffi-lib"libgslcblas"'("0"#f)#:global?#t))
(define-ffi-definergsl(ffi-lib"libgsl"'("27"#f)#:global?#t))

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

(define-ffi-definergsl(ffi-lib"/usr/lib/x86_64-linux-gnu/libgslcblas"#:global?#t))
(define-ffi-definergslcblas(ffi-lib"/usr/lib/x86_64-linux-gnu/libgsl"#:global?#t))

3.3DrRacket freeze/crash🔗 i

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

top ← prev up next →