Commit 5f7b97c1 authored by liam's avatar liam

Added series acceleration; example works.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3249 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 37c9bc7e
......@@ -58,3 +58,7 @@ e.g. SBCL.
- reading and writing functions will require handling of C streams,
not provided.
- Would be nice to have common callback definition and usage.
- I already have a generic alloc and free; need to apply that to all
the later chapters. But this would require making gratuitous CLOS
objects?
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Sat Nov 17 2007 - 22:37
; modified: Wed Nov 21 2007 - 21:21
;********************************************************
;;; $Id: $
......@@ -171,4 +171,5 @@
;; numerical-integration to something in init, the
;; numerical-integration dependency won't be necessary.
(:file "numerical-differentiation" :depends-on (init numerical-integration))
(:file "chebyshev" :depends-on (init numerical-integration))))
(:file "chebyshev" :depends-on (init numerical-integration))
(:file "series-acceleration" :depends-on (init))))
......@@ -237,7 +237,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
</tr>
<tr>
<td><a href="Fast-Fourier-Transforms.html">FFTs</a></td>
<td>Not done</td>
<td>Not started (needs knowledge of C usage)</td>
</tr>
<tr>
<td><a href="Numerical-Integration.html">Numerical Integration</a></td>
......@@ -294,7 +294,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
</tr>
<tr>
<td><a href="Series-Acceleration.html">Series Acceleration</a></td>
<td>Not started</td>
<td>Done</td>
</tr>
<tr>
<td><a href="Wavelet-Transforms.html">Wavelet Transforms</a></td>
......
;********************************************************
; file: series-acceleration.lisp
; description: Series acceleration.
; date: Wed Nov 21 2007 - 18:41
; author: Liam Healy
; modified: Wed Nov 21 2007 - 21:21
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; Needs with- macro to allocate, initialize, and free object.
;;;;****************************************************************************
;;;; 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))
(defun-gsl allocate-levin (order)
"gsl_sum_levin_u_alloc"
((order :size))
:c-return :pointer
:documentation
"Allocate a workspace for a Levin @math{u}-transform of @var{n}
terms. The size of the workspace is @math{O(2n^2 + 3n)}.")
(defun-gsl free-levin (levin)
"gsl_sum_levin_u_free"
((levin :pointer))
:documentation
"Free a previously allocated Levin acceleration.")
(defun-gsl accelerate (array levin)
"gsl_sum_levin_u_accel"
(((gsl-array array) :pointer) ((dim0 array) :size) (levin :pointer)
(accelerated-sum :double) (absolute-error :double))
:documentation
"From the terms of a series in @var{array}, compute the extrapolated
limit of the series using a Levin @math{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
@code{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
;;;;****************************************************************************
(defun-gsl allocate-levin-truncated (order)
"gsl_sum_levin_utrunc_alloc"
((order :size))
:c-return :pointer
:documentation
"Allocate a workspace for a Levin @math{u}-transform of @var{n}
terms, without error estimation. The size of the workspace is
@math{O(3n)}.")
(defun-gsl free-levin-truncated (levin)
"gsl_sum_levin_utrunc_free"
((levin :pointer))
:documentation
"Free a previously allocated Levin acceleration.")
(defun-gsl 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
"From the terms of a series in @var{array}, compute the extrapolated
limit of the series using a Levin @math{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 @code{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
@var{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)
(levin (allocate-levin maxterms))
(sum 0.0d0)
(zeta2 (/ (expt pi 2) 6)))
(with-data (array vector-double maxterms)
(dotimes (n maxterms)
(setf (gsl-aref array n) (coerce (/ (expt (1+ n) 2)) 'double-float))
(incf sum (gsl-aref 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))))
(free-levin levin)))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment