series-acceleration.lisp 5.31 KB
Newer Older
1 2
;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
3
;; Time-stamp: <2012-01-13 12:01:09EST series-acceleration.lisp>
4
;;
5
;; Copyright 2007, 2008, 2009, 2011 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22

(in-package :gsl)

23 24
;;; /usr/include/gsl/gsl_sum.h

25 26 27 28
;;;;****************************************************************************
;;;; Creation and calculation of Levin series acceleration
;;;;****************************************************************************

29
(defun levin-value (levin slot)
30
  (cffi:foreign-slot-value (mpointer levin) '(:struct levin-c) slot))
31

32
(defmobject levin "gsl_sum_levin_u"
33
  ((order :sizet))
34 35
  "Levin u-transform"
  :documentation			; FDL
36 37
  "Make a workspace for a Levin u-transform of n
   terms.  The size of the workspace is O(2n^2 + 3n).")
38

liam committed
39
(defmfun accelerate (array levin)
40
  "gsl_sum_levin_u_accel"
41
  (((grid:foreign-pointer array) :pointer) ((dim0 array) :sizet) ((mpointer levin) :pointer)
42
   (accelerated-sum (:pointer :double)) (abserr (:pointer :double)))
43
  :inputs (array)
44 45 46
  :documentation			; FDL
  "From the terms of a series in array, compute the extrapolated
   limit of the series using a Levin u-transform.  Additional
47 48 49
   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 
50
   w->sum_plain. The algorithm calculates the truncation error
51 52 53 54 55 56 57 58
   (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
;;;;****************************************************************************

59
(defmobject levin-truncated "gsl_sum_levin_utrunc"
60
  ((order :sizet))
61 62
  "truncated Levin u-transform"
  :documentation			; FDL
63 64 65
  "Make a workspace for a Levin u-transform of n
   terms, without error estimation.  The size of the workspace is
   O(3n).")
66

liam committed
67
(defmfun accelerate-truncated (array levin)
68
  "gsl_sum_levin_utrunc_accel"
69
  (((grid:foreign-pointer array) :pointer) ((dim0 array) :sizet) ((mpointer levin) :pointer)
70
   (accelerated-sum (:pointer :double)) (abserr (:pointer :double)))
71
  :inputs (array)
72 73 74
  :documentation			; FDL
  "From the terms of a series in array, compute the extrapolated
   limit of the series using a Levin u-transform.  Additional
75 76
   working space must be provided in levin.  The extrapolated sum is
   returned with an estimate of the absolute error.  The actual
77
   term-by-term sum is returned in w->sum_plain. The algorithm
78 79 80
   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
81
   abserr_trunc.  To improve the reliability of the algorithm the
82 83 84 85 86 87 88 89 90
   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.

91
(defun acceleration-example (&optional (print-explanation t))
92 93
  (let ((maxterms 20)
	(sum 0.0d0)
94
	(zeta2 (/ (expt dpi 2) 6)))
95
    (let ((levin (make-levin maxterms))
96
	  (array (grid:make-foreign-array 'double-float :dimensions maxterms)))
97
      (dotimes (n maxterms)
98 99
	(setf (grid:aref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
	(incf sum (grid:aref array n)))
100 101
      (multiple-value-bind (accelerated-sum error)
	  (accelerate array levin)
102 103 104 105 106 107 108 109 110 111 112 113 114 115
	(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)
116
		accelerated-sum
117 118
		error
		(- accelerated-sum zeta2))))))
119 120 121 122 123 124 125

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

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