numerical-integration.lisp 10.5 KB
Newer Older
liam's avatar
liam committed
1 2
;; Numerical integration
;; Liam Healy, Wed Jul  5 2006 - 23:14
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2008-10-25 11:35:53EDT numerical-integration.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5

6 7
;;; To do: QAWS, QAWO, QAWF, more tests

8 9
(in-package :gsl)

10 11 12 13
;;;;****************************************************************************
;;;; QNG non-adaptive Gauss-Kronrod integration
;;;;****************************************************************************

liam's avatar
liam committed
14 15 16
(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
17
  ;; what these are if they are too large, it will do a minimum number
18
  ;; of points anyway.
19
  "gsl_integration_qng"
20
  ((function :pointer)
21
   (a :double) (b :double)
liam's avatar
liam committed
22
   (absolute-error :double) (relative-error :double)
Liam Healy's avatar
Liam Healy committed
23
   (result :double) (abserr :double) (neval sizet))
liam's avatar
liam committed
24
  :documentation			; FDL
25 26
  "Apply the Gauss-Kronrod 10-point, 21-point, 43-point and
   87-point integration rules in succession until an estimate of the
liam's avatar
liam committed
27 28
   integral of f over (a,b) is achieved within the desired
   absolute and relative error limits, absolute-error and relative-error.  The
29 30
   function returns the final approximation, an estimate of
   the absolute error, and the number of function evaluations
liam's avatar
liam committed
31
   used.  The Gauss-Kronrod rules are designed in such a way
32 33 34
   that each rule uses all the results of its predecessors, in order to
   minimize the total number of function evaluations.")

35 36 37 38
;;;;****************************************************************************
;;;; QAG adaptive Gauss-Kronrod integration
;;;;****************************************************************************

39
(defgo-s (integration-workspace size)
40
	 integration-workspace-alloc integration-workspace-free)
41

liam's avatar
liam committed
42
(defmfun integration-workspace-alloc (size)
Liam Healy's avatar
Liam Healy committed
43
  "gsl_integration_workspace_alloc" ((size sizet))
44
  :c-return :pointer
45
  :export nil
46
  :index (letm integration-workspace)
liam's avatar
liam committed
47 48
  :documentation			; FDL
  "Allocate a workspace sufficient to hold n double
49 50
  precision intervals, their integration results and error estimates.")

liam's avatar
liam committed
51
(defmfun integration-workspace-free (pointer)
Liam Healy's avatar
Liam Healy committed
52
  "gsl_integration_workspace_free"
Liam Healy's avatar
Liam Healy committed
53
  ((pointer :pointer))
54
  :c-return :void
55
  :export nil
56
  :index (letm integration-workspace)
liam's avatar
liam committed
57 58
  :documentation			; FDL
  "Free the memory associated with the workspace.")
59

60 61 62 63
(cffi:defcenum integrate-method
  :gauss15 :gauss21 :gauss31
  :gauss41 :gauss51 :gauss61)

liam's avatar
liam committed
64
(defmfun integration-QAG
65
    (function a b method limit workspace
liam's avatar
liam committed
66 67
	      &optional (absolute-error 1.0d0) (relative-error 1.0d0))
  ;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
68 69 70
  ;; what these are if they are too large, it will do a minimum number
  ;; of points anyway.
  "gsl_integration_qag"
71
  ((function :pointer)
72
   (a :double) (b :double)
liam's avatar
liam committed
73
   (absolute-error :double) (relative-error :double)
Liam Healy's avatar
Liam Healy committed
74
   (limit sizet) (method integrate-method) (workspace :pointer)
75
   (result :double) (abserr :double))
liam's avatar
liam committed
76
  :documentation			; FDL
77
  "Apply an integration rule adaptively until an estimate
liam's avatar
liam committed
78 79 80
  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,
81
  and an estimate of the absolute error.  The integration rule
liam's avatar
liam committed
82
  is determined by the value of method, which should
83 84 85 86 87 88 89 90
  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
liam's avatar
liam committed
91 92
  stored in the memory provided by workspace.  The maximum number of
  subintervals is given by limit, which may not exceed the allocated
93 94
  size of the workspace.")

95 96 97 98
;;;;****************************************************************************
;;;; QAGS adaptive integration with singularity
;;;;****************************************************************************

liam's avatar
liam committed
99 100 101
(defmfun integration-QAGS
    (function a b limit workspace
	      &optional (absolute-error 1.0d0) (relative-error 1.0d0))
102
  "gsl_integration_qags"
103
  ((function :pointer)
104
   (a :double) (b :double)
Liam Healy's avatar
Liam Healy committed
105
   (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
106
   (result :double) (abserr :double))
liam's avatar
liam committed
107
  :documentation			; FDL
108
  "Apply the Gauss-Kronrod 21-point integration rule
liam's avatar
liam committed
109 110 111
   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
112 113 114
   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
liam's avatar
liam committed
115 116 117 118
   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
119 120 121 122 123 124
   workspace.")

;;;;****************************************************************************
;;;; QAGP adaptive integration with known singular points
;;;;****************************************************************************

liam's avatar
liam committed
125 126
(defmfun integration-QAGP
    (function points limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
127
  "gsl_integration_qagp"
128
  ((function :pointer)
Liam Healy's avatar
Liam Healy committed
129
   ((mpointer points) :pointer) ((dim0 points) sizet)
Liam Healy's avatar
Liam Healy committed
130
   (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
131
   (result :double) (abserr :double))
liam's avatar
liam committed
132
  :documentation			; FDL
133 134
  "Apply the adaptive integration algorithm QAGS taking
   account of the user-supplied locations of singular points.  The array
liam's avatar
liam committed
135
   pts of length npts should contain the endpoints of the
136 137
   integration ranges defined by the integration region and locations of
   the singularities.  For example, to integrate over the region
liam's avatar
liam committed
138 139
   (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
140 141
   (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
liam's avatar
liam committed
142
   region then this routine will be faster than #'integration-QAGS.")
143 144 145 146 147

;;;;****************************************************************************
;;;; QAGI adaptive integration on infinite intervals
;;;;****************************************************************************

liam's avatar
liam committed
148 149
(defmfun integration-QAGi
    (function limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
150
  "gsl_integration_qagi"
151
  ((function :pointer)
Liam Healy's avatar
Liam Healy committed
152
   (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
153
   (result :double) (abserr :double))
liam's avatar
liam committed
154 155 156 157
  :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,
158 159 160 161 162 163 164
   \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.")

liam's avatar
liam committed
165 166
(defmfun integration-QAGiu
    (function a limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
167
  "gsl_integration_qagiu"
168
  ((function :pointer) (a :double)
Liam Healy's avatar
Liam Healy committed
169
   (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
170
   (result :double) (abserr :double))
liam's avatar
liam committed
171 172 173 174 175
  :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
176 177
   and then integrated using the QAGS algorithm.")

liam's avatar
liam committed
178 179
(defmfun integration-QAGil
    (function b limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
180
  "gsl_integration_qagil"
181
  ((function :pointer) (b :double)
Liam Healy's avatar
Liam Healy committed
182
   (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
183
   (result :double) (abserr :double))
liam's avatar
liam committed
184 185 186 187 188
  :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
189 190 191 192 193 194
   and then integrated using the QAGS algorithm.")

;;;;****************************************************************************
;;;; QAWC adaptive integration for Cauchy principal values
;;;;****************************************************************************

liam's avatar
liam committed
195 196 197
(defmfun integration-QAWC
    (function a b c limit workspace
	      &optional (absolute-error 1.0d0) (relative-error 1.0d0))
198
  "gsl_integration_qawc"
199
  ((function :pointer)
200
   (a :double) (b :double) (c :double)
Liam Healy's avatar
Liam Healy committed
201
   (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)
202
   (result :double) (abserr :double))
liam's avatar
liam committed
203
  :documentation			; FDL
204
  "Compute the Cauchy principal value of the integral of
liam's avatar
liam committed
205 206 207 208
   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}}
209
   The adaptive bisection algorithm of QAG is used, with modifications to
liam's avatar
liam committed
210 211
   ensure that subdivisions do not occur at the singular point x = c.
   When a subinterval contains the point x = c or is close to
212
   it then a special 25-point modified Clenshaw-Curtis rule is used to control
liam's avatar
liam committed
213 214
   the singularity.  Further away from the singularity the algorithm
   uses an ordinary 15-point Gauss-Kronrod integration rule.")
215

216 217 218 219 220

;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

liam's avatar
liam committed
221
(defun-single one-sine (x) (sin x))
222

223
#|
224 225
;;; Parameters may be defined through the lexical environment:
(let ((mult 2.0d0))
liam's avatar
liam committed
226
  (defun-single two-sine (x) (sin (* mult x))))
227
|#
228

Liam Healy's avatar
Liam Healy committed
229
(save-test numerical-integration
liam's avatar
liam committed
230 231 232 233 234 235
  (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)))