;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2008-10-25 11:35:53EDT numerical-integration.lisp>
;; $Id$
;;; To do: QAWS, QAWO, QAWF, more tests
(in-package :gsl)
;;;;****************************************************************************
;;;; QNG non-adaptive Gauss-Kronrod integration
;;;;****************************************************************************
(defmfun integration-QNG
(function a b &optional (absolute-error 1.0d0) (relative-error 1.0d0))
;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
;; what these are if they are too large, it will do a minimum number
;; of points anyway.
"gsl_integration_qng"
((function :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double)
(result :double) (abserr :double) (neval sizet))
:documentation ; FDL
"Apply the Gauss-Kronrod 10-point, 21-point, 43-point and
87-point integration rules in succession until an estimate of the
integral of f over (a,b) is achieved within the desired
absolute and relative error limits, absolute-error and relative-error. The
function returns the final approximation, an estimate of
the absolute error, and the number of function evaluations
used. The Gauss-Kronrod rules are designed in such a way
that each rule uses all the results of its predecessors, in order to
minimize the total number of function evaluations.")
;;;;****************************************************************************
;;;; QAG adaptive Gauss-Kronrod integration
;;;;****************************************************************************
(defgo-s (integration-workspace size)
integration-workspace-alloc integration-workspace-free)
(defmfun integration-workspace-alloc (size)
"gsl_integration_workspace_alloc" ((size sizet))
:c-return :pointer
:export nil
:index (letm integration-workspace)
:documentation ; FDL
"Allocate a workspace sufficient to hold n double
precision intervals, their integration results and error estimates.")
(defmfun integration-workspace-free (pointer)
"gsl_integration_workspace_free"
((pointer :pointer))
:c-return :void
:export nil
:index (letm integration-workspace)
:documentation ; FDL
"Free the memory associated with the workspace.")
(cffi:defcenum integrate-method
:gauss15 :gauss21 :gauss31
:gauss41 :gauss51 :gauss61)
(defmfun integration-QAG
(function a b method limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
;; what these are if they are too large, it will do a minimum number
;; of points anyway.
"gsl_integration_qag"
((function :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double)
(limit sizet) (method integrate-method) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Apply an integration rule adaptively until an estimate
of the integral of f over (a,b) is achieved within the
desired absolute and relative error limits, absolute-error and
relative-error. The function returns the final approximation,
and an estimate of the absolute error. The integration rule
is determined by the value of method, which should
be chosen from the following symbolic names,
:gauss15 :gauss21 :gauss31 :gauss41 :gauss51 :gauss61
corresponding to the 15, 21, 31, 41, 51 and 61 point Gauss-Kronrod
rules. The higher-order rules give better accuracy for smooth functions,
while lower-order rules save time when the function contains local
difficulties, such as discontinuities.
On each iteration the adaptive integration strategy bisects the interval
with the largest error estimate. The subintervals and their results are
stored in the memory provided by workspace. The maximum number of
subintervals is given by limit, which may not exceed the allocated
size of the workspace.")
;;;;****************************************************************************
;;;; QAGS adaptive integration with singularity
;;;;****************************************************************************
(defmfun integration-QAGS
(function a b limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
"gsl_integration_qags"
((function :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Apply the Gauss-Kronrod 21-point integration rule
adaptively until an estimate of the integral of f over
(a,b) is achieved within the desired absolute and relative error
limits, absolute-error and relative-error. The results are extrapolated
using the epsilon-algorithm, which accelerates the convergence of the
integral in the presence of discontinuities and integrable
singularities. The function returns the final approximation from the
extrapolation, and an estimate of the absolute error. The subintervals
and their results are stored in the
memory provided by workspace. The maximum number of subintervals
is given by limit, which may not exceed the allocated size of the
workspace.")
;;;;****************************************************************************
;;;; QAGP adaptive integration with known singular points
;;;;****************************************************************************
(defmfun integration-QAGP
(function points limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
"gsl_integration_qagp"
((function :pointer)
((mpointer points) :pointer) ((dim0 points) sizet)
(absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Apply the adaptive integration algorithm QAGS taking
account of the user-supplied locations of singular points. The array
pts of length npts should contain the endpoints of the
integration ranges defined by the integration region and locations of
the singularities. For example, to integrate over the region
(a,b) with break-points at x_1, x_2, x_3 (where
a < x_1 < x_2 < x_3 < b) then an array with
(setf (data array) #(a x_1 x_2 x_3 b)) should be used.
If you know the locations of the singular points in the integration
region then this routine will be faster than #'integration-QAGS.")
;;;;****************************************************************************
;;;; QAGI adaptive integration on infinite intervals
;;;;****************************************************************************
(defmfun integration-QAGi
(function limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
"gsl_integration_qagi"
((function :pointer)
(absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the integral of the function f over the
infinite interval (-\infty,+\infty). The integral is mapped onto the
semi-open interval (0,1] using the transformation x = (1-t)/t,
\int_{-\infty}^{+\infty} dx \, f(x)
= \int_0^1 dt \, (f((1-t)/t) + f(-(1-t)/t))/t^2.
It is then integrated using the QAGS algorithm. The normal 21-point
Gauss-Kronrod rule of QAGS is replaced by a 15-point rule, because the
transformation can generate an integrable singularity at the origin. In
this case a lower-order rule is more efficient.")
(defmfun integration-QAGiu
(function a limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
"gsl_integration_qagiu"
((function :pointer) (a :double)
(absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the integral of the function f over the
semi-infinite interval (a,+\infty). The integral is mapped onto the
semi-open interval (0,1] using the transformation x = a + (1-t)/t,
int_{a}^{+\infty} dx, f(x) = \int_0^1 dt f(a + (1-t)/t)/t^2
and then integrated using the QAGS algorithm.")
(defmfun integration-QAGil
(function b limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
"gsl_integration_qagil"
((function :pointer) (b :double)
(absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the integral of the function f over the
semi-infinite interval (-\infty,b). The integral is mapped onto the
semi-open interval (0,1] using the transformation x = b - (1-t)/t,
\int_{-\infty}^{b} dx, f(x) = \int_0^1 dt, f(b - (1-t)/t)/t^2
and then integrated using the QAGS algorithm.")
;;;;****************************************************************************
;;;; QAWC adaptive integration for Cauchy principal values
;;;;****************************************************************************
(defmfun integration-QAWC
(function a b c limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
"gsl_integration_qawc"
((function :pointer)
(a :double) (b :double) (c :double)
(absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the Cauchy principal value of the integral of
f over (a,b), with a singularity at c,
I = \int_a^b dx, {f(x)/x - c} = lim_{epsilon -> 0}
{\int_a^{c-epsilon} dx, {f(x)/x - c} + \int_{c+epsilon}^b dx,
{f(x) \over x - c}}
The adaptive bisection algorithm of QAG is used, with modifications to
ensure that subdivisions do not occur at the singular point x = c.
When a subinterval contains the point x = c or is close to
it then a special 25-point modified Clenshaw-Curtis rule is used to control
the singularity. Further away from the singularity the algorithm
uses an ordinary 15-point Gauss-Kronrod integration rule.")
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(defun-single one-sine (x) (sin x))
#|
;;; Parameters may be defined through the lexical environment:
(let ((mult 2.0d0))
(defun-single two-sine (x) (sin (* mult x))))
|#
(save-test numerical-integration
(integration-qng one-sine 0.0d0 pi)
(letm ((ws (integration-workspace 20)))
(integration-QAG one-sine 0.0d0 pi :gauss15 20 ws))
(letm ((ws (integration-workspace 20)))
(integration-QAG one-sine 0.0d0 pi :gauss15 50 ws)))