Commit 2217c408 by Liam Healy

### Added integration methods QAWS, QAWO, QAWF; new tests for numerical integration

Integration methods
integration-QAWS: adaptive integration for singular functions
integration-QAWO: adaptive integration for oscillatory functions
integration-QAWF: adaptive integration for Fourier integrals
have been added in a file numerical-integration-with-tables.lisp.
These require tables with are new mobjects.  All of the tests from GSL
in gsl-1.11/integration/test.c and gsl-1.11/integration/tests.c for
ported functions are now in either this file or
numerical-integration.lisp.
SBCL 64: TOTAL: 1522 assertions passed, 5 failed, 0 execution errors.
CCL 64:  TOTAL: 1525 assertions passed, 2 failed, 0 execution errors.
parent 45ca8e88
 ;; Numerical integration techniques that require tables ;; Liam Healy 2009-04-04 15:24:05EDT ;; Time-stamp: <2009-04-04 18:09:47EDT numerical-integration-with-tables.lisp> ;; $Id:$ (in-package :gsl) ;;; /usr/include/gsl/gsl_integration.h ;;;;**************************************************************************** ;;;; QAWS adaptive integration for singular functions ;;;;**************************************************************************** (defmobject qaws-table "gsl_integration_qaws_table" ((alpha :double) (beta :double) (mu :int) (nu :int)) "table for QAWS numerical integration method" :documentation ; FDL "Make and initialize a table for the QAWS adaptive integration method for singular functions. It a singular weight function W(x) with the parameters (alpha, beta, mu, nu), W(x) = (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x) where alpha > -1, beta > -1, and mu = 0, 1, nu = 0, 1. The weight function can take four different forms depending on the values of mu and nu, W(x) = (x-a)^alpha (b-x)^beta (mu = 0, nu = 0) W(x) = (x-a)^alpha (b-x)^beta log(x-a) (mu = 1, nu = 0) W(x) = (x-a)^alpha (b-x)^beta log(b-x) (mu = 0, nu = 1) W(x) = (x-a)^alpha (b-x)^beta log(x-a) log(b-x) (mu = 1, nu = 1) The singular points (a,b) do not have to be specified until the integral is computed, where they are the endpoints of the integration range." :initialize-suffix "set" ;; This defines the reinitialize-instance but causes the make-qaws-table ;; to initialize twice. :initialize-args ((alpha :double) (beta :double) (mu :int) (nu :int))) (defmfun integration-QAWS (function a b alpha beta mu nu &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (table (make-qaws-table alpha beta mu nu)) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qaws" ((callback :pointer) (a :double) (b :double) ((mpointer table) :pointer) (absolute-error :double) (relative-error :double) (limit sizet) ((mpointer workspace) :pointer) (result :double) (abserr :double)) :callbacks (callback gsl-function nil (function :double (:input :double) :slug)) :callback-dynamic (nil (function)) :documentation ; FDL "Compute the integral of the function f(x) over the interval (a,b) with the singular weight function (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x). The parameters of the weight function (alpha, beta, mu, nu) are used to make the default table. The integral is I = int_a^b dx f(x) (x-a)^alpha (b-x)^beta log^mu (x-a) log^nu (b-x). The adaptive bisection algorithm of QAG is used. When a subinterval contains one of the endpoints then a special 25-point modified Clenshaw-Curtis rule is used to control the singularities. For subintervals which do not include the endpoints an ordinary 15-point Gauss-Kronrod integration rule is used.") ;;;;**************************************************************************** ;;;; QAWO adaptive integration for oscillatory functions ;;;;**************************************************************************** (cffi:defcenum integrate-sine-cosine :cosine :sine) (defmobject qawo-table "gsl_integration_qawo_table" ((omega :double) (L :double) (trig integrate-sine-cosine) (n sizet)) "table for QAWO numerical integration method" :documentation ; FDL "Make a table describing a sine or cosine weight function W(x) with the parameters (omega, L), W(x) = sin(omega x) W(x) = cos(omega x) The parameter L must be the length of the interval over which the function will be integrated L = b - a. The choice of sine or cosine is made with the parameter trig which should be one of :cosine or :sine. This makes a table of the trigonometric coefficients required in the integration process. The parameter n determines the number of levels of coefficients that are computed. Each level corresponds to one bisection of the interval L, so that n levels are sufficient for subintervals down to the length L/2^n. An error of class 'table-limit-exceeded is signalled if the number of levels is insufficient for the requested accuracy." :initialize-suffix "set" ;; This defines the reinitialize-instance but causes the make-qawo-table ;; to initialize twice. :initialize-args ((omega :double) (L :double) (trig integrate-sine-cosine))) (defmfun integration-QAWO (function a omega L trig n &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (table (make-qawo-table omega L trig n)) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qawo" ((callback :pointer) (a :double) (absolute-error :double) (relative-error :double) (limit sizet) ((mpointer workspace) :pointer) ((mpointer table) :pointer) (result :double) (abserr :double)) :callbacks (callback gsl-function nil (function :double (:input :double) :slug)) :callback-dynamic (nil (function)) :documentation ; FDL "Use an adaptive algorithm to compute the integral of f over (a,b) with the weight function \sin(\omega x) or \cos(\omega x) defined by the table wf, I = \int_a^b dx f(x) sin(omega x) I = \int_a^b dx f(x) cos(omega x) The results are extrapolated using the epsilon-algorithm to accelerate the convergence of the integral. The function returns the final approximation from the extrapolation, result, and an estimate of the absolute error, abserr. 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. Those subintervals with large widths d where d\omega > 4 are computed using a 25-point Clenshaw-Curtis integration rule, which handles the oscillatory behavior. Subintervals with a small widths where d\omega < 4 are computed using a 15-point Gauss-Kronrod integration.") ;;;;**************************************************************************** ;;;; QAWF adaptive integration for Fourier integrals ;;;;**************************************************************************** (defmfun integration-QAWF (function a omega L trig n &optional (absolute-error *default-absolute-error*) (table (make-qawo-table omega L trig n)) (limit 1000) (workspace (make-integration-workspace limit)) (cycle-workspace (make-integration-workspace limit))) "gsl_integration_qawf" ((callback :pointer) (a :double) (absolute-error :double) (limit sizet) ((mpointer workspace) :pointer) ((mpointer cycle-workspace) :pointer) ((mpointer table) :pointer) (result :double) (abserr :double)) :callbacks (callback gsl-function nil (function :double (:input :double) :slug)) :callback-dynamic (nil (function)) :documentation ; FDL "This function attempts to compute a Fourier integral of the function f over the semi-infinite interval [a,+\infty). I = \int_a^{+\infty} dx f(x) sin(omega x) I = \int_a^{+\infty} dx f(x) cos(omega x) The parameter \omega and choice of \sin or \cos is taken from the table wf (the length L can take any value, since it is overridden by this function to a value appropriate for the fourier integration). The integral is computed using the QAWO algorithm over each of the subintervals, C_1 = [a, a + c] C_2 = [a + c, a + 2 c] ... = ... C_k = [a + (k-1) c, a + k c] where c = (2 floor(|\omega|) + 1) \pi/|\omega|. The width c is chosen to cover an odd number of periods so that the contributions from the intervals alternate in sign and are monotonically decreasing when f is positive and monotonically decreasing. The sum of this sequence of contributions is accelerated using the epsilon-algorithm. This function works to an overall absolute tolerance of abserr. The following strategy is used: on each interval C_k the algorithm tries to achieve the tolerance TOL_k = u_k abserr where u_k = (1 - p)p^{k-1} and p = 9/10. The sum of the geometric series of contributions from each interval gives an overall tolerance of abserr. If the integration of a subinterval leads to difficulties then the accuracy requirement for subsequent intervals is relaxed, TOL_k = u_k max(abserr, max_{i
 ;; Numerical integration ;; Liam Healy, Wed Jul 5 2006 - 23:14 ;; Time-stamp: <2009-03-31 22:09:43EDT numerical-integration.lisp> ;; Time-stamp: <2009-04-04 19:35:50EDT numerical-integration.lisp> ;; $Id$ ;;; To do: QAWS, QAWO, QAWF, more tests (in-package :gsl) ;;; /usr/include/gsl/gsl_integration.h ;;;;**************************************************************************** ;;;; Default error values ;;;;**************************************************************************** (export '(*default-absolute-error* *default-relative-error*)) (defparameter *default-absolute-error* 1.0d-5 "The default absolute error used in numerical integration.") (defparameter *default-relative-error* 0.0d0 "The default relative error used in numerical integration.") ;;;;**************************************************************************** ;;;; QNG non-adaptive Gauss-Kronrod integration ;;;;**************************************************************************** (defmfun integration-QNG (function a b &optional (absolute-error 1.0d0) (relative-error 1.0d0)) (function a b &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*)) ;; 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. ... ... @@ -26,7 +38,7 @@ :callbacks (callback gsl-function nil (function :double (:input :double) :slug)) :callback-dynamic (nil (function)) :documentation ; FDL :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 ... ... @@ -53,9 +65,11 @@ :gauss41 :gauss51 :gauss61) (defmfun integration-QAG (function a b method limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function a b method &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) ;; 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. ... ... @@ -92,9 +106,10 @@ ;;;;**************************************************************************** (defmfun integration-QAGS (function a b limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function a b &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qags" ((callback :pointer) (a :double) (b :double) ... ... @@ -122,12 +137,14 @@ ;;;;**************************************************************************** (defmfun integration-QAGP (function points limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function points &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qagp" ((callback :pointer) ((mpointer points) :pointer) ((dim0 points) sizet) ((c-pointer points) :pointer) ((dim0 points) sizet) (absolute-error :double) (relative-error :double) (limit sizet) ((mpointer workspace) :pointer) (result :double) (abserr :double)) :inputs (points) ... ... @@ -151,9 +168,11 @@ ;;;;**************************************************************************** (defmfun integration-QAGi (function limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qagi" ((callback :pointer) (absolute-error :double) (relative-error :double) (limit sizet) ... ... @@ -173,9 +192,11 @@ this case a lower-order rule is more efficient.") (defmfun integration-QAGiu (function a limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function a &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qagiu" ((callback :pointer) (a :double) (absolute-error :double) (relative-error :double) (limit sizet) ... ... @@ -191,9 +212,11 @@ and then integrated using the QAGS algorithm.") (defmfun integration-QAGil (function b limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function b &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qagil" ((callback :pointer) (b :double) (absolute-error :double) (relative-error :double) (limit sizet) ... ... @@ -213,9 +236,11 @@ ;;;;**************************************************************************** (defmfun integration-QAWC (function a b c limit &optional (absolute-error 1.0d0) (relative-error 1.0d0) (workspace (make-integration-workspace limit))) (function a b c &optional (absolute-error *default-absolute-error*) (relative-error *default-relative-error*) (limit 1000) (workspace (make-integration-workspace limit))) "gsl_integration_qawc" ((callback :pointer) (a :double) (b :double) (c :double) ... ... @@ -241,7 +266,97 @@ ;;;; Examples and unit test ;;;;**************************************************************************** ;;; CCL 1.2 returns a complex number when the exponent is a double ;;; float even if it's positive whole number value. This is a ;;; workaround. #+ccl (defun nn-expt (base power) (if (= power (floor power)) (expt base (floor power)) (expt base power))) (defun integration-test-f1 (alpha) (lambda (x) (* (expt x alpha) (log (/ x))))) (defun integration-test-f3 (alpha) (lambda (x) (cos (* (expt 2 alpha) (sin x))))) (defun integration-test-f11 (alpha) (lambda (x) (#+ccl nn-expt #-ccl expt (log (/ x)) (1- alpha)))) (defun integration-test-f15 (alpha) (lambda (x) (* x x (exp (* (- (expt 2 (- alpha))) x))))) (defun integration-test-f16 (alpha) (lambda (x) (cond ((and (= alpha 1.0d0) (zerop x)) 1.0d0) ((and (> alpha 1.0d0) (zerop x)) 0.0d0) (t (/ (#+ccl nn-expt #-ccl expt x (1- alpha)) (expt (1+ (* 10 x)) 2)))))) (defun integration-test-f454 (x) (* (expt x 3) (* (log (abs (* (- (expt x 2) 1.0d0) (- (expt x 2) 2.0d0))))))) (defun integration-test-f455 (x) (/ (log x) (1+ (* 100 x x)))) (defun integration-test-f459 (x) (/ (+ (* 5.0d0 (expt x 3)) 6.0d0))) (defun integration-test-myfn1 (x) (exp (- (- x) (expt x 2)))) (defun integration-test-myfn2 (alpha) (lambda (x) (exp (* alpha x)))) (save-test numerical-integration (integration-qng 'sin 0.0d0 pi) (integration-QAG 'sin 0.0d0 pi :gauss15 20) (integration-QAG 'sin 0.0d0 pi :gauss21 40)) (integration-qng 'sin 0.0d0 pi) (integration-QAG 'sin 0.0d0 pi :gauss15 20) (integration-QAG 'sin 0.0d0 pi :gauss21 40) ;; Tests from gsl-1.11/integration/test.c ;; Functions defined in gsl-1.11/integration/tests.c (integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.1d0 0.0d0) (integration-QNG (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.1d0 0.0d0) (integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.0d0 1.0d-9) (integration-QNG (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.0d0 1.0d-9) (integration-QNG (integration-test-f3 1.3d0) 0.3d0 2.71d0 0.0d0 1d-12) (integration-QNG (integration-test-f3 1.3d0) 2.71d0 0.3d0 0.0d0 1d-12) (integration-QNG (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.0d0 1.0d-13) (integration-QNG (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.0d0 1.0d-13) (integration-QNG (integration-test-f1 -0.9d0) 0.0d0 1.0d0 0.0d0 1.0d-3) ; error (integration-QNG (integration-test-f1 -0.9d0) 1.0d0 0.0d0 0.0d0 1.0d-3) ; error (integration-QAG (integration-test-f1 2.6d0) 0.0d0 1.0d0 :gauss15 0.0d0 1.0d-10 1000) (integration-QAG (integration-test-f1 2.6d0) 1.0d0 0.0d0 :gauss15 0.0d0 1.0d-10 1000) (integration-QAG (integration-test-f1 2.6d0) 0.0d0 1.0d0 :gauss21 1.0d-14 0.0d0 1000) (integration-QAG (integration-test-f1 2.6d0) 1.0d0 0.0d0 :gauss21 1.0d-14 0.0d0 1000) (integration-QAG ; roundoff error (integration-test-f3 1.3d0) 0.3d0 2.71d0 :gauss31 1.0d-14 0.0d0 1000) (integration-QAG ; roundoff error (integration-test-f3 1.3d0) 2.71d0 0.3d0 :gauss31 1.0d-14 0.0d0 1000) (integration-QAG ; singularity error (integration-test-f16 2.0d0) -1.0d0 1.0d0 :gauss51 1.0d-14 0.0d0 1000) (integration-QAG ; singularity error (integration-test-f16 2.0d0) 1.0d0 -1.0d0 :gauss51 1.0d-14 0.0d0 1000) (integration-QAG ; iteration limit error (integration-test-f16 2.0d0) -1.0d0 1.0d0 :gauss61 1.0d-14 0.0d0 3) (integration-QAG ; iteration limit error (integration-test-f16 2.0d0) 1.0d0 -1.0d0 :gauss61 1.0d-14 0.0d0 3) (integration-QAGS (integration-test-f1 2.6d0) 0.0d0 1.0d0 0.0d0 1d-10 1000) (integration-QAGS (integration-test-f1 2.6d0) 1.0d0 0.0d0 0.0d0 1d-10 1000) (integration-QAGS (integration-test-f11 2.0d0) 1.0d0 1000.0d0 1d-7 0.0d0 1000) (integration-QAGS (integration-test-f11 2.0d0) 1000.0d0 1.0d0 1d-7 0.0d0 1000) (integration-QAGiu 'integration-test-f455 0.0d0 0.0d0 1d-3 1000) (integration-QAGiu (integration-test-f15 5.0d0) 0.0d0 0.0d0 1d-7 1000) (integration-QAGiu (integration-test-f16 1.0d0) 99.9d0 1d-7 0.0d0 1000) (integration-QAGi 'integration-test-myfn1 1.0d-7 0.0d0 1000) (integration-QAGil (integration-test-myfn2 1.0d0) 1.0d0 1.0d-7 0.0d0 1000) (integration-QAGp 'integration-test-f454 (copy (vector 0.0d0 1.0d0 (sqrt 2.0d0) 3.0d0) 'vector-double-float) 0.0d0 1.0d-3 1000) (integration-QAWc 'integration-test-f459 -1.0d0 5.0d0 0.0d0 0.0d0 1.0d-3 1000))
 ... ... @@ -96,7 +96,7 @@ complex conjugate scalar product of two complex vectors of length 3: #2m(49.27d0 -13.49d0 32.5d0 42.73d0 -17.24d0 43.31d0)) #C(-2940.2118d0 1861.9380999999998d0)

There are over 1200 examples available from within GSLL with the

There are over 1500 examples available from within GSLL with the function examples. These examples also serve as a test suite for GSLL. ... ... @@ -334,6 +334,7 @@ histogram histogram2d histogram-pdf histogram2d-pdf basis-spline chebyshev hankel wavelet wavelet-workspace random-number-generator quasi-random-number-generator discrete-random polynomial-complex-workspace integration-workspace qaws-table qawo-table eigen-symm eigen-symmv eigen-herm eigen-hermv eigen-nonsymm eigen-nonsymmv eigen-gensymm eigen-gensymmv eigen-gen eigen-genv ... ... @@ -385,7 +386,7 @@ and arrays used internally or for function return. Time-stamp: <2009-04-04 09:58:19EDT index.html> Time-stamp: <2009-04-04 18:01:52EDT index.html> ... ...

 ;; Definition of GSLL system ;; Liam Healy ;; Time-stamp: <2009-04-04 09:55:40EDT gsll.asd> ;; Time-stamp: <2009-04-04 17:56:17EDT gsll.asd> ;; $Id$ (asdf:defsystem "gsll" ... ... @@ -178,6 +178,8 @@ :depends-on (init data random) :components ((:file "numerical-integration") (:file "numerical-integration-with-tables" :depends-on (numerical-integration)) (:file "monte-carlo") (:file "numerical-differentiation"))) (:module ordinary-differential-equations ... ...
 ;; Helpers that define a single GSL function interface ;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp ;; Time-stamp: <2009-04-01 21:36:13EDT defmfun-single.lisp> ;; Time-stamp: <2009-04-04 18:47:47EDT defmfun-single.lisp> ;; $Id:$ (in-package :gsl) ... ... @@ -99,12 +99,12 @@ (append before after (callback-symbol-set callback-dynamic callbacks (first callback-dynamic-variables)) (let ((auxstart (position '&aux arglist))) ;; &aux bindings are checked ;; &optional/&key/&aux defaults are checked (let ((auxstart (after-llk arglist))) (when auxstart (apply 'append (mapcar 'rest (subseq arglist (1+ auxstart)))))))))))) (mapcar 'rest (remove-if 'atom auxstart))))))))))) (first callback-dynamic-variables)) ,@(when documentation (list documentation)) ,(funcall body-maker name arglist gsl-name c-arguments key-args)) ... ...
 ;; Lisp forms ;; Liam Healy 2009-03-07 15:49:25EST forms.lisp ;; Time-stamp: <2009-04-02 22:34:03EDT forms.lisp> ;; Time-stamp: <2009-04-04 19:08:43EDT forms.lisp> (in-package :gsl) ... ... @@ -40,3 +40,20 @@ (when pos (nth pos cats))))) (defun after-llk (arglist) "The portion of the arglist from the first llk on." (loop for elt in arglist with seen = nil when (or seen (member elt *defmfun-llk*)) do (setf seen t) when seen collect elt)) #| ;;; Oddly, this gives a warning in SBCL but works the same. (defun after-llk (arglist) "The portion of the arglist from the first llk on." (when (intersection *defmfun-llk* arglist) (subseq arglist (some (lambda (itm) (position itm arglist)) *defmfun-llk*)))) |#
 ... ... @@ -2,17 +2,220 @@ (in-package :gsl) (LISP-UNIT:DEFINE-TEST NUMERICAL-INTEGRATION (LISP-UNIT::ASSERT-NUMERICAL-EQUAL