series-acceleration.lisp 4.67 KB
 liam committed Jan 28, 2008 1 2 ``````;; Series acceleration. ;; Liam Healy, Wed Nov 21 2007 - 18:41 `````` Liam Healy committed Jan 21, 2009 3 ``````;; Time-stamp: <2009-01-21 22:48:11EST series-acceleration.lisp> `````` lhealy committed Jul 24, 2008 4 ``````;; \$Id\$ `````` liam committed Nov 22, 2007 5 6 7 `````` (in-package :gsl) `````` 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 Jan 21, 2009 26 27 28 ``````(defun levin-value (levin slot) (cffi:foreign-slot-value (mpointer levin) 'levin-c slot)) `````` Liam Healy committed Dec 23, 2008 29 30 31 32 33 ``````(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 34 `````` `````` liam committed Feb 18, 2008 35 ``````(defmfun accelerate (array levin) `````` liam committed Nov 22, 2007 36 `````` "gsl_sum_levin_u_accel" `````` Liam Healy committed Dec 26, 2008 37 `````` (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer) `````` liam committed Nov 22, 2007 38 `````` (accelerated-sum :double) (absolute-error :double)) `````` liam committed Feb 04, 2008 39 40 41 `````` :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 42 43 44 `````` 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 45 `````` w->sum_plain. The algorithm calculates the truncation error `````` liam committed Nov 22, 2007 46 47 48 49 50 51 52 53 `````` (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 54 55 56 57 58 59 ``````(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 60 `````` `````` liam committed Feb 18, 2008 61 ``````(defmfun accelerate-truncated (array levin) `````` liam committed Nov 22, 2007 62 `````` "gsl_sum_levin_utrunc_accel" `````` Liam Healy committed Dec 26, 2008 63 `````` (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer) `````` liam committed Nov 22, 2007 64 `````` (accelerated-sum :double) (absolute-error :double)) `````` liam committed Feb 04, 2008 65 66 67 `````` :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 68 69 `````` 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 70 `````` term-by-term sum is returned in w->sum_plain. The algorithm `````` liam committed Nov 22, 2007 71 72 73 `````` 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 74 `````` abserr_trunc. To improve the reliability of the algorithm the `````` liam committed Nov 22, 2007 75 76 77 78 79 80 81 82 83 `````` 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. `````` Liam Healy committed Jan 21, 2009 84 ``````(defun acceleration-example (&optional (print-explanation t)) `````` liam committed Jan 28, 2008 85 86 87 `````` (let ((maxterms 20) (sum 0.0d0) (zeta2 (/ (expt pi 2) 6))) `````` Liam Healy committed Dec 26, 2008 88 89 `````` (let ((levin (make-levin maxterms)) (array (make-marray 'double-float :dimensions maxterms))) `````` liam committed Feb 04, 2008 90 `````` (dotimes (n maxterms) `````` liam committed Feb 24, 2008 91 92 `````` (setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float)) (incf sum (maref array n))) `````` liam committed Feb 04, 2008 93 94 `````` (multiple-value-bind (accelerated-sum error) (accelerate array levin) `````` Liam Healy committed Jan 21, 2009 95 96 97 98 99 100 101 102 103 104 105 106 107 108 `````` (when print-explanation (format t "term-by-term sum =~f using ~d terms~&" sum maxterms) (format t "term-by-term sum =~f using ~d terms~&" (levin-value levin 'sum-plain) (levin-value levin 'terms-used)) (format t "exact value = ~f~&" zeta2) (format t "accelerated sum = ~f using ~d terms~&" accelerated-sum (levin-value levin 'terms-used)) (format t "estimated error = ~f~&" error) (format t "actual error = ~f ~&" (- accelerated-sum zeta2))) (values sum maxterms (levin-value levin 'sum-plain) (levin-value levin 'terms-used) `````` liam committed Feb 04, 2008 109 `````` accelerated-sum `````` Liam Healy committed Jan 21, 2009 110 111 `````` error (- accelerated-sum zeta2)))))) `````` liam committed Jan 28, 2008 112 113 114 115 116 117 118 `````` ;; 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 `````` Liam Healy committed Jan 21, 2009 119 120 `````` (save-test series-acceleration (acceleration-example nil))``````