2.5. Acceleration of Convergence#

2.5.1. Aitken’s \(\Delta^2\) Method#

Assume that \(\{c_n\}_{n=0}^{\infty}\) converges to a point \(c\) linearly. Since the rate of convergence is linear,

\[\begin{align*} \frac{c_{n+1} - c}{c_{n} - c} \approx \frac{c_{n+2} - c}{c_{n+1} - c}. \end{align*}\]

It follows that,

\[\begin{align*} (c_{n+1} - c)^2 \approx (c_{n+2} - c)(c_{n} - c). \end{align*}\]

Expanding both sides,

\[\begin{align*} & \quad c_{n+1}^2 -2 c_{n+1}c + c^2 \approx c_{n}c_{n+2} - (c_{n} + c_{n+2})c + c^2\\ \Rightarrow & \quad (c_{n+2} -2c_{n+1} + c_{n})c \approx c_{n}c_{n+2} - c_{n+1}^2\\ \Rightarrow & \quad c \approx \frac{c_{n}c_{n+2} - c_{n+1}^2}{(c_{n+2} -2c_{n+1} + c_{n})}. \end{align*}\]

Therefore,

\[\begin{align*} c & \approx \frac{c_{n}c_{n+2} - c_{n+1}^2}{(c_{n+2} -2c_{n+1} + c_{n})}\\ & = \frac{c_{n}c_{n+2} - c_{n+1}^2 +2c_{n} c_{n+1} + c_{n}^2 - 2c_{n} c_{n+1}- c_{n}^2}{(c_{n+2} -2c_{n+1} + c_{n})}\\ & = \frac{(c_{n}c_{n+2} - 2c_{n} c_{n+1} + c_{n}^2) - (c_{n+1}^2 -2c_{n} c_{n+1} + c_{n}^2)} {(c_{n+2} -2c_{n+1} + c_{n})}\\ & = \frac{c_{n}(c_{n+2} - 2 c_{n+1} + c_{n}) - (c_{n+1}^2 -2c_{n} c_{n+1} + c_{n}^2)}{(c_{n+2} -2c_{n+1} + c_{n})}\\ & = c_{n} - \frac{(c_{n+1} - c_{n})^2}{(c_{n+2} -2c_{n+1} + c_{n})} \end{align*}\]

  • Forward difference: \(\Delta c_{n} = c_{n+1} - c_{n}\)

  • Backward difference: \(\nabla c_{n} = c_{n} - c_{n-1}\)

\[\begin{align*} c \approx c_{n} - \frac{(\Delta c_{n})^2}{\Delta^2 c_{n}} \end{align*}\]

The following theorem is Theorem 2.14 from [Burden and Faires, 2005].

Theorem

Assume that \(\{c\}_{n}^{\infty}\) is a sequence that converges linearly to the limit \(c\) and that

(2.63)#\[\begin{equation} \lim_{n\to \infty} \frac{c_{n+1} - c}{c_{n} - c}<1 \end{equation}\]

The Aitken’s \(\Delta^2\) sequence \(\{\hat{c}\}_{n}^{\infty}\) converges to \(c\) faster than \(\{c\}_{n}^{\infty}\) in the sense that

(2.64)#\[\begin{equation} \lim_{n\to \infty} \frac{\hat{c}_{n+1} - c}{c_{n} - c} = 0. \end{equation}\]

2.5.2. Steffensen’s Method#

We can use Aitken’s \(\Delta^2\) method for finding fixed points and then increase the rate of convergence from linear to quadratic. This procedure is known as Steffensen’s Method.

We can accelerate the convergence of fixed-point iteration from linear to quadratic by a modification of Aitken’s \(\Delta^2\) method. This process is known as Steffensen’s method and it is slightly different from Aitken’s \(\Delta^2\) method. Aitken’s \(\Delta^2\) method generates the terms in order

\[\begin{align*} c^{0}_{0}, ~c^{0}_{1} = g\left(c^{0}_{0}\right), ~c^{0}_{2} = g\left(c^{0}_{1}\right), ~c^{1}_{0} = \{ \Delta^2 \}\left(c^{0}_{0}\right), ~c^{1}_{1} = g\left(c^{1}_{0}\right),\ldots \end{align*}\]

Observe that every third term of the Steffensen sequence is generated by

(2.65)#\[\begin{equation} \text{Every 3rd term }= c_{0}^{n} - \frac{(c_{1}^{n} - c_{0}^{n})^2}{(c_{2}^{n} -2c_{1}^{n} + c_{0}^{n})} \end{equation}\]

and the rest of the points are generated by fixed-point iteration on the previous term.

import pandas as pd 
import numpy as np

def Aitken_method(f, x0, TOL = 1e-4, Nmax = 100):
    '''
    Parameters
    ----------
    f : function
        DESCRIPTION.
    x0 : float
        DESCRIPTION. Initial estimate
    TOL : TYPE, optional
        DESCRIPTION. Tolerance (epsilon). The default is 1e-4.
    Nmax : Int, optional
        DESCRIPTION. Maximum number of iterations. The default is 100.

    Returns
    -------
    a dataframe

    '''
    xn = []
    xn.append(x0)
    En = [0]
    for n in range(0, Nmax-1):
        x1 = f(x0)
        x2 = f(x1)
        
        denominator = (x2 - x1) - (x1 - x0)
        if abs(denominator) < 10e-16:
            print('Denominator is too small')
            break
        
        newx = x2 - ( (x2 - x1)**2 ) / denominator
        xn.append(newx)
        
        if abs(newx - x2) < TOL:
            break
        
        x0 = newx
        
    En.extend(np.diff(xn))
    return pd.DataFrame({'xn': xn, 'En': En})

function [xn, En] = Aitken_method(f, x0, TOL, Nmax)
%{
Parameters
----------
f : function
    DESCRIPTION.
x0 : float
    DESCRIPTION. Initial estimate
TOL : TYPE, optional
    DESCRIPTION. Tolerance (epsilon). The default is 1e-4.
Nmax : Int, optional
    DESCRIPTION. Maximum number of iterations. The default is 100.

Returns
-------
a dataframe

Example
f =@(x) (1 - 0.5*x*x)
x0 = 1
%}

switch nargin
    case 3
        Nmax = 100;
    case 2
        TOL = 1e-4; Nmax = 100;
end

xn = zeros(Nmax, 1);
xn(1) = x0;
En = zeros(Nmax, 1);
En(1) = 0;
for n = 2:(Nmax-1)
    % x1
    x1 = f(x0);
    %
    x2 = f(x1);
    denominator = (x2 - x1) - (x1 - x0);

    if abs(denominator) < 10e-16
        disp('Denominator is too small')
        break
    end

    newx = x2 - ( (x2 - x1)^2 ) / denominator;

    En(n) = abs(newx - x2);
    xn(n) = newx;
    if abs(newx - x2) < TOL
        xn = xn(1:n);
        En = En(1:n);
        break
    end

    x0 = newx; 
end

Example: Find the fixed point of \(g(x) = (x^2 - 1)/3\) on \([-1,1]\) using the fixed-point method and Steffensen’s method with \(c_0 = 0\). The exact real root here is

\[\begin{equation*} c = \frac{3}{2} - \frac{\sqrt{13}}{2} \approx -0.3027756377319947. \end{equation*}\]

Solution:

Fixed-Point method: We can generate the entries of series \(\{c_{0}, c_{1}, c_{2}, \ldots\}\) using \(c_{n} = g(c_{n-1})\) and \(c_{0} = 0\). We have,

\[\begin{equation*} \begin{array}{lcl} c_{0} = 0 & \Rightarrow & c_{1} = g(c_{0}) = g(0) = -0.333333\\ c_{1} = -0.333333 & \Rightarrow & c_{2} = g(c_{1}) = g(-0.333333) = -0.296296 \\ c_{2} = -0.296296 & \Rightarrow & c_{3} = g(c_{2}) = g(-0.296296) = -0.304070 \\ c_{3} = -0.304070 & \Rightarrow & c_{4} = g(c_{3}) = g(-0.304070) = -0.302514 \\ \vdots & \vdots & \vdots \end{array} \end{equation*}\]
# This part is used for producing tables and figures
import sys
sys.path.insert(0,'..')
import hd_tools as hd
Loading BokehJS ... 
from hd_RootFinding_Algorithms import Fixed_Point
import numpy as np
g = lambda x : (x**2 - 1)/3
data = Fixed_Point(g = g, x0 = 0, TOL=1e-6, Nmax=20)
display(data.style.format({'fxn': "{:.4e}"}))
hd.it_method_plot(data, title = """Aitken’s Method""", ycol = 'En', ylabel = 'En)', color = 'OrangeRed')
  xn En
0 0.000000 1.000000
1 -0.333333 0.333333
2 -0.296296 0.037037
3 -0.304070 0.007773
4 -0.302514 0.001556
5 -0.302828 0.000315
6 -0.302765 0.000063
7 -0.302778 0.000013
8 -0.302775 0.000003
9 -0.302776 0.000001

The iterations was stopped when \(E_n < 10^{-6}\).

Steffensen’s method:

In this method, Every 3rd term is approximated as follows, and the rest of the points are approximated through the fixed-point method.

(2.66)#\[\begin{equation}\label{Steff_eq} \text{Every 3rd term }= c_{0}^{n} - \dfrac{(c_{1}^{n} - c_{0}^{n})^2}{(c_{2}^{n} -2c_{1}^{n} + c_{0}^{n})} \end{equation}\]

Note that, for the first step, set \(c_{0}^{0} = c_0\) (here \(c_0 = 0\)), and to approximate, \(c_{0}^{1}\), we need two additional points, \(c_{1}^{0}\) and \(c_{2}^{0}\). These points can be approximated through the fixed-point as follows,

\[\begin{align*} c_{1}^{0} &= g(c_{0}^{0}) = g(0) = -0.333333,\\ c_{2}^{0} &= g(c_{1}^{0}) = g(-0.333333) = -0.296296. \end{align*}\]

Now, the next approximation is done through \eqref{Steff_eq}:

\[\begin{equation*} c_{0}^{1}= c_{0}^{0} - \frac{(c_{1}^{0} - c_{0}^{0})^2}{(c_{2}^{0} -2c_{1}^{0} + c_{0}^{0})} = -0.300000. \end{equation*}\]

The next two point needs to be found through the fixed point, thus,

\[\begin{align*} c_{1}^{1} &= g(c_{0}^{1}) = g(-0.3) = -0.303333,\\ c_{2}^{1} &= g(c_{1}^{1}) = g(-0.303333) = -0.302663. \end{align*}\]

Similarly,

\[\begin{equation*} c_{0}^{2}= c_{0}^{1} - \frac{(c_{1}^{1} - c_{0}^{1})^2}{(c_{2}^{1} -2c_{1}^{1} + c_{0}^{1})} = -0.302775. \end{equation*}\]

Then, the next two points need to be found through the fixed point method,

\[\begin{align*} c_{1}^{2} &= g(c_{0}^{2}) = g(-0.302775) = -0.302776,\\ c_{2}^{2} &= g(c_{1}^{2}) = g(-0.302776) = -0.302776. \end{align*}\]

Constituting this procedure, we can present this in the following table.

Furthermore, using Steffensen’s method:

from hd_RootFinding_Algorithms import Steffensen_method
g = lambda x : (x**2 - 1)/3
data = Steffensen_method(g = g, c0 = 0, TOL = 1e-6, Nmax = 100)
display(data.style.format({'fxn': "{:.4e}"}))
  cn0 cn1 c2n En
0 0.000000 -0.333333 -0.296296 0.000000
1 -0.300000 -0.303333 -0.302663 0.300000
2 -0.302775 -0.302776 -0.302776 0.002775
3 -0.302776 nan nan 0.000000