series-acceleration.lisp 5.02 KB
Newer Older
1 2
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2008-08-23 11:30:37EDT series-acceleration.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7 8 9 10 11 12 13

(in-package :gsl)

;;;;****************************************************************************
;;;; Creation and calculation of Levin series acceleration
;;;;****************************************************************************

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

24
(defgo-s (levin order) allocate-levin free-levin)
25

liam's avatar
liam committed
26
(defmfun allocate-levin (order)
27
  "gsl_sum_levin_u_alloc"
Liam Healy's avatar
Liam Healy committed
28
  ((order sizet))
29
  :c-return :pointer
30
  :export nil
31
  :index (letm levin)
32 33 34
  :documentation			; FDL
  "Allocate a workspace for a Levin u-transform of n
   terms.  The size of the workspace is O(2n^2 + 3n).")
35

liam's avatar
liam committed
36
(defmfun free-levin (levin)
37 38
  "gsl_sum_levin_u_free"
  ((levin :pointer))
39 40
  :c-return :void		; Error in GSL manual, should be void?
  :export nil
41
  :index (letm levin)
42
  :documentation			; FDL
43 44
  "Free a previously allocated Levin acceleration.")

liam's avatar
liam committed
45
(defmfun accelerate (array levin)
46
  "gsl_sum_levin_u_accel"
Liam Healy's avatar
Liam Healy committed
47
  (((c-pointer array) :pointer) ((dim0 array) sizet) (levin :pointer)
48
   (accelerated-sum :double) (absolute-error :double))
49 50 51
  :documentation			; FDL
  "From the terms of a series in array, compute the extrapolated
   limit of the series using a Levin u-transform.  Additional
52 53 54
   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 
55
   w->sum_plain. The algorithm calculates the truncation error
56 57 58 59 60 61 62 63
   (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
;;;;****************************************************************************

64
(defgo-s (levin-truncated order) allocate-levin-truncated free-levin-truncated)
65

liam's avatar
liam committed
66
(defmfun allocate-levin-truncated (order)
67
  "gsl_sum_levin_utrunc_alloc"
Liam Healy's avatar
Liam Healy committed
68
  ((order sizet))
69
  :c-return :pointer
70
  :export nil
71
  :index (letm levin-truncated)
72
  :documentation			; FDL
liam's avatar
liam committed
73
  "Allocate a workspace for a Levin u-transform of n
74
   terms, without error estimation.  The size of the workspace is
75
   O(3n).")
76

liam's avatar
liam committed
77
(defmfun free-levin-truncated (levin)
78 79
  "gsl_sum_levin_utrunc_free"
  ((levin :pointer))
80
  :export nil
81
  :index (letm levin-truncated)
82
  :documentation			; FDL
83 84
  "Free a previously allocated Levin acceleration.")

liam's avatar
liam committed
85
(defmfun accelerate-truncated (array levin)
86
  "gsl_sum_levin_utrunc_accel"
Liam Healy's avatar
Liam Healy committed
87
  (((c-pointer array) :pointer) ((dim0 array) sizet) (levin :pointer)
88
   (accelerated-sum :double) (absolute-error :double))
89 90 91
  :documentation			; FDL
  "From the terms of a series in array, compute the extrapolated
   limit of the series using a Levin u-transform.  Additional
92 93
   working space must be provided in levin.  The extrapolated sum is
   returned with an estimate of the absolute error.  The actual
94
   term-by-term sum is returned in w->sum_plain. The algorithm
95 96 97
   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
98
   abserr_trunc.  To improve the reliability of the algorithm the
99 100 101 102 103 104 105 106 107 108
   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 ()
109 110 111
  (let ((maxterms 20)
	(sum 0.0d0)
	(zeta2 (/ (expt pi 2) 6)))
112
    (letm ((levin (levin maxterms))
113
	   (array (vector-double-float maxterms)))
114
      (dotimes (n maxterms)
115 116
	(setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
	(incf sum (maref array n)))
117 118 119 120 121 122 123 124 125 126 127 128
      (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"
		(cffi:foreign-slot-value levin 'levin 'sum-plain)
		(cffi:foreign-slot-value levin 'levin 'terms-used))
	(format t "~&exact value     = ~f" zeta2)
	(format t "~&accelerated sum = ~f using ~d terms"
		accelerated-sum
		(cffi:foreign-slot-value levin 'levin 'terms-used))
	(format t "~&estimated error = ~f" error)
	(format t "~&actual error = ~f" (- accelerated-sum zeta2))))))
129 130 131 132 133 134 135

;; 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