r"""This module provides some techniques for sequence acceleration.
Most of the theory and notation follows D. Laurie's discussion in
Appendix A. of [2]_. The general model is
that the sequence has the form:
.. math::
s_k = s + \sum_{j=1}^{m} c_{j} \phi_{k,j} + \eta_k
where :math:`s` is the converged result, :math:`\phi_{k,j}` is some set
of known residual dependencies with unknown coefficients :math:`c_j`
and :math:`\eta_k` is a residual that will be taken to be negligible
[1]_.
The acceleration will be summarized in a table of coefficients
:math:`s_{k,n}` as follows:
.. math::
\newcommand{\myarrows}{\;\mathllap{\rightarrow\kern-1.2em
\lower0.6em\hbox{\ensuremath{\nearrow}}}}
\mat{S} = \begin{pmatrix}
s_{1,0} & \myarrows s_{1,1} & \myarrows s_{1,2} & \myarrows s_{1,3}\\
s_{2,0} & \myarrows s_{2,1} & \myarrows s_{2,2}\\
s_{3,0} & \myarrows s_{3,1}\\
s_{4,0}
\end{pmatrix}
where the arrows show the dependencies: Each term :math:`s_{n,k}` is a
function of the terms :math:`\{s_{n}, s_{n+1}, \cdots s_{n+k}\}`. The
first column is thus the original sequence.
Linear Accelerations
====================
Here we consider linear sequence transformations :math:`X` such that
:math:`X(as_{k} + bt_{k}) = aX(s_{k}) + bX(t_{k})`. In particular, we
consider each column of :math:`\mat{S}` to be a linear transformation
of the previous column: :math:`s_{k,j} = X_{j}(s_{k,j-1})`. In this
case, we may provide an explicit construction using `n` terms that
eliminates the first `n` residuals by demanding:
.. math::
\begin{aligned}
s_{k+\delta} &= s_{k,n} + \sum_{j=1}^{n}c_{j}\phi_{k+\delta,j}, &
\delta &\in[0,1,\cdots,n].
\end{aligned}
:label: linear_model
This defines the entry :math:`s_{k,n}` in terms of the original
elements of the sequence below and to the left by defining :math:`n +
1` equations for the :math:`n+1` unknowns :math:`\{s_{k,n},
c_{1,\cdots,n}\}`. Note that for each :math:`s_{n,k}`, the set of
coefficients :math:`c_{j}` is *different* because this is not the
original sequence (for example, we have truncated and ignored the
error terms :math:`\eta_{k}`).
Naively, one may directly solve these systems of equations for each of
entries :math:`s_{k,n}`, but this takes :math:`\order(n^5)`
operations. There is a general procedure called the "E-algorithm"
that reduces this to :math:`\order(n^3)` which we now present. This
is implemented by :func:`e_algorithm`. For special forms of the
residuals :math:`\phi_{k,j}` specialized methods exist which are often
much faster. Some of these are described later.
.. note:: For linear models, one should obtain convergence either
proceeding down or to the right, but the key feature is that the
convergence of the rows to the right should be much faster: this
indicates that the underlying model is good. If the rows do not
converge much faster than the columns, one should consider a
different acceleration technique.
General E-Algorithm
-------------------
As we discussed above, each successive column of `S` may be expressed
in terms of this as an affine combination of the previous two entries
:math:`s_{k,j-1}` and :math:`s_{k+1,j-1}`:
.. math::
\begin{aligned}
s_{k,j} &= f_{k,j}s_{k,j-1} + (1-f_{k,j})s_{k+1, j-1}\\
&= s_{k+1, j-1} + \frac{r_{k,j}}{1-r_{k,j}}(s_{k+1,j-1}-s_{k,j-1})
\end{aligned}
(The affinity is required to preserve the linearity, for example, so
that :math:`s_{:,j} + x` and :math:`s_{:,j+1} + x` converge to the
same value :math:`s + x`.)
Start with the original sequence :math:`s_{k,0}` as the first column
of :math:`\mat{S}`. Consider this, along with the new sequence
:math:`s_{k,1} = X_{1}(s_{k,0})` and the sequence of residuals
:math:`\phi_{k,1}`. The transformation :math:`X` can be represented
as a matrix:
.. math::
\mat{X}_{1} = \begin{pmatrix}
f_{1,1} & 1 - f_{1,1}\\
& f_{2,1} & 1 - f_{2,1}\\
& & \ddots & \ddots
\end{pmatrix}
Equations :eq:`linear_model` for :math:`n=1` are then satisfied if
:math:`X_{1}(\phi_{k,1}) = 0` which defines the coefficients:
.. math::
f_{k,1} = \frac{\phi_{k,1}}{\phi_{k+1,1} - \phi_{k,1}}
= \frac{\phi_{k,1}}{\Delta\phi_{k,1}}
The E-algorithm is then to apply this transformation :math:`X_{1}` to
the rest of the columns of :math:`\phi_{k,j}` to obtain a new set of
residuals (the first column having been annihilated). Now treat
:math:`s_{k,1}` as the new original sequence and repeat the process.
This is implemented by :func:`e_algorithm`. In certain special cases,
one can directly derive formula for the coefficients :math:`r_{k,j}`.
The table :math:`\mat{S}` can them be directly constructed. This is
done in :func:`direct_algorithm`.
To estimate the error-sensitivity, we consider each entry
:math:`s_{k,j}(s_{k}, s_{k+1}, \cdots, s_j)` to be a function of the
original sequence entries. Given errors in the input :math:`s_j \pm
\delta_j`, an estimate of the error in the result is
.. math::
\delta s_{k,j} \approx \sqrt{
\sum_{n} \left(\delta_n\pdiff{s_{k,j}}{s_n}\right)^2}.
For linear algorithms, this is easy to compute since :math:`s_{k,j} =
\left[\mat{K}_{j}\cdot \vect{s}\right]_k` where :math:`\mat{K}_{j} =
\prod_{n=1}^{j}\mat{X}_{n}`. The squares of the errors :math:`(\delta
s_{k,j})^2` will be the diagonals of the matrix
:math:`\mat{K}\cdot\diag{\delta_n^2}\cdot\mat{K}^{T}`.
For some of the faster specialized algorithms, it is a bit cumbersome
to compute this exactly, so we directly propagate the errors as either
.. math::
\delta s_{k,j} \approx
\sqrt{(1-f_{k,j})^2\delta s_{k+1,j})^2
+ f_{k,j}^2(\delta s_{k+1,j})^2}
or
.. math::
\delta s_{k,j} \approx
\abs{1-f_{k,j}}\delta s_{k+1,j}
+ \abs{f_{k,j}}\delta s_{k+1,j})
This will overestimate the errors as cancellations that occur at
various stages are not taken into account, but in practise, this is
not very significant.
Richardson Extrapolation
------------------------
The Richardson extrapolation method considers models of the form
.. math::
\phi_{k,j} = h_k^{P_j}
:label: richardson
These most typically arise in problems such as differentiation,
integration, or discretization where there is a step-size :math:`h_k`
that is decreased at successive stages. The most general model
:eq:`richardson` requires the full E-algorithm, but there are several
special cases that admit analytic solutions:
1) Geometric progression of step sizes with constant `r`:
.. math::
\begin{aligned}
h_{k+1}/h_{k} &= r, & r_{k,j} &= r^{P_j}
\end{aligned}
2) Arithmetic progression of exponents:
.. math::
\begin{aligned}
P_{j} &= c j, &
r_{k,j} &= \left(\frac{h_{k+1}}{h_{k}}\right)^{c}
\end{aligned}
Modified Salzer's Extrapolation
-------------------------------
The modified Salzer's extrapolation is based on a model of the form
.. math:: \phi_{k,j} = \frac{\psi_k}{(k+k_0)^{j-1}}.
This provides another efficient method which can be derived by
dividing :eq:`linear_model` by :math:`\psi_{k+\delta}`:
.. math::
\begin{aligned}
\frac{s_{k+\delta}}{\psi_{k+\delta}} &=
\frac{s_{k,n}}{\psi_{k+\delta}} + \sum_{j=0}^{n-1}
\frac{c_{j+1}}{(k+\delta+k_{0})^j}, &
\delta \in [0,1,\cdots,n].
\end{aligned}
From this we can construct the solution in terms of the divided
difference operators :math:`D^{n}`: at :math:`\delta=0` which acts as
a difference of the sequence :math:`s_k` as a function of :math:`t_k =
(k + k_0)^{-1}`. The idea here is that the sum on the right is a
polynomial of order :math:`n-1` in `t_k`, and hence will be
annihilated by the `n` th difference operator at :math:`\delta=0`.
.. math::
D^0 s_k &= s_k\\
D^{n} s_k &= \frac{D^{n-1} s_{k+1} - D^{n-1} s_{k}}{t_{k+n} - t_k}\\
s_{k,n} &= \frac{D^{n}\left(\frac{s_k}{\psi_k}\right)}
{D^{n}\left(\frac{1}{\psi_k}\right)}
Semi-linear Accelerations
-------------------------
The previous linear acceleration techniques are fairly robust, and one
can generally prove that if the original sequence converges, then the
accelerated sequence will also converge to the same limit. However,
to obtain convergence, one must put good inputs into the model
:eq:`linear_model`.
Perhaps the most famous of these methods are those due to Levin.
These are generalizations of the Salzer expansion where the function
:math:`\psi_k` is expressed in terms of the residuals :math:`a_k =
s_{k+1} - s_k` of the sequence (as opposed to linear methods where
:math:`\psi_k` is determined a priori independently of the actual
sequence). Semi-linear methods often work well, even if no a priori
information about the sequence is know, however, they are also more
difficult to characterize. (For example, even if the transformation
is exact for two sequences, it may not be exact for them sum as the
transformations are not completely linear.)
To understand where the Levin transformation works, note that
:math:`a_k \sim`
The three most common Levin accelerations are:
*T*: :math:`\psi_k = a_k`
*U*: :math:`\psi_k = (k + k_0)a_k`
*W*: :math:`\psi_k = \frac{a_k^2}{a_{k+1} - a_{k}}`
Notes and References
====================
.. [1] Determining the form of :math:`\phi_{k,j}` is somewhat of an
art. It is not just sufficient to choose an accurate sequence: the
residual :math:`\eta_{k}` that remains after summing `m` terms must
really be negligible. For example, consider :math:`s =
\sum_{k=1}^{\infty} k^{\alpha}` with partial sums :math:`s_k`.
From the residual, one might be tempted to use :math:`\phi_{k,j} =
{k+j}^{\alpha}` which will ultimately be exact. However, after
summing only `n` terms, the residual will still be :math:`\eta_{k}
= \sum_{j=2m}^{\infty}{k+j}^{\alpha} \approx
\order[(2m)^{\alpha+1}]` which is certainly not negligible if the
series needs convergence acceleration!
Instead, one must find a good set of approximation that fall off
sufficiently rapidly which is why these types of sequences work
well with Salzar's modified algorithm.
.. [2] F. Bornemann, D. Laurie, S. Wagon, and J. Waldvogel, `The SIAM
100-Digit Challenge: A Study in High-Accuracy Numerical Computing
<http://www-m3.ma.tum.de/m3old/bornemann/challengebook/index.html>`_,
SIAM, 2004.
"""
from __future__ import division
__all__ = ['e_algorthm', 'direct_alorithm', 'salzar', 'levin', 'Aitkin']
import numpy as np
def e_algorithm(s, phi=None, r=None, ds=None):
r"""Return the triangular (`n` x `n`) acceleration table `a` for
`n` terms of the sequence `s[k]` with residuals `phi[k, j]` or
multipliers `r[k, j]`.
Parameters
----------
s : 1-d array of length `n`
`s[k]` should be the k'th term in the sequence.
phi : 2-d triangular array (`n` x `n`)
`phi[k, j]` should described the residuals according to the
model
.. math::
s_k = s + \sum_{j=1}^{n} c_{j} \phi_{k,j} + \eta_k.
Only the upper triangular portion where `k + j < n` will be
used. Only one of `phi` or `r` should be specified.
r : 2-d triangular array (`n` x `n`)
`r[k, j]` should be the multipliers defining the sequence
transformation through
.. math::
\begin{aligned}
s_{k,j}
= s_{k+1, j-1} + \frac{r_{k,j}}{1-r_{k,j}}(s_{k+1,j-1}-s_{k,j-1})
\end{aligned}
Only the upper triangular portion where `k + j < n` will be
used. Only one of `phi` or `r` should be specified.
ds : `None` or float or 1-d array of length `n`
This should be an estimate of the errors in the terms of the
sequence. If it is a single number, it is interpreted as the
relative error in the terms, otherwise it is taken to be the
array of absolute errors.
If it is not `None`, then the error estimate array `dS` will
also be returned.
Returns
-------
S : array (`n` x `n`)
Triangular acceleration array. The rows `S[k,:]` should
converge faster than the columns `S[:,n]`, otherwise the model
should not be trusted. Without round-off error, `S[0,-1]`
would be the best estimate, but at some point, round-off error
becomes significant.
dS : array (`n` x `n`)
If `ds` is provided and not `None`, then this array of error
estimates will also be returned.
Examples
--------
Consider :math:`s = \sum_{k=1}^{\infty} k^{\alpha}` with partial
sums :math:`s_k`:
.. math::
s_{k} = s - \sum_{j=1}^{m} {k+j}^{\alpha}\\
\phi_{k,j} = {k+j}^{\alpha}
>>> alpha = -3/2
>>> exact = 2.61237534868549
>>> N = 20
>>> k, j = map(np.double, np.ogrid[1:N+1, 1:N+1])
>>> s = np.cumsum(k**alpha)
>>> phi = (k+j)**alpha
>>> phi = (k)**(1+alpha)/k**(j-1)
>>> S = e_algorithm(s, phi)
Notes
-----
"""
s = np.asarray(s)
if r is None:
if phi is None:
raise ValueError("Must provide one of 'r' or 'phi'")
phi = np.asarray(phi)
else:
if phi is not None:
raise ValueError("Must provide only one of 'r' or 'phi'")
r = np.asarray(r)
N = len(s)
s.shape = (N, 1)
S = np.zeros((N, N), dtype=s.dtype)
S[:, 0] = s[:, 0]
if ds is not None:
ds = np.asarray(ds)
if ds.shape is ():
# Convert relative error into absolute errors if ds is a
# number
ds = abs(s)*ds
ds.shape = (N, 1)
dS = np.zeros((N, N), dtype=s.dtype)
dS[:, 0] = ds[:, 0]
K = np.eye(N, dtype=s.dtype) # Error matrix
for j in xrange(1, N):
if r is None:
f = phi[1:, 0]/np.diff(phi[:, 0])
else:
r_ = r[:-1, 0]
f = r_/(1 - r_)
X = np.diag(f, -1)[1:, :] + np.diag(1 - f, 1)[:-1, :]
if r is None:
assert np.all(abs(np.dot(X, phi[:, 0])) < 1e-14)
phi = np.dot(X, phi[:, 1:])
else:
r = r[:-1, 1:]
s = np.dot(X, s)
if ds is not None:
K = np.dot(X, K)
dS[:-j, j] = np.sqrt(np.diag(np.dot(K,ds**2*K.T)))
S[:-j, j] = s[:, 0]
if ds is None:
return S
else:
return S, dS
def direct_algorithm(s, r, ds=None):
r"""Return the triangular (`n` x `n`) acceleration table `a` for `n`
terms of the sequence `s[k]` with multipliers `r[k, j]`.
Provides the same answer as `e_algorithm` with `r` specified, but
a slightly more pessimistic error estimate (using the faster algorithm).
Parameters
----------
s : 1-d array of length `n`
`s[k]` should be the k'th term in the sequence.
r : 2-d triangular array (`n` x `n`)
`r[k, j]` should be the multipliers defining the sequence
transformation through
.. math::
\begin{aligned}
s_{k,j}
= s_{k+1, j-1} + \frac{r_{k,j}}{1-r_{k,j}}(s_{k+1,j-1}-s_{k,j-1})
\end{aligned}
Only the upper triangular portion where `k + j < n` will be
used.
ds : `None` or float or 1-d array of length `n`
This should be an estimate of the errors in the terms of the
sequence. If it is a single number, it is interpreted as the
relative error in the terms, otherwise it is taken to be the
array of absolute errors.
If it is not `None`, then the error estimate array `dS` will
also be returned.
Returns
-------
S : array (`n` x `n`)
Triangular acceleration array. The rows `S[k,:]` should
converge faster than the columns `S[:,n]`, otherwise the model
should not be trusted. Without round-off error, `S[0,-1]`
would be the best estimate, but at some point, round-off error
becomes significant.
dS : array (`n` x `n`)
If `ds` is provided and not `None`, then this array of error
estimates will also be returned.
Examples
--------
Consider numerically differentiating a function through the
sequence `s_h = [f(x+h) - f(x)]/h`:
.. math::
f'(x) = s_h + \sum_{j=1} \frac{h^{j}f^{(n)}(x)}{(1+j)!}
We can use the Richardson extrapolation here with the model
:math:`\phi_{h,j} = h^{j}`. If we decrease the step-size `h`
geometrically :math:`h_k = h_0 l^{-k}`, then this falls into case
1) with :math:`r = h_{k+1}/h_{k} = l^{-1}` and :math:`P_j = j` and
hence :math:`r_{k,j} = r^{P_j} = l^{-j}`.
>>> l = 2
>>> h0 = 0.1
>>> N = 7
>>> f = np.sin
>>> x = 1.0
>>> k, j = map(np.double, np.ogrid[1:N+1, 1:N+1])
>>> r = l**(-j) + 0*k
>>> h = h0*l**(-k)
>>> s = (f(x + h) - f(x))/h
>>> S = direct_algorithm(s, r)
>>> digits = np.floor(-np.log10(abs(S/np.cos(x) - 1))).astype(int)
>>> print digits
[[ 1 3 5 9 12 12 12]
[ 1 4 6 10 12 12 0]
[ 2 4 7 11 12 0 0]
[ 2 5 8 12 0 0 0]
[ 2 6 9 0 0 0 0]
[ 2 6 0 0 0 0 0]
[ 3 0 0 0 0 0 0]]
We see that the convergence across the top row is much better than
down the first column, indicating that the model is good, and we
achieve 13 places of accuracy, which is about all we can expect
because of round-off error (the decrease in accuracy in the first
row is a result of round-off error). Note that, even with an
optimal choice of `h ~ 1e-8`---balancing the round-off error with
the truncation error---direct application of this finite
difference method would only yield at most 8 correct digits.
Now consider a slightly problematic example of integrating a
function (Romberg integration)
.. math::
s = \int_0^1 f(x) \d{x} = -1
using the trapezoidal rule with :math:`h = 2^{-k}`
.. math::
s_h = h \sum_{n=1}^{2^k} f[h(n+\tfrac{1}{2})]
Typically, the error scales with even powers of `h`, so we would
typically use the Richardson extrapolation 1) with `r=0.5` and
`P_j = 2j`:
>>> N = 7
>>> f = np.sin
>>> exact = 1-np.cos(1.0)
>>> k, j = map(np.double, np.ogrid[1:N+1, 1:N+1])
>>> h = 2**(-k)
>>> s = [np.sum(h_*f((np.arange(1, 2**k_+1)-0.5)*h_))
... for (h_, k_) in zip(h, k)]
>>> Pj = 2*j
>>> r = 0.5**Pj + 0*k
>>> S = direct_algorithm(s, r)
>>> digits = np.floor(- np.log10(abs(S/exact - 1) + 1e-16)).astype(int)
>>> print digits
[[ 1 4 8 12 16 15 15]
[ 2 5 9 14 15 15 0]
[ 3 7 11 15 15 0 0]
[ 3 8 13 15 0 0 0]
[ 4 9 15 0 0 0 0]
[ 4 10 0 0 0 0 0]
[ 5 0 0 0 0 0 0]]
We see that this works quite well. Now consider `f(x) = ln(x)`
which has a singularity at `x=0`:
>>> f = np.log
>>> exact = -1
>>> s = [np.sum(h_*f((np.arange(1, 2**k_+1)-0.5)*h_))
... for (h_, k_) in zip(h, k)]
>>> S = direct_algorithm(s, r)
>>> digits = np.floor(-np.log10(abs(S/exact - 1))).astype(int)
>>> print digits
[[0 1 1 1 2 2 2]
[1 1 1 2 2 2 0]
[1 1 2 2 2 0 0]
[1 2 2 2 0 0 0]
[1 2 2 0 0 0 0]
[2 2 0 0 0 0 0]
[2 0 0 0 0 0 0]]
The crucial point here is that the rows are *not* converging
faster than the columns. This indicates that the models wrong.
The logarithmic singularity introduces an additional error of
linear in `h`, thus we must include this power in
:math:`P_j \in [1, 2, 4, 6, \cdots]`
>>> Pj = 2*(j-1)
>>> Pj[0, 0] = 1
>>> r = 0.5**Pj + 0*k
>>> S = direct_algorithm(s, r)
>>> digits = np.floor(-np.log10(abs(S/exact - 1))).astype(int)
>>> print digits
[[ 0 2 4 6 9 12 15]
[ 1 2 5 8 11 14 0]
[ 1 3 6 10 14 0 0]
[ 1 4 8 12 0 0 0]
[ 1 4 9 0 0 0 0]
[ 2 5 0 0 0 0 0]
[ 2 0 0 0 0 0 0]]
Notice now that the rows are converging much faster than the
columns... We now have a correct model.
Notes
-----
"""
s = np.asarray(s).ravel()
N = len(s)
S = np.zeros((N, N), dtype=s.dtype)
S[:,0] = s
if ds is not None:
ds = np.asarray(ds)
if ds.shape is ():
# Convert relative error into absolute errors if ds is a
# number
ds = abs(s)*ds
ds.shape = N
dS = np.zeros((N, N), dtype=s.dtype)
dS[:, 0] = ds
for j in xrange(1, N):
r_ = r[:-1,0]
f = r_/(r_ - 1)
s = s[1:] - f*np.diff(s)
if ds is not None:
ds = abs(1 - f)*ds[1:] + abs(f)*ds[:-1]
# Quadratic error estimate
# ds = np.sqrt(((1 - f)*ds[1:])**2 +
# (f*ds[:-1])**2)
dS[:-j, j] = ds
r = r[:-1,1:]
S[:-j, j] = s
if ds is None:
return S
else:
return S, dS
def salzer(s, psi, k0=0, ds=None, ds_psi=None, dpsi=None):
r"""Return the triangular (`n` x `n`) acceleration table `a` for `n`
terms of the sequence `s[k]` with Salzer factors `psi[k]`.
Parameters
----------
s : 1-d array of length `n`
`s[k]` should be the k'th term in the sequence.
psi : 1-d array of length `n`
`psi[k]` should be the multipliers in the model
.. math::
s_k = s + \sum_{j=1}^{n} c_{j}\frac{\psi_k}{(k + k_0)^{j-1}}
+ \eta_k
:label: Salzer
where `k` starts from 1. For example, when evaluating a sum
:math:`s = \sum_{k=1}^{\infty} f(x)`, then a good choice of
residual estimate is often the tail approximated as an integral
:math:`\psi_k \approx \int_{k}^{\infty} f(x)\d{x}`.
k0 : float
Constant in :eq:`Salzer`.
ds : `None` or float or 1-d array of length `n`
This should be an estimate of the errors in the terms of the
sequence. If it is a single number, it is interpreted as the
relative error in the terms, otherwise it is taken to be the
array of absolute errors.
If it is not `None`, then the error estimate array `dS` will
also be returned.
ds_psi, dpsi : `None` or float or 1-d array of length `n`
These arrays should contain the absolute errors in `s/psi` and
`psi` in cases where `psi` is expressed in terms of the
sequence `s` (such as the Levin transformations). Note, these
are only used if `ds` is also provided.
Returns
-------
S : array (`n` x `n`)
Triangular acceleration array. The rows `S[k,:]` should
converge faster than the columns `S[:,n]`, otherwise the model
should not be trusted. Without round-off error, `S[0,-1]`
would be the best estimate, but at some point, round-off error
becomes significant.
dS : array (`n` x `n`)
If `ds` is provided and not `None`, then this array of error
estimates will also be returned.
Examples
--------
Consider :math:`s = \sum_{k=1}^{\infty} k^{-\alpha-1}` with partial
sums :math:`s_k`:
.. math::
s_{k} = s - \sum_{j=k}^{\infty} k^{-\alpha-1}.
For large `k`, the tail can be well approximated by the integral
:math:`\psi_k \approx \int_k^\infty k^{-\alpha-1} =
k^{-\alpha}/\alpha`.
>>> alpha = 1/2
>>> exact = 2.61237534868549
>>> N = 12
>>> k = np.arange(1, N+1)
>>> s = np.cumsum(k**(-alpha-1))
>>> psi = k**(-alpha)
>>> S = salzer(s, psi, k0=0)
>>> digits = np.floor(-np.log10(abs(S/exact - 1))).astype(int)
>>> print digits
[[ 0 0 1 2 4 4 5 6 7 9 9 9]
[ 0 1 2 3 4 5 7 7 8 9 9 0]
[ 0 1 2 3 4 5 7 7 9 9 0 0]
[ 0 1 2 4 5 6 7 8 11 0 0 0]
[ 0 1 2 4 5 6 7 8 0 0 0 0]
[ 0 1 3 4 5 7 8 0 0 0 0 0]
[ 0 1 3 5 5 7 0 0 0 0 0 0]
[ 0 1 3 5 6 0 0 0 0 0 0 0]
[ 0 1 3 5 0 0 0 0 0 0 0 0]
[ 0 1 3 0 0 0 0 0 0 0 0 0]
[ 0 2 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0]]
"""
N = len(s)
s, psi = map(np.asarray, [s, psi])
s.shape = N
psi.shape = N
k = np.arange(1, N+1, dtype=np.double)
t = 1/(k + k0)
numer = s/psi
denom = 1/psi
S = np.zeros((N,N), dtype=s.dtype)
S[:,0] = s
t0 = t[:-1]
t1 = t[1:]
if ds is not None:
ds = np.asarray(ds)
if ds.shape is ():
# Convert relative error into absolute errors if ds is a
# number
ds = abs(s)*ds
ds.shape = N
if dpsi is None:
dpsi = 0*ds
else:
dpsi = np.asarray(dpsi)
dpsi.shape = N
d_psi = dpsi/abs(psi)**2 # d(1/psi)
if ds_psi is None:
ds_psi = np.sqrt(abs(ds/psi)**2 + abs(d_psi*s)**2)
else:
ds_psi = np.asarray(dpsi)
ds_psi.shape = N
dS = np.zeros((N, N), dtype=s.dtype)
dS[:, 0] = ds
dnumer = ds_psi
ddenom = d_psi
for n in xrange(1,N):
dt = (t1 - t0)
numer = np.diff(numer)/dt
denom = np.diff(denom)/dt
if ds is not None:
dnumer = dnumer[:-1] + abs((dnumer[1:] + dnumer[:-1])/dt)
ddenom = ddenom[:-1] + abs((ddenom[1:] + ddenom[:-1])/dt)
dS[:-n,n] = np.sqrt(abs(dnumer/denom)**2
+ abs(ddenom*numer/denom**2)**2)
t0 = t0[:-1]
t1 = t1[1:]
S[:-n,n] = numer/denom
if ds is None:
return S
else:
return S, dS
[docs]def levin(s, type='u', k0=0, degree=None, ds=None):
r"""Return the triangular acceleration table `a` for `n` terms of
the sequence `s[k]` using the Levin transformation `type`.
Parameters
----------
s : 1-d array of length `n`
`s[k]` should be the k'th term in the sequence.
type : `'t'`, `'u'`, or `'w'`
Type `'u'` is typically the most generally useful as it works
whenever type `'t'` works, however, `'t'` works well for sequences
:math:`s_k \sim r^k` and usually requires one less term, so for
these types of sequences, it should be used to reduce roundoff
error. For types `'t'` and `'u'` the returned table
acceleration table `a` is (`n-1` x `n-1`). For type `w`, the
table is (`n-2` x `n-2`).
k0 : float
Constant for the `'u'` transform.
degree : int
Degree of extrapolation to use to determine last difference `a_k`
ds : `None` or float or 1-d array of length `n`
This should be an estimate of the errors in the terms of the
sequence. If it is a single number, it is interpreted as the
relative error in the terms, otherwise it is taken to be the
array of absolute errors.
If it is not `None`, then the error estimate array `dS` will
also be returned.
Returns
-------
S : array (`n` x `n`)
Triangular acceleration array. The rows `S[k,:]` should
converge faster than the columns `S[:,n]`, otherwise the model
should not be trusted. Without round-off error, `S[0,-1]`
would be the best estimate, but at some point, round-off error
becomes significant.
dS : array (`n` x `n`)
If `ds` is provided and not `None`, then this array of error
estimates will also be returned.
"""
type = type.lower()
s = np.asarray(s)
a = np.diff(s)
s_ = s[:-1]
N = len(s_)
k = np.arange(1, N+1, dtype=np.double)
if 't' == type:
psi = a
elif 'u' == type:
psi = (k + k0)*a
elif 'w' == type:
s_ = s_[:-1]
da = np.diff(a)
a = a[:-1]
psi = a*a/da
dpsi = None
ds_psi = None
ds_ = ds
if ds is not None:
ds = np.asarray(ds)
if ds.shape is ():
# Convert relative error into absolute errors if ds is a
# number
ds = abs(s)*ds
ds.shape = len(s)
if 't' == type or 'u' == type:
# In this case psi ~ a and
# s_psi = s/psi ~ (1/(s_{k+1}/s_{k} - 1))
# ds_psi = ds_s*psi**2 where ds_s = d(s_{k+1}/s_{k})
dpsi = np.sqrt(ds[:-1]**2 + ds[1:]**2)
ds_s = np.sqrt(abs(ds[:-1]/s[1:])**2 +
abs(ds[1:]*s[:-1]/s[1:]**2)**2)
ds_psi = ds_s*abs(psi)**2
if 'u' == type:
# Include factor of (k + k_0)
dpsi *= (k + k0)
ds_psi /= (k + k0)
ds_ = ds[:-1]
elif 'w' == type:
# In this case psi ~ a^2/(da) and we approximate...
dp
# s_psi = s/psi ~ (1/(s_{k+1}/s_{k} - 1))
# ds_psi = ds_s*psi**2 where ds_s = d(s_{k+1}/s_{k})
da_ = np.sqrt(ds[:-1]**2 + ds[1:]**2)[:-1]
dda_ = np.sqrt(ds[:-2]**2 + (2*ds[1:-1])**2 + ds[2:]**2)
dpsi = np.sqrt(abs(2*da_*a/da)**2
+ abs(dda_*psi/da)**2)
ds_psi = None # Compute this from ds and dpsi
ds_ = ds[:-2]
return salzer(s=s_, psi=psi, ds=ds_, dpsi=dpsi, ds_psi=ds_psi)
def Aitken(s):
r"""Return the triangular acceleration table `a` for `n` terms of
the sequence `s[k]` using the Aitken transformation.
Parameters
----------
s : 1-d array of length `n`
`s[k]` should be the k'th term in the sequence.
Returns
-------
S : array (`n` x `n`)
Triangular acceleration array. The rows `S[k,:]` should
converge faster than the columns `S[:,n]`, otherwise the model
should not be trusted. Without round-off error, `S[0,-1]`
would be the best estimate, but at some point, round-off error
becomes significant.
"""
s = np.asarray(s).ravel()
N = len(s)
S = np.zeros((N, int(N/2)), dtype=s.dtype)
S[:,0] = s
for i in xrange(1,int(N/2)):
s = s[2:] - (s[2:] - s[1:-1])**2/(s[2:] - 2*s[1:-1] + s[:-2])
S[:-2*i,i] = s
return S
def _diff(s, k=2):
"""Return the array of differences `d` of the sequence `s` with
the last term extrapolated so that `len(s) == len(d)`.
Parameters
----------
s : array
Sequence
k : int
Order of extrapolation. The k'th difference of the sequence
will have a repeated last term.
Examples
--------
>>> s = [1, 4, 9, 16, 25]
>>> print _diff(s, k=1)
[ 3. 5. 7. 9. 9.]
>>> print _diff(s, k=2)
[ 3. 5. 7. 9. 11.]
"""
d = np.diff(np.asarray(s)).tolist()
d.append(np.polyval(np.polyfit(np.arange(k,0,-1), d[-k:],k-1),0))
return np.asarray(d)
try:
from pygsl.sum import levin_sum as _levin_sum
def levin_sum(s, n0=0):
"""Return an acceleration of the sequence `s` with an error
estimate using the levin summation technique.
Parameters
----------
's' : list
Initial sequence to be extrapolated.
'n0' : int, 'best', optional
Don't include the first `n0` terms in the extrapolation
portion. This is usefull if the first few terms do not
follow the decay pattern.
If `'best'`, then try truncaing the first few terms to
determine which tail will yeild the levin expansion with
the smallest error estimate. This makes the algorithm
`O(len(s)**2)`.
.. note:: The error estimate will not include possible
roundoff accumulation in the first `n0` terms.
Example
-------
>>> s = 1.0/(1+np.arange(30))**2/(np.pi**2/6.0)
>>> levin_sum(s)
(0.9999999999..., ...e-11)
>>> levin_sum(s, n0='best')
(0.9999999999..., ...e-11)
>>> s[:4] = 1, -1, 1, -1 # Mess up first few terms
>>> s[1] += 205/144/(np.pi**2/6.0)
>>> levin_sum(s)
(0.9999..., ...e-05)
>>> levin_sum(s, n0=4)
(0.999999999..., ...e-10)
>>> levin_sum(s, n0='best')
(0.999999999..., ...e-10)
"""
if 'best' == n0:
# Find the truncation point that minimizes the error
errs = [_levin_sum(s[n0:])[1] for n0 in xrange(len(s)-2)]
n0 = np.argmin(errs)
res, err = _levin_sum(s[n0:])
return np.sum(s[:n0]) + res, err
def levin_seq(s, n0=0):
"""Return an acceleration of the sequence `s` with an error
estimate using the levin summation technique.
Parameters
----------
's' : list
Initial sequence to be extrapolated.
'n0' : int, 'best', optional
Don't include the first `n0` terms in the extrapolation
portion. This is usefull if the first few terms do not
follow the decay pattern.
If `'best'`, then try truncaing the first few terms to
determine which tail will yeild the levin expansion with
the smallest error estimate. This makes the algorithm
`O(len(s)**2)`.
Example
-------
>>> s = 1.0 + 1.0/(1+np.arange(30))
>>> levin_seq(s)
(0.99999999999999..., ...e-15)
>>> levin_seq(s, n0='best')
(0.99999999999999..., ...e-15)
>>> s[:4] = 9, -199, 39, 4 # Mess up first few terms
>>> levin_seq(s)
(1.0003..., 0.0004...)
>>> levin_seq(s, n0='best')
(1.00000000000000..., ...e-14)
"""
ds = np.diff(s)
if 'best' == n0:
# Find the truncation point that minimizes the error
errs = [_levin_sum(ds[n0:])[1] for n0 in xrange(len(s)-2)]
n0 = np.argmin(errs)
res, err = _levin_sum(ds[n0:])
return s[n0] + res, err
except ImportError:
pass