root/trunk/series-acceleration.lisp

;; Series acceleration.
;; Liam Healy, Wed Nov 21 2007 - 18:41
;; Time-stamp: <2008-03-09 19:21:47EDT series-acceleration.lisp>
;; $Id$

(in-package :gsl)

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

(cffi:defcstruct levin
  "The definition of Levin series acceleration for GSL."
  (size size)
  (position-in-array size)
  (terms-used size)
  (sum-plain :double)
  (q-num :pointer)
  (q-den :pointer)
  (dq-num :pointer)
  (dq-den :pointer)
  (dsum :pointer))

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

(defmfun allocate-levin (order)
  "gsl_sum_levin_u_alloc"
  ((order size))
  :c-return :pointer
  :export nil
  :index (letm levin)
  :documentation                        ; FDL
  "Allocate a workspace for a Levin u-transform of n
   terms.  The size of the workspace is O(2n^2 + 3n).")

(defmfun free-levin (levin)
  "gsl_sum_levin_u_free"
  ((levin :pointer))
  :c-return :void               ; Error in GSL manual, should be void?
  :export nil
  :index (letm levin)
  :documentation                        ; FDL
  "Free a previously allocated Levin acceleration.")

(defmfun accelerate (array levin)
  "gsl_sum_levin_u_accel"
  (((gsl-array array) :pointer) ((dim0 array) size) (levin :pointer)
   (accelerated-sum :double) (absolute-error :double))
  :documentation                        ; FDL
  "From the terms of a series in array, compute the extrapolated
   limit of the series using a Levin u-transform.  Additional
   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 
   w->sum_plain. The algorithm calculates the truncation error
   (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
;;;;****************************************************************************

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

(defmfun allocate-levin-truncated (order)
  "gsl_sum_levin_utrunc_alloc"
  ((order size))
  :c-return :pointer
  :export nil
  :index (letm levin-truncated)
  :documentation                        ; FDL
  "Allocate a workspace for a Levin u-transform of n
   terms, without error estimation.  The size of the workspace is
   O(3n).")

(defmfun free-levin-truncated (levin)
  "gsl_sum_levin_utrunc_free"
  ((levin :pointer))
  :export nil
  :index (letm levin-truncated)
  :documentation                        ; FDL
  "Free a previously allocated Levin acceleration.")

(defmfun accelerate-truncated (array levin)
  "gsl_sum_levin_utrunc_accel"
  (((gsl-array array) :pointer) ((dim0 array) size) (levin :pointer)
   (accelerated-sum :double) (absolute-error :double))
  :documentation                        ; FDL
  "From the terms of a series in array, compute the extrapolated
   limit of the series using a Levin u-transform.  Additional
   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 w->sum_plain. The algorithm
   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
   abserr_trunc.  To improve the reliability of the algorithm the
   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 ()
  (let ((maxterms 20)
        (sum 0.0d0)
        (zeta2 (/ (expt pi 2) 6)))
    (letm ((levin (levin maxterms))
           (array (vector-double-float maxterms)))
      (dotimes (n maxterms)
        (setf (maref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
        (incf sum (maref array n)))
      (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))))))

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