Aitken's delta-squared process
In numerical analysis, Aitken's delta-squared process or Aitken extrapolation is a series acceleration method, used for accelerating the rate of convergence of a sequence. It is named after Alexander Aitken, who introduced this method in 1926.[1] Its early form was known to Seki Kōwa (end of 17th century) and was found for rectification of the circle, i.e. the calculation of π. It is most useful for accelerating the convergence of a sequence that is converging linearly.
Definition
Given a sequence [math]\displaystyle{ X = {(x_n)}_{n\in\N} }[/math], one associates with this sequence the new sequence [math]\displaystyle{ A X = {\left(\frac{x_n\,x_{n+2}-x_{n+1}^2}{x_n+x_{n+2}-2\,x_{n+1}}\right)}_{n\in\Z^*}, }[/math] which can, with improved numerical stability, also be written as [math]\displaystyle{ (A X)_n = x_n-\frac{(\Delta x_n)^2}{\Delta^2 x_n}, }[/math] or equivalently as [math]\displaystyle{ (A X)_n = x_{n+2} - \frac{(\Delta x_{n+1})^2}{\Delta^2 x_n} = x_{n+2} - \frac{(x_{n+2}-x_{n+1})^2}{(x_{n+2}-x_{n+1})-(x_{n+1}-x_{n})} }[/math] where [math]\displaystyle{ \Delta x_{n}={(x_{n+1}-x_{n})},\ \Delta x_{n+1}={(x_{n+2}-x_{n+1})}, }[/math] and [math]\displaystyle{ \Delta^2 x_n=x_n -2x_{n+1} + x_{n+2}=\Delta x_{n+1}-\Delta x_{n},\ }[/math] for [math]\displaystyle{ n = 0, 1, 2, 3, \dots }[/math].
Obviously, [math]\displaystyle{ A X }[/math] is ill-defined if [math]\displaystyle{ \Delta^2 x }[/math] contains a zero element, or equivalently, if the sequence of first differences has a repeating term.
From a theoretical point of view, if that occurs only for a finite number of indices, one could easily agree to consider the sequence [math]\displaystyle{ A X }[/math] restricted to indices [math]\displaystyle{ n \gt n_0 }[/math] with a sufficiently large [math]\displaystyle{ n_0 }[/math]. From a practical point of view, one does in general rather consider only the first few terms of the sequence, which usually provide the needed precision. Moreover, when numerically computing the sequence, one has to take care to stop the computation when rounding errors in the denominator become too large, where the Δ2 operation may cancel too many significant digits. (It would be better for numerical calculation to use [math]\displaystyle{ \Delta x_{n+1} - \Delta x_{n}\ = (x_{n+2}-x_{n+1})-(x_{n+1}-x_{n})\ }[/math] rather than [math]\displaystyle{ x_n - 2x_{n+1} + x_{n+2}\ }[/math].)
Properties
Aitken's delta-squared process is a method of acceleration of convergence, and a particular case of a nonlinear sequence transformation.
Convergence of [math]\displaystyle{ \ \{ x_n \}_{n=1}^\infty\ }[/math] to limit [math]\displaystyle{ \ \ell\ }[/math] is called "linear" if there is some number μ ∈ (0, 1) for which [math]\displaystyle{ \lim_{n\to \infty} \frac{|x_{n+1} - \ell|}{|x_n - \ell|} = \mu\ . }[/math] Which means that the distance between the sequence and its limit shrinks by nearly the same proportion on every step, and that rate of reduction becomes closer to being constant with every step. (This is also called "geometric convergence"; this form of convergence is common for power series.)
Aitken's method will accelerate the sequence [math]\displaystyle{ x_n }[/math] if [math]\displaystyle{ \ \lim_{n \to \infty}\frac{(A x)_n-\ell}{x_n-\ell} = 0\ . }[/math]
[math]\displaystyle{ A }[/math] is not a linear operator, but a constant term drops out, viz: [math]\displaystyle{ \ A[x-\ell] = Ax - \ell\ , }[/math] if [math]\displaystyle{ \ell }[/math] is a constant. This is clear from the expression of [math]\displaystyle{ \ Ax\ }[/math] in terms of the finite difference operator [math]\displaystyle{ \ \Delta\ . }[/math]
Although the new process does not in general converge quadratically, it can be shown that for a fixed point process, that is, for an iterated function sequence [math]\displaystyle{ \ x_{n+1} = f(x_n)\ }[/math] for some function [math]\displaystyle{ \ f\ , }[/math] converging to a fixed point, the convergence is quadratic. In this case, the technique is known as Steffensen's method.
Empirically, the A-operation eliminates the "most important error term". One can check this by considering a sequence of the form [math]\displaystyle{ x_n=\ell+a^n+b^n }[/math], where [math]\displaystyle{ 0 \lt b \lt a \lt 1 }[/math]: The sequence [math]\displaystyle{ Ax }[/math] will then go to the limit like [math]\displaystyle{ b^n }[/math] goes to zero.
Geometrically, the graph of an exponential function [math]\displaystyle{ f(t) }[/math] that satisfies [math]\displaystyle{ f(n) = x_n }[/math], [math]\displaystyle{ f(n+1) = x_{n+1} }[/math] and [math]\displaystyle{ f(n+2)=x_{n+2} }[/math] has an horizontal asymptote at [math]\displaystyle{ \frac{x_{n} x_{n+2}-x_{n+1}^2}{x_{n}-2x_{n+1}+x_{n+2}} }[/math] (if [math]\displaystyle{ x_{n} - 2x_{n+1} + x_{n+2} \neq 0 }[/math]).
One can also show that if [math]\displaystyle{ x }[/math] goes to its limit [math]\displaystyle{ \ell }[/math] at a rate strictly greater than 1, [math]\displaystyle{ Ax }[/math] does not have a better rate of convergence. (In practice, one rarely has e.g. quadratic convergence which would mean over 30 resp. 100 correct decimal places after 5 resp. 7 iterations (starting with 1 correct digit); usually no acceleration is needed in that case.)
In practice, [math]\displaystyle{ Ax }[/math] converges much faster to the limit than [math]\displaystyle{ x }[/math] does, as demonstrated by the example calculations below. Usually, it is much cheaper to calculate [math]\displaystyle{ Ax }[/math] (involving only calculation of differences, one multiplication and one division) than to calculate many more terms of the sequence [math]\displaystyle{ x }[/math]. Care must be taken, however, to avoid introducing errors due to insufficient precision when calculating the differences in the numerator and denominator of the expression.
Example calculations
Example 1: The value of [math]\displaystyle{ \sqrt{2} \approx 1.4142136 }[/math] can be approximated by assuming an initial value for [math]\displaystyle{ a_0 }[/math] and iterating the following: [math]\displaystyle{ a_{n+1} = \frac{a_n + \frac{2}{a_n}}{2}. }[/math] Starting with [math]\displaystyle{ a_0 = 1: }[/math]
| n | x = iterated value | Ax |
|---|---|---|
| 0 | 1 | 1.4285714 |
| 1 | 1.5 | 1.4141414 |
| 2 | 1.4166667 | 1.4142136 |
| 3 | 1.4142157 | -- |
| 4 | 1.4142136 | -- |
It is worth noting here that Aitken's method does not save two iteration steps; computation of the first three Ax values required the first five x values. Also, the second Ax value is decidedly inferior to the 4th x value, mostly due to the fact that Aitken's process assumes linear, rather than quadratic, convergence.[citation needed]
Example 2: The value of [math]\displaystyle{ \frac{\pi}{4} }[/math] may be calculated as an infinite sum:
[math]\displaystyle{ \frac{\pi}{4} = \sum_{n=0}^\infty \frac{(-1)^n}{2n+1} \approx 0.785398 }[/math]
| n | term | x = partial sum | Ax |
|---|---|---|---|
| 0 | 1 | 1 | 0.79166667 |
| 1 | −0.33333333 | 0.66666667 | 0.78333333 |
| 2 | 0.2 | 0.86666667 | 0.78630952 |
| 3 | −0.14285714 | 0.72380952 | 0.78492063 |
| 4 | 0.11111111 | 0.83492063 | 0.78567821 |
| 5 | −9.0909091×10−2 | 0.74401154 | 0.78522034 |
| 6 | 7.6923077×10−2 | 0.82093462 | 0.78551795 |
| 7 | -6.6666667×10−2 | 0.75426795 | -- |
| 8 | 5.8823529×10−2 | 0.81309148 | -- |
In this example, Aitken's method is applied to a sublinearly converging series, accelerating convergence considerably. It is still sublinear, but much faster than the original convergence: the first Ax value, whose computation required the first three x values, is closer to the limit than the eighth x value.
Example pseudocode for Aitken extrapolation
The following is an example of using the Aitken extrapolation to help find the limit of the sequence [math]\displaystyle{ x_{n+1} = f(x_n) }[/math] when given some initial [math]\displaystyle{ x_0, }[/math] where the limit of this sequence is assumed to be a fixed point [math]\displaystyle{ f }[/math] (say [math]\displaystyle{ \alpha = f(\alpha) }[/math]). For instance, if the sequence is given by [math]\displaystyle{ x_{n+1} = \frac{1}{2} \left(x_n + \frac{2}{x_n}\right) }[/math] with starting point [math]\displaystyle{ x_0 = 1, }[/math] then the function will be [math]\displaystyle{ f(x) := \frac{1}{2}\left(x + \frac{2}{x}\right), }[/math] which has [math]\displaystyle{ \alpha := \sqrt{2} }[/math] as a fixed point (see Methods of computing square roots); it is this fixed point whose value will be approximated.
This pseudo code also computes the Aitken approximation to [math]\displaystyle{ f^{\prime}(\alpha) }[/math]. The Aitken extrapolates will be denoted by aitkenX. During the computation of the extrapolate, it is important to check if the denominator becomes too small, which could happen if we already have a large amount of accuracy; without this check, a large amount of error could be introduced by the division. This small number will be denoted by epsilon. Because the binary representation of the fixed point could be infinite (or at least too large to fit in the available memory), the calculation will stop once the approximation is within tolerance of the true value.
%These choices depend on the problem being solved
x0 = 1 %The initial value
f(x) = (1/2)*(x + 2/x) %The function that finds the next element in the sequence
tolerance = 10^-10 %10 digit accuracy is desired
epsilon = 10^-16 %Do not divide by a number smaller than this
maxIterations = 20 %Do not allow the iterations to continue indefinitely
haveWeFoundSolution = false %Were we able to find the solution to within the desired tolerance? not yet
for i = 1 : maxIterations
x1 = f(x0)
x2 = f(x1)
if (x1 ~= x0)
lambda = absoluteValue((x2 - x1)/(x1 - x0)) %OPTIONAL: Computes an approximation of |f'(fixedPoint)|, which is denoted by lambda
end
denominator = (x2 - x1) - (x1 - x0);
if (absoluteValue(denominator) < epsilon) %To avoid greatly increasing error, do not divide by too small of a number
print('WARNING: denominator is too small')
break %Leave the loop
end
aitkenX = x2 - ( (x2 - x1)^2 )/denominator
if (absoluteValue(aitkenX - x2) < tolerance) %If the value is within tolerance
print("The fixed point is ", aitkenX)) %Display the result of the Aitken extrapolation
haveWeFoundSolution = true
break %Done, so leave the loop
end
x0 = aitkenX %Update x0 to start again
end
if (haveWeFoundSolution == false) %If we were not able to find a solution to within the desired tolerance
print("Warning: Not able to find solution to within the desired tolerance of ", tolerance)
print("The last computed extrapolate was ", aitkenX)
end
See also
- Rate of convergence
- Limit of a sequence
- Fixed point iteration
- Richardson extrapolation
- Sequence transformation
- Shanks transformation
- Steffensen's method
Notes
- ↑ "On Bernoulli's numerical solution of algebraic equations". Proceedings of the Royal Society of Edinburgh 46: 289–305. 1926. doi:10.1017/S0370164600022070.
References
- William H. Press, et al., Numerical Recipes in C, (1987) Cambridge University Press, ISBN:0-521-43108-5 (See section 5.1)
- Abramowitz and Stegun, Handbook of Mathematical Functions, section 3.9.7
- Kendall E. Atkinson, An Introduction to Numerical Analysis, (1989) John Wiley & Sons, Inc, ISBN:0-471-62489-6
 |  |
