Skip to content

Instantly share code, notes, and snippets.

@pv
Last active November 10, 2020 13:12
Show Gist options
  • Save pv/756d55c3ceed89291a2ab68302ffee51 to your computer and use it in GitHub Desktop.
Save pv/756d55c3ceed89291a2ab68302ffee51 to your computer and use it in GitHub Desktop.
Wynn epsilon table algorithm
# -*- coding:utf-8 -*-
from __future__ import division, absolute_import, print_function
import sys
import math
import collections
try:
from math import inf, nan
except ImportError:
inf = float('inf')
nan = float('nan')
eps = sys.float_info.epsilon
class WynnEpsilon(object):
r"""Wynn's epsilon algorithm.
Extrapolates the limit of the sequence::
{s0, s1, s2, ...}
Parameters
----------
sequence : iterable of array_like, optional
The sequence or its initial part.
max_rows : int, optional
Maximum number of table rows to retain.
Must be >= 2. Default: all
max_cols : int, optional
Maximum number of table columns to compute.
Default: all
irregular_tol : float, optional
Tolerance threshold for detecting irregular behavior
in the table.
Attributes
----------
max_rows
max_cols
irregular_tol
limit_estimate
error_estimate
epsilon
Methods
-------
append
extend
Notes
-----
The algorithm computes the even-epsilon table
.. math::
\begin{array}{c|ccccc}
p & k=-1 & k=0 & k=1 & k=2 & \cdots \\
\hline%& -------&-------&------------------&------------------&----------
0 & \infty & s_0 & & & \\
1 & \infty & s_1 & & & \\
2 & \infty & s_2 & \epsilon_2^{(0)} & & \\
3 & \infty & s_3 & \epsilon_2^{(1)} & \epsilon_4^{(0)} & \\
\vdots
\end{array}
Note that usually in literature this is expressed in a shifted
form, where the rows above are diagonals.
The values :math:`\epsilon` are associated with Padé approximants
related to the sequence. [1]_
Only last two rows are required for the recurrence relation::
e3 e0 ##
## e1 ##
## e2 new
d1 = e1 - e3
d2 = e2 - e1
d3 = e1 - e0
new = e1 + 1/(1/d1 + 1/d2 - 1/d3)
When divisors become small, this can signal convergence or
irregular behavior of the extrapolant. [2]_ Handling of these
situations in this implementation follows closely that in
QUADPACK/QELG. [3]_
References
----------
[1] P. Wynn, "Upon Systems of Recursions which Obtain among
the Quotients of the Pade Table,"
Numerische Mathematik 8 (3), 264 (1966).
[2] R. P. Eddy, "The Even-Rho and Even-Epsilon Algorithms
for Accelerating Convergence of a Numerical Sequence,"
DTNSRDC-81/083, David W Taylor Naval Ship Research and
development center (Bethesda MD, 1981).
[3] R. Piessens, E. de Doncker, QUADPACK (1983).
"""
def __init__(self, sequence=None, max_rows=None, max_cols=None, irregular_tol=1e-4):
self.limit_estimate = nan
self.error_estimate = nan
self._prev_results = collections.deque([], 3)
self.irregular_tol = irregular_tol
if max_rows is not None and max_rows != inf:
if not (max_rows >= 2):
raise ValueError("max_rows must be >= 2")
self.epsilon = collections.deque([], max_rows)
else:
self.epsilon = collections.deque([])
self.max_cols = max_cols
if sequence is not None:
self.extend(sequence)
def append(self, s):
"""
Add the next known element in the sequence.
"""
# Compute new elements
if len(self.epsilon) < 2:
n = 0
else:
n = min(len(self.epsilon[-1]), len(self.epsilon[-2]))
if self.max_cols is not None:
n = min(n, self.max_cols)
best_error = 0*s + inf
best_result = s
new_row = [s]
for k in range(n):
if k == 0:
e3 = inf
else:
e3 = self.epsilon[-2][k-1]
e0 = self.epsilon[-2][k]
e1 = self.epsilon[-1][k]
e2 = new_row[k]
d1 = e1 - e3 # = 1 / H_k^{p-1}
d2 = e2 - e1 # = 1 / V_k^p
d3 = e1 - e0 # = 1 / V_k^{p-1}
err1 = abs(d1)
tol1 = max(abs(e1), abs(e3)) * eps
err2 = abs(d2)
tol2 = max(abs(e2), abs(e1)) * eps
err3 = abs(d3)
tol3 = max(abs(e1), abs(e0)) * eps
if err2 <= tol2 and err3 <= tol3:
# Close to convergence
best_result = e2
best_error = err2 + err3
break
if (err1 <= tol1 and k > 0) or err2 <= tol2 or err3 <= tol3:
# Two elements close to each other: drop columns
break
ss = 1/d1 + 1/d2 - 1/d3
epsinf = abs(ss * e1)
if epsinf < self.irregular_tol:
# Irregular behavior in the table: drop columns
break
new = e1 + 1/ss
error = err2 + abs(new - e2) + abs(err3)
if error < best_error:
best_result = new
best_error = error
new_row.append(new)
self.epsilon.append(new_row)
self.limit_estimate = best_result
self._prev_results.append(best_result)
if len(self._prev_results) < 3:
self.error_estimate = inf
else:
self.error_estimate = max(abs(best_result - x)
for x in self._prev_results)
self.error_estimate = max(self.error_estimate,
eps*abs(self.limit_estimate))
def extend(self, terms):
for term in terms:
self.append(term)
def __str__(self):
ncols = max(len(row) for row in self.epsilon)
header = " ".join(["="*15]*ncols)
text = header + "\n"
text += "\n".join(" ".join("{:^15g}".format(x) for x in row)
for row in self.epsilon)
text += "\n" + header
text += "\n-> {} +/- {}".format(self.limit_estimate, self.error_estimate)
return text
def test_pi_limit():
import pytest
def s(n):
# s(n) -> 4 atan(1)
return sum((-1)**j * 4 / (2 * j + 1) for j in range(n + 1))
etab = WynnEpsilon()
etab.extend([s(j) for j in range(5)])
# Compare to [Graves-Morris et al., J. Comput. Appl. Math. 122, 51 (2000)]
# for diagonals.
expected = [[4.0],
[2.667],
[3.467, 3.167],
[2.895, 3.133],
[3.340, 3.145, 3.142]]
for j, row in enumerate(expected):
assert etab.epsilon[j] == pytest.approx(row, abs=1e-3)
assert etab.limit_estimate == pytest.approx(math.pi, abs=etab.error_estimate)
print(etab)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment