series-acceleration.lisp 4.37 KB
 liam committed Jan 28, 2008 1 2 ``````;; Series acceleration. ;; Liam Healy, Wed Nov 21 2007 - 18:41 `````` Liam Healy committed Dec 26, 2008 3 ``````;; Time-stamp: <2008-12-26 18:22:40EST series-acceleration.lisp> `````` lhealy committed Jul 24, 2008 4 ``````;; \$Id\$ `````` liam committed Nov 22, 2007 5 6 7 `````` (in-package :gsl) `````` Liam Healy committed Dec 26, 2008 8 9 ``````;;; /usr/include/gsl/gsl_sum.h `````` liam committed Nov 22, 2007 10 11 12 13 ``````;;;;**************************************************************************** ;;;; Creation and calculation of Levin series acceleration ;;;;**************************************************************************** `````` Liam Healy committed Dec 26, 2008 14 ``````(cffi:defcstruct levin-c `````` liam committed Nov 22, 2007 15 `````` "The definition of Levin series acceleration for GSL." `````` Liam Healy committed Aug 23, 2008 16 17 18 `````` (size sizet) (position-in-array sizet) (terms-used sizet) `````` liam committed Nov 22, 2007 19 20 21 22 23 24 25 `````` (sum-plain :double) (q-num :pointer) (q-den :pointer) (dq-num :pointer) (dq-den :pointer) (dsum :pointer)) `````` Liam Healy committed Dec 23, 2008 26 27 28 29 30 ``````(defmobject levin "gsl_sum_levin_u" ((order sizet)) "Levin u-transform" ; FDL "Make a workspace for a Levin u-transform of n terms. The size of the workspace is O(2n^2 + 3n).") `````` liam committed Nov 22, 2007 31 `````` `````` liam committed Feb 18, 2008 32 ``````(defmfun accelerate (array levin) `````` liam committed Nov 22, 2007 33 `````` "gsl_sum_levin_u_accel" `````` Liam Healy committed Dec 26, 2008 34 `````` (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer) `````` liam committed Nov 22, 2007 35 `````` (accelerated-sum :double) (absolute-error :double)) `````` liam committed Feb 04, 2008 36 37 38 `````` :documentation ; FDL "From the terms of a series in array, compute the extrapolated limit of the series using a Levin u-transform. Additional `````` liam committed Nov 22, 2007 39 40 41 `````` working space must be provided in levin. The extrapolated sum is returned with an estimate of the absolute error. The actual term-by-term sum is returned in `````` liam committed Feb 04, 2008 42 `````` w->sum_plain. The algorithm calculates the truncation error `````` liam committed Nov 22, 2007 43 44 45 46 47 48 49 50 `````` (the difference between two successive extrapolations) and round-off error (propagated from the individual terms) to choose an optimal number of terms for the extrapolation.") ;;;;**************************************************************************** ;;;; Acceleration with error estimation from truncation ;;;;**************************************************************************** `````` Liam Healy committed Dec 23, 2008 51 52 53 54 55 56 ``````(defmobject levin-truncated "gsl_sum_levin_utrunc" ((order sizet)) "truncated Levin u-transform" ; FDL "Make a workspace for a Levin u-transform of n terms, without error estimation. The size of the workspace is O(3n).") `````` liam committed Nov 22, 2007 57 `````` `````` liam committed Feb 18, 2008 58 ``````(defmfun accelerate-truncated (array levin) `````` liam committed Nov 22, 2007 59 `````` "gsl_sum_levin_utrunc_accel" `````` Liam Healy committed Dec 26, 2008 60 `````` (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer) `````` liam committed Nov 22, 2007 61 `````` (accelerated-sum :double) (absolute-error :double)) `````` liam committed Feb 04, 2008 62 63 64 `````` :documentation ; FDL "From the terms of a series in array, compute the extrapolated limit of the series using a Levin u-transform. Additional `````` liam committed Nov 22, 2007 65 66 `````` working space must be provided in levin. The extrapolated sum is returned with an estimate of the absolute error. The actual `````` liam committed Feb 04, 2008 67 `````` term-by-term sum is returned in w->sum_plain. The algorithm `````` liam committed Nov 22, 2007 68 69 70 `````` terminates when the difference between two successive extrapolations reaches a minimum or is sufficiently small. The difference between these two values is used as estimate of the error and is stored in `````` liam committed Feb 04, 2008 71 `````` abserr_trunc. To improve the reliability of the algorithm the `````` liam committed Nov 22, 2007 72 73 74 75 76 77 78 79 80 81 `````` extrapolated values are replaced by moving averages when calculating the truncation error, smoothing out any fluctuations.") ;;;;**************************************************************************** ;;;; Example ;;;;**************************************************************************** ;;; From Sec. 29.3 in the GSL manual. (defun acceleration-example () `````` liam committed Jan 28, 2008 82 83 84 `````` (let ((maxterms 20) (sum 0.0d0) (zeta2 (/ (expt pi 2) 6))) `````` Liam Healy committed Dec 26, 2008 85 86 `````` (let ((levin (make-levin maxterms)) (array (make-marray 'double-float :dimensions maxterms))) `````` liam committed Feb 04, 2008 87 `````` (dotimes (n maxterms) `````` liam committed Feb 24, 2008 88 89 `````` (setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float)) (incf sum (maref array n))) `````` liam committed Feb 04, 2008 90 91 92 93 `````` (multiple-value-bind (accelerated-sum error) (accelerate array levin) (format t "~&term-by-term sum =~f using ~d terms" sum maxterms) (format t "~&term-by-term sum =~f using ~d terms" `````` Liam Healy committed Dec 26, 2008 94 95 `````` (cffi:foreign-slot-value levin 'levin-c 'sum-plain) (cffi:foreign-slot-value levin 'levin-c 'terms-used)) `````` liam committed Feb 04, 2008 96 97 98 `````` (format t "~&exact value = ~f" zeta2) (format t "~&accelerated sum = ~f using ~d terms" accelerated-sum `````` Liam Healy committed Dec 26, 2008 99 `````` (cffi:foreign-slot-value levin 'levin-c 'terms-used)) `````` liam committed Feb 04, 2008 100 101 `````` (format t "~&estimated error = ~f" error) (format t "~&actual error = ~f" (- accelerated-sum zeta2)))))) `````` liam committed Jan 28, 2008 102 103 104 105 106 107 108 `````` ;; term-by-term sum =1.5961632439130233 using 20 terms ;; term-by-term sum =1.5759958390005426 using 13 terms ;; exact value = 1.6449340668482264 ;; accelerated sum = 1.6449340669228176 using 13 terms ;; estimated error = 0.00000000008883604962761638 ;; actual error = 0.00000000007459122208786084``````