Iteration is a fundamental principle in computer science. As the name suggests, it is a process that is repeated until an answer is achieved or stopped. In this section, we study the process of iteration using repeated substitution.

More specifically, given a function g defined on the real numbers with real values and given a point x0 in the domain of g, the fixed point iteration is

\[ x_{i+1} = g(x_i ) \quad i =0, 1, 2, \ldots , \]
which gives rise to the sequence \( \left\{ x_i \right\}_{i\ge 0} . \) If this sequence converges to a point x, then one can prove that the obtained x is a fixed point of g, namely,
\[ x = g(x ) . \]

Example. Consider the convergent iteration

\[ p_0 = 0.5 \qquad \mbox{and} \qquad p_{k+1} = e^{-2p_k} \quad \mbox{for} \quad k=0,1,2,\ldots . \]
The first 10 terms are obtained by the calculations
\begin{align*} p_1 &= e^{-1} \approx 0.367879 , \\ p_2 &= e^{-2*p_1} \approx 0.479142 , \\ p_3 &= e^{-2*p_2} \approx 0.383551 , \\ \vdots & \qquad \vdots \\ p_9 &= e^{-2*p_8} \approx 0.409676 , \\ p_{10} &= e^{-2*p_9} \approx 0.440717 . \\ \end{align*}
The sequence is converging, and Mathematica reveals that
\[ \lim_{k\to \infty} p_k = 0.426302751 \ldots . \]
p[0] = 0.5
Do[Print["k= ", k, " " , "p= " , p[k] = Exp[-2*p[k - 1]]], {k, 1, 10}]
FindRoot[p == Exp[-2*p], {p, 0.5}]

Example. Consider \( g(x) = 10/ (x^3 -1) \) and the fixed point iterative scheme \( x_{i+1} = 10/ (x^3_i -1) ,\) with the initial guess x0 = 2.

If we repeat the same procedure, we will be surprised that the iteration is gone into an infinite loop without converging.

Theorem: Assume that the function g is continuous on the interval [a,b].

  • If the range of the mapping y = g(x) satisfies \( y \in [a,b] \) for all \( x \in [a,b] , \) then g has a fixed point in [a,b].
  • Furthermore, suppose that the derivative g'(x) is defined over (a,b) and that a positive constant (called Lipschitz constant) K < 1 exists with \( |g' (x) | \le K < 1 \) for all \( x \in (a,b) , \) then g has a unique fixed point P in [a,b].

Theorem: Assume that the function g and its derivative are continuous on the interval [a,b]. Suppose that \( g(x) \in [a,b] \) for all \( x \in [a,b] , \) and the initial approximation x0 also belongs to the interval [a,b].

  • If \( |g' (x) | \le K < 1 \) for all \( x \in (a,b) , \) then the iteration \( x_{i+1} = g(x_i ) \) will converge to the unique fixed point \( P \in [a,b] . \) In this case, P is said to be an attractive fixed point: \( P= g(P) . \)
  • If \( |g' (x) | > 1 \) for all \( x \in (a,b) , \) then the iteration \( x_{i+1} = g(x_i ) \) will not converge to P. In this case, P is said to be a repelling fixed point and the iteration exhibits local divergence. ■
In practice, it is often difficult to check the condition \( g([a,b]) \subset [a,b] \) given in the previous theorem. We now present a variant of it.

Theorem: Let P be a fixed point of g(x), that is, \( P= g(P) . \) Suppose g(x) is differentiable on \( \left[ P- \varepsilon , P+\varepsilon \right] \quad\mbox{for some} \quad \varepsilon > 0 \) and g(x) satisfies the condition \( |g' (x) | \le K < 1 \) for all \( x \in \left[ P - \varepsilon , P+\varepsilon \right] . \) Then the sequence \( x_{i+1} = g(x_i ) , \) with starting point \( x_0 \in \left[ P- \varepsilon , P+\varepsilon \right] , \) converges to P. ■

Remark: If g is invertible then P is a fixed point of g if and only if P is a fixed point of g-1.

Remark: The above therems provide only sufficient conditions. It is possible for a function to violate one or more of the hypotheses, yet still have a (possibly unique) fixed point.

Example: The function \( g(x) = 2x\left( 1-x \right) \) violates the hypothesis of the theorem because it is continuous everywhere \( (-\infty , \infty ) . \) Indeed, g(x) clearly does not map the interval \( [0.5, \infty ) \) into itself; moreover, its derivative \( |g' (x)| \to + \infty \) for large x. However, g(x) has fixed points at x = 0 and x = 1/2.

Example: Consider the equation

\[ x = 1 + 0.4\, \sin x , \qquad \mbox{with} \quad g(x) = 1 + 0.4\, \sin x . \]
Note that g(x) is a continuous functions everywhere and \( 0.6 \le g(x) \le 1.4 \) for any \( x \in \mathbb{R} . \) Its derivative \( \left\vert g' (x) \right\vert = \left\vert 0.4\,\cos x \right\vert \le 0.4 < 1 . \) Therefore, we can apply the theorem and conclude that the fixed point iteration
\[ x_{k+1} = 1 + 0.4\, \sin x_k , \qquad k=0,1,2,,\ldots \]
will converge for arbitrary initial point from the interval [0.6 , 1.4].

Example: Consider a similar equation

\[ x = 1 + 2\, \sin x , \qquad \mbox{with} \quad g(x) = 1 + 2\, \sin x . \]
Since \( -1 \le g(x) \le 3 , \) we are looking for a fixed point from this interval, [-1,3]. Since for its derivative we have
\[ g'(x) = 2\, \cos x \qquad \Longrightarrow \qquad \max_{x\in [-1,3]} \,\left\vert g' (x) \right\vert =2 > 1, \]
we do not expect convergence of the fixed point iteration \( x_{k+1} = 1 + 2\, \sin x_k . \)
function [k, p, err, P] = fixpt(g,p0,tol,max1)
% Input:   g is the iteration function input as a sring 'g' 
%          p0 is the initial guess for the fixed point
%          tol is the tolerance 
%          max1 is the maximum number of iterations
% Output:  k is the number of iterations that were carried out 
%          p is the approximation to the fixed point 
%          err is the error in the approximation
%          P contains the sequence {pn}
P(1)=p0;
for k=2:max1 
   P(k)=feval(g,P(k-1)); 
   err=abs(P(k)-P(k-1)); 
   relerr=err/(abs(P(k))+eps);
   p=P(k);
   if (err<: tol) | (relerr< tol), break;end
   end
if k==max1;
   disp('maximum number of iterations exceeded')
end
P=P';
Another matlab function:
% Fixed-Point Iteration Numerical Method for finding the x root of f(x) to make f(x) = 0

function [xR,err,n,xRV,errV,AFD1,AFD2] = FixedPointNM(AF,xi,ed)
% Inputs: with examples
% AF = anonymous function equation: AF = @(x) 1-((20^2)./(9.81*(((3*x)+((x.^2)/2)).^3))).*(3+x);
% xi = initial guess x = xR, where xR = x root: xi = 0.5;
% ed = desired approximate relative error = |(current - previous)/current|: ed = 0.01;
% Outputs
% xR = x root
% err = approximate relative error
% n = number of iterations
% xRV = x root vector
% errV = approximate relative error vector
% AFD1 = anonymous function 1st derivative
% AFD2 = anonymous function 2nd derivative
%% Computations
% Derivatives
AFSYM = sym(AF); % symbolic math function
AFD1SYM = diff(AFSYM); % symbolic math function 1st derivative
AFD2SYM = diff(AFD1SYM); % symbolic math function 2nd derivative
AFD1 = matlabFunction(AFD1SYM); % anonymous function 1st derivative
AFD2 = matlabFunction(AFD2SYM); % anonymous function 2nd derivative
% Fixed-Point Iteration Method
err = 1; % initial approximate relative error = 100%
k = 1; % initial counter
while ed < err % compares desired versus calculated error
    xR = xi+AF(xi); % x root using fixed point iteration method
    xRV(k+1) = xR; % stores the x root per iteration in a vector
    err = abs((xRV(k+1) - xRV(k))/xRV(k+1)); % approximate relative error
    errV(k+1) = err; % stores the error into a vector
    xi = xR; % new guess x = xR, where xR = x root
    k = k+1; % increase counter
end
n = k - 1; % number of iterations
xRV = xRV(2:end); % readjust x root vector
errV = errV(2:end); % readjust approximate relative error vector
Another code:
function FixedPointIteration
%Fixed point iteration example.
%4x^2 = 10 - x^3 can be rewritten as x^2 = 0.25*(10 - x^3) 
%===> x = 0.5*(10- x^3)^0.5 if we take the postive root only.
%now we can use fixed point iteration.

i = 1;         
p0 = 1.5;           %our guess is p = 1.5.
error = 0.0000001;  %the precision we want to find our solution in.
N = 100;            %max number of iterations.

syms 'x' 
g(x) = 0.5*((10 - (x^3))^0.5); %i.e. what we are trying to solve

while i <= N
    p = g(p0);                  %the next p value
    if abs((p - p0)) < error    %stopping criterion using the required precision
        fprintf('Solution is %f \n', double(p))     %the output p is symbolic so need ot force it into a decimal
        return
    end
    
i = i + 1;   
p0 = p;         %update the value of i and p to continue the loop with

end
fprintf('Solution did not converge to within precision = %d in %d iterations \n', error, N)
end
 

 

II. Aitken's Algorithm


This algorithm was proposed by one of New Zealand's greatest mathematicians Alexander Craig "Alec" Aitken (1895--1967). It consists in the following stages:
  • select x0;
  • calculate
    \[ x_1 = g(x_0 ) , \qquad x_2 = g(x_1 ) ; \]
  • calculate
    \[ x_3 = x_2 + \frac{\lambda_2}{1- \lambda_2} \left( x_2 - x_1 \right) , \qquad \mbox{where} \quad \lambda_2 = \frac{x_2 - x_1}{x_1 - x_0} ; \]
  • calculate
    \[ x_4 = g(x_3 ) , \qquad x_5 = g(x_4 ) ; \]
  • calculate x6 as the extrapolate of \( \left\{ x_3 , x_4 , x_5 \right\} . \) Continue this procedure, ad infinatum.

Theorem (Aitken's Acceleration): Assume that the sequence \( \{ p_n \}_{n\ge 0} \) converges linearly to the limit p  and that \( p_n \ne p \) for all \( n \ge 0. \) If there exists a real number A < 1 such that  

\[ \lim_{n \to \infty} \, \frac{p- p_{n+1}}{p- p_n} =A, \]
then the sequence \( \{ q_n \}_{n\ge 0} \) defined by  
\[ q_n = p_n - \frac{\left( \Delta p_n \right)^2}{\Delta^2 p_n} = p_n - \frac{\left( p_{n+1} - p_n \right)^2}{p_{n+2} - 2p_{n+1} + p_n} \]
converges to p faster than \( \{ p_n \}_{n\ge 0} \) in the sense that \( \lim_{n \to \infty} \, \left\vert \frac{p - q_n}{p- p_n} \right\vert =0 . \) ■ 

Mathworks code:
% function f1=aitken(x,f,x1)
%  f1= aitken(x,f,x1)
%  f = corresponding function of x
%  f1= corresponding function of x1
   
    if ~isequal(length(x),length(f))
        disp('Error: Incomplete data points')
        f1=[];
        return
    end
    
    n=length(x); %length of available data
    D=x1-x;
    A(:,1)=f;
    for i=2:n
        for j=2:i
            A(i,j)=(D(j-1)*A(i,j-1)-D(i)*A(j-1,j-1))/(x(i)-x(j-1));
        end
    end
        
    f1=A(n,n);

 

III. Steffensen's Acceleration


Johan Frederik Steffensen (1873--1961) was a Danish mathematician, statistician, and actuary who did research in the fields of calculus of finite differences and interpolation. He was professor of actuarial science at the University of Copenhagen from 1923 to 1943. Steffensen's inequality and Steffensen's iterative numerical method are named after him.

When Aitken's process is combined with the fixed point iteration in Newton's method, the result is called Steffensen's acceleration.  Starting with p0, two steps of Newton's method are used to compute \( p_1 = p_0 - \frac{f(p_0 )}{f' (p_0 )} \) and \( p_2 = p_1 - \frac{f(p_1 )}{f' (p_1 )} , \) then Aitken's process is used to compute \( q_0 = p_0 -  \frac{\left( \Delta p_0 \right)^2}{\Delta^2 p_0} . \)  To continue the iteration set \( q_0 = p_0 \) and repeat the previous steps. This means that we have a fixed-point iteration:  

\[ p_0 , \qquad p_1 = g(p_0 ), \qquad p_2 = g(p_1 ). \]
Then we compute
\[ q_0 = p_0 -  \frac{\left( \Delta p_0 \right)^2}{\Delta^2 p_0} = p_0 - \frac{\left( p_1 - p_0 \right)^2}{p_2 - 2p_1 +p_0} . \]
At this point, we ``restart'' the fixed point iteration with \( p_0 = q_0 , \) e.g.,
\[ p_3 = q_0 , \qquad p_4 = g(p_3 ), \qquad p_5 = g(p_4 ). \]
and compute
\[ q_3 = p_3 -  \frac{\left( \Delta p_3 \right)^2}{\Delta^2 p_3} = p_3 - \frac{\left( p_4 - p_3 \right)^2}{p_5 - 2p_4 +p_3} . \]
If at some point \( \Delta^2 p_n =0 \) (which appears in the denominator), then we stop and select the current value of pn+2 as our approximate answer.

Steffensen's acceleration is used to quickly find a solution of the fixed-point equation x = g(x) given an initial approximation p0. It is assumed that both g(x) and its derivative are continuous, \( | g' (x) | < 1, \) and that ordinary fixed-point iteration converges slowly (linearly) to p.

function [p,Q] = steff(g,dg,p0,delta,epsilon,max1)
% Input:   g is the iteration function input as a sring 'g' 
%          dg is the derivative of g input as a string 'dg'
%          p0 is the initial guess for the fixed point
%          delta is the tolerance for p0
%          epsilon is the tolerance for the function values y
%          max1 is the maximum number of iterations
% Output:  Q is the matrix containing the Steffensen sequence 
%          p is the Steffensen approximation to the fixed point 

% Initialize the matrix R
R=zeros(max1,3);
R(1,1)=p0;
for k=1:max1 
    for j=2:3 
       % denominator in Newton method is calculated
       nrdenom=feval(dg,R(k,j-1));
       
       % calculate Newton approximations
       if nrdenom==0
          'division by zero in Newton method'
          break
       else
          R(k,j)=R(k,j-1)-feval(g,R(k,j-1))/nrdenom;
       end
    
    % denominator in Aitken's acceleration process 
       aadenom=R(k,3)-2*R(k,2)+R(k,1);
       
       % calculate Aitken's acceleration approximations
       if aadenom==0
          'division by zero in Aitken's acceleration'
          break
       else
          R(k+1,1)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;
       end
    
    end
    % end program if division by zero occured
    if (nrdenom==0)|(aadenom==0)  
       break
    end
      
    % stopping criteria are evaluated
    err=abs(R(k,1)-R(k+1,1));
    relerr=err/abs(R(k+1,1)+delta);
    y=feval(g,R(k+1,1)) ;
    if (err < delta)|(relerr < delta)|(y < epsilon)
       % p and the matrix Q are determined
       p=R(k+1,1);
       Q=R(1:k+1,:);
       break
    end
end

 

V. Wegstein's Algorithm


To solve the fixed point problem \( g(x) =x , \) the following algorithm, presented in
J.H. Wegstein, Accelerating convergence of iterative processes, Communications of the Association for the Computing Machinery (ACM), 1, 9--13, 1958.
can be used
\[ x_{k+1} = \frac{x_{k-1} g(x_k ) - x_k g(x_{k-1})}{g(x_k ) + x_{k-1} -x_k - g(x_{k-1})} , \quad k=1,2,\ldots . \]
Till some extend, Wegstein's method is a version of predictor-corrector method because it can be broken into two stages:
\[ \begin{split} p^{(n+1)} = g \left( x^{(n)} \right) , \quad x^{(n+1)} = q\, x^{(n)} + \left( 1-q \right) p^{(n+1)} , \quad n=1,2,\ldots ; \\ q= \frac{b}{b-1} , \quad b= \frac{x^{(n)} - p^{(n+1)}}{x^{(n-1)} - x^{(n)}} , \end{split} \]
for \( n=0,1,2,\ldots . \) Wegstein (1958) found that the iteration converge differently depending on the value of q.
q < 0  Convergence is monotone 
0 < q < 0.5  Convergence is oscillatory 
0.5 < q < 1  Convegence 
q > 1  Divergence 
Wegstein's method depends on two first guesses x0 and x1 and poor start may cause the process to fail, converge to another root, or, at best, add unnecessary number of iterations.
function  Complete


Complete