series-acceleration.lisp 4.37 KB
Newer Older
1 2
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
3
;; Time-stamp: <2008-12-26 18:22:40EST 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 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).")
31

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

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

liam's avatar
liam committed
58
(defmfun accelerate-truncated (array levin)
59
  "gsl_sum_levin_utrunc_accel"
60
  (((c-pointer array) :pointer) ((dim0 array) sizet) ((mpointer levin) :pointer)
61
   (accelerated-sum :double) (absolute-error :double))
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
65 66
   working space must be provided in levin.  The extrapolated sum is
   returned with an estimate of the absolute error.  The actual
67
   term-by-term sum is returned in w->sum_plain. The algorithm
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
71
   abserr_trunc.  To improve the reliability of the algorithm the
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 ()
82 83 84
  (let ((maxterms 20)
	(sum 0.0d0)
	(zeta2 (/ (expt pi 2) 6)))
85 86
    (let ((levin (make-levin maxterms))
	  (array (make-marray 'double-float :dimensions maxterms)))
87
      (dotimes (n maxterms)
88 89
	(setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
	(incf sum (maref array n)))
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"
94 95
		(cffi:foreign-slot-value levin 'levin-c 'sum-plain)
		(cffi:foreign-slot-value levin 'levin-c 'terms-used))
96 97 98
	(format t "~&exact value     = ~f" zeta2)
	(format t "~&accelerated sum = ~f using ~d terms"
		accelerated-sum
99
		(cffi:foreign-slot-value levin 'levin-c 'terms-used))
100 101
	(format t "~&estimated error = ~f" error)
	(format t "~&actual error = ~f" (- accelerated-sum zeta2))))))
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