numerical-integration.lisp 10.5 KB
 liam committed Jan 21, 2008 1 2 ;; Numerical integration ;; Liam Healy, Wed Jul 5 2006 - 23:14  Liam Healy committed Oct 25, 2008 3 ;; Time-stamp: <2008-10-25 11:35:53EDT numerical-integration.lisp>  lhealy committed Jul 24, 2008 4 ;; $Id$  liam committed Jul 06, 2006 5   liam committed Jul 12, 2006 6 7 ;;; To do: QAWS, QAWO, QAWF, more tests  liam committed Jul 06, 2006 8 9 (in-package :gsl)  liam committed Jul 07, 2006 10 11 12 13 ;;;;**************************************************************************** ;;;; QNG non-adaptive Gauss-Kronrod integration ;;;;****************************************************************************  liam committed Feb 18, 2008 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  liam committed Jul 06, 2006 17  ;; what these are if they are too large, it will do a minimum number  liam committed Jul 07, 2006 18  ;; of points anyway.  liam committed Jul 06, 2006 19  "gsl_integration_qng"  liam committed Feb 11, 2007 20  ((function :pointer)  liam committed Jul 06, 2006 21  (a :double) (b :double)  liam committed Feb 18, 2008 22  (absolute-error :double) (relative-error :double)  Liam Healy committed Aug 16, 2008 23  (result :double) (abserr :double) (neval sizet))  liam committed Feb 18, 2008 24  :documentation ; FDL  liam committed Jul 06, 2006 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 committed Feb 18, 2008 27 28  integral of f over (a,b) is achieved within the desired absolute and relative error limits, absolute-error and relative-error. The  liam committed Jul 07, 2006 29 30  function returns the final approximation, an estimate of the absolute error, and the number of function evaluations  liam committed Feb 18, 2008 31  used. The Gauss-Kronrod rules are designed in such a way  liam committed Jul 06, 2006 32 33 34  that each rule uses all the results of its predecessors, in order to minimize the total number of function evaluations.")  liam committed Jul 07, 2006 35 36 37 38 ;;;;**************************************************************************** ;;;; QAG adaptive Gauss-Kronrod integration ;;;;****************************************************************************  liam committed Feb 03, 2008 39 (defgo-s (integration-workspace size)  liam committed Feb 03, 2008 40  integration-workspace-alloc integration-workspace-free)  liam committed Jan 28, 2008 41   liam committed Feb 18, 2008 42 (defmfun integration-workspace-alloc (size)  Liam Healy committed Aug 16, 2008 43  "gsl_integration_workspace_alloc" ((size sizet))  liam committed Jul 07, 2006 44  :c-return :pointer  liam committed Jan 28, 2008 45  :export nil  liam committed Feb 03, 2008 46  :index (letm integration-workspace)  liam committed Feb 18, 2008 47 48  :documentation ; FDL "Allocate a workspace sufficient to hold n double  liam committed Jul 07, 2006 49 50  precision intervals, their integration results and error estimates.")  liam committed Feb 18, 2008 51 (defmfun integration-workspace-free (pointer)  Liam Healy committed Aug 23, 2008 52  "gsl_integration_workspace_free"  Liam Healy committed Sep 20, 2008 53  ((pointer :pointer))  liam committed Jul 07, 2006 54  :c-return :void  liam committed Jan 28, 2008 55  :export nil  liam committed Feb 03, 2008 56  :index (letm integration-workspace)  liam committed Feb 18, 2008 57 58  :documentation ; FDL "Free the memory associated with the workspace.")  liam committed Jul 07, 2006 59   liam committed Jul 08, 2006 60 61 62 63 (cffi:defcenum integrate-method :gauss15 :gauss21 :gauss31 :gauss41 :gauss51 :gauss61)  liam committed Feb 18, 2008 64 (defmfun integration-QAG  liam committed Jul 08, 2006 65  (function a b method limit workspace  liam committed Feb 18, 2008 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  liam committed Jul 08, 2006 68 69 70  ;; what these are if they are too large, it will do a minimum number ;; of points anyway. "gsl_integration_qag"  liam committed Feb 11, 2007 71  ((function :pointer)  liam committed Jul 08, 2006 72  (a :double) (b :double)  liam committed Feb 18, 2008 73  (absolute-error :double) (relative-error :double)  Liam Healy committed Aug 16, 2008 74  (limit sizet) (method integrate-method) (workspace :pointer)  liam committed Jul 08, 2006 75  (result :double) (abserr :double))  liam committed Feb 18, 2008 76  :documentation ; FDL  liam committed Jul 08, 2006 77  "Apply an integration rule adaptively until an estimate  liam committed Feb 18, 2008 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,  liam committed Jul 08, 2006 81  and an estimate of the absolute error. The integration rule  liam committed Feb 18, 2008 82  is determined by the value of method, which should  liam committed Jul 08, 2006 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 committed Feb 18, 2008 91 92  stored in the memory provided by workspace. The maximum number of subintervals is given by limit, which may not exceed the allocated  liam committed Jul 08, 2006 93 94  size of the workspace.")  liam committed Jul 09, 2006 95 96 97 98 ;;;;**************************************************************************** ;;;; QAGS adaptive integration with singularity ;;;;****************************************************************************  liam committed Feb 18, 2008 99 100 101 (defmfun integration-QAGS (function a b limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))  liam committed Jul 09, 2006 102  "gsl_integration_qags"  liam committed Feb 11, 2007 103  ((function :pointer)  liam committed Jul 09, 2006 104  (a :double) (b :double)  Liam Healy committed Aug 16, 2008 105  (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)  liam committed Jul 09, 2006 106  (result :double) (abserr :double))  liam committed Feb 18, 2008 107  :documentation ; FDL  liam committed Jul 09, 2006 108  "Apply the Gauss-Kronrod 21-point integration rule  liam committed Feb 18, 2008 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  liam committed Jul 09, 2006 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 committed Feb 18, 2008 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  liam committed Jul 09, 2006 119 120 121 122 123 124  workspace.") ;;;;**************************************************************************** ;;;; QAGP adaptive integration with known singular points ;;;;****************************************************************************  liam committed Feb 18, 2008 125 126 (defmfun integration-QAGP (function points limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))  liam committed Jul 09, 2006 127  "gsl_integration_qagp"  liam committed Feb 11, 2007 128  ((function :pointer)  Liam Healy committed Aug 23, 2008 129  ((mpointer points) :pointer) ((dim0 points) sizet)  Liam Healy committed Aug 16, 2008 130  (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)  liam committed Jul 09, 2006 131  (result :double) (abserr :double))  liam committed Feb 18, 2008 132  :documentation ; FDL  liam committed Jul 09, 2006 133 134  "Apply the adaptive integration algorithm QAGS taking account of the user-supplied locations of singular points. The array  liam committed Feb 18, 2008 135  pts of length npts should contain the endpoints of the  liam committed Jul 09, 2006 136 137  integration ranges defined by the integration region and locations of the singularities. For example, to integrate over the region  liam committed Feb 18, 2008 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  liam committed Jul 09, 2006 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 committed Feb 18, 2008 142  region then this routine will be faster than #'integration-QAGS.")  liam committed Jul 09, 2006 143 144 145 146 147  ;;;;**************************************************************************** ;;;; QAGI adaptive integration on infinite intervals ;;;;****************************************************************************  liam committed Feb 18, 2008 148 149 (defmfun integration-QAGi (function limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))  liam committed Jul 09, 2006 150  "gsl_integration_qagi"  liam committed Feb 11, 2007 151  ((function :pointer)  Liam Healy committed Aug 16, 2008 152  (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)  liam committed Jul 09, 2006 153  (result :double) (abserr :double))  liam committed Feb 18, 2008 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,  liam committed Jul 09, 2006 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 committed Feb 18, 2008 165 166 (defmfun integration-QAGiu (function a limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))  liam committed Jul 09, 2006 167  "gsl_integration_qagiu"  liam committed Feb 11, 2007 168  ((function :pointer) (a :double)  Liam Healy committed Aug 16, 2008 169  (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)  liam committed Jul 09, 2006 170  (result :double) (abserr :double))  liam committed Feb 18, 2008 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  liam committed Jul 09, 2006 176 177  and then integrated using the QAGS algorithm.")  liam committed Feb 18, 2008 178 179 (defmfun integration-QAGil (function b limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))  liam committed Jul 09, 2006 180  "gsl_integration_qagil"  liam committed Feb 11, 2007 181  ((function :pointer) (b :double)  Liam Healy committed Aug 16, 2008 182  (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)  liam committed Jul 09, 2006 183  (result :double) (abserr :double))  liam committed Feb 18, 2008 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  liam committed Jul 09, 2006 189 190 191 192 193 194  and then integrated using the QAGS algorithm.") ;;;;**************************************************************************** ;;;; QAWC adaptive integration for Cauchy principal values ;;;;****************************************************************************  liam committed Feb 18, 2008 195 196 197 (defmfun integration-QAWC (function a b c limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))  liam committed Jul 09, 2006 198  "gsl_integration_qawc"  liam committed Feb 11, 2007 199  ((function :pointer)  liam committed Jul 09, 2006 200  (a :double) (b :double) (c :double)  Liam Healy committed Aug 16, 2008 201  (absolute-error :double) (relative-error :double) (limit sizet) (workspace :pointer)  liam committed Jul 09, 2006 202  (result :double) (abserr :double))  liam committed Feb 18, 2008 203  :documentation ; FDL  liam committed Jul 09, 2006 204  "Compute the Cauchy principal value of the integral of  liam committed Feb 18, 2008 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}}  liam committed Jul 09, 2006 209  The adaptive bisection algorithm of QAG is used, with modifications to  liam committed Feb 18, 2008 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  liam committed Jul 09, 2006 212  it then a special 25-point modified Clenshaw-Curtis rule is used to control  liam committed Feb 18, 2008 213 214  the singularity. Further away from the singularity the algorithm uses an ordinary 15-point Gauss-Kronrod integration rule.")  liam committed Jul 09, 2006 215   liam committed Jul 07, 2006 216 217 218 219 220  ;;;;**************************************************************************** ;;;; Examples and unit test ;;;;****************************************************************************  liam committed Jan 21, 2008 221 (defun-single one-sine (x) (sin x))  liam committed Jul 07, 2006 222   liam committed Jan 21, 2008 223 #|  liam committed Dec 30, 2007 224 225 ;;; Parameters may be defined through the lexical environment: (let ((mult 2.0d0))  liam committed Jan 21, 2008 226  (defun-single two-sine (x) (sin (* mult x))))  liam committed Jan 21, 2008 227 |#  liam committed Jul 07, 2006 228   Liam Healy committed Oct 25, 2008 229 (save-test numerical-integration  liam committed Feb 18, 2008 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)))