series-acceleration.lisp 4.67 KB
Newer Older
1 2
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
3
;; Time-stamp: <2009-01-21 22:48:11EST series-acceleration.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_sum.h

10 11 12 13
;;;;****************************************************************************
;;;; Creation and calculation of Levin series acceleration
;;;;****************************************************************************

14
(cffi:defcstruct levin-c
15
  "The definition of Levin series acceleration for GSL."
Liam Healy's avatar
Liam Healy committed
16 17 18
  (size sizet)
  (position-in-array sizet)
  (terms-used sizet)
19 20 21 22 23 24 25
  (sum-plain :double)
  (q-num :pointer)
  (q-den :pointer)
  (dq-num :pointer)
  (dq-den :pointer)
  (dsum :pointer))

26 27 28
(defun levin-value (levin slot)
  (cffi:foreign-slot-value (mpointer levin) 'levin-c slot))

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).")
34

liam's avatar
liam committed
35
(defmfun accelerate (array levin)
36
  "gsl_sum_levin_u_accel"
37
  (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer)
38
   (accelerated-sum :double) (absolute-error :double))
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
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 
45
   w->sum_plain. The algorithm calculates the truncation error
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
;;;;****************************************************************************

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).")
60

liam's avatar
liam committed
61
(defmfun accelerate-truncated (array levin)
62
  "gsl_sum_levin_utrunc_accel"
63
  (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer)
64
   (accelerated-sum :double) (absolute-error :double))
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
68 69
   working space must be provided in levin.  The extrapolated sum is
   returned with an estimate of the absolute error.  The actual
70
   term-by-term sum is returned in w->sum_plain. The algorithm
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
74
   abserr_trunc.  To improve the reliability of the algorithm the
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.

84
(defun acceleration-example (&optional (print-explanation t))
85 86 87
  (let ((maxterms 20)
	(sum 0.0d0)
	(zeta2 (/ (expt pi 2) 6)))
88 89
    (let ((levin (make-levin maxterms))
	  (array (make-marray 'double-float :dimensions maxterms)))
90
      (dotimes (n maxterms)
91 92
	(setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
	(incf sum (maref array n)))
93 94
      (multiple-value-bind (accelerated-sum error)
	  (accelerate array levin)
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)
109
		accelerated-sum
110 111
		error
		(- accelerated-sum zeta2))))))
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
119 120

(save-test series-acceleration (acceleration-example nil))