#!/usr/bin/env python
# coding: utf-8

# <style>
# @import url(https://www.numfys.net/static/css/nbstyle.css);
# </style>
# <a href="https://www.numfys.net"><img class="logo" /></a>
# 
# # Fixed-Point Iteration
# 
# ### Modules - Root Finding
# <section class="post-meta">
# By Eilif Sommer Øyre, Jonas Themsland and Jon Andreas Støvneng  
# </section>
# 
# Last edited: March 11th 2018
# 
# ---

# ## Introduction
# 
# An equation $f(x) = 0$ can always be written as $g(x) = x$. Fixed-point iteration can be applied to approximate the fixed number $r = g(r)$.
# 
# After selecting an initial guess $x_0$, the fixed-point iteration is $x_{n+1} = g(x_n)$. That is,
# 
# \begin{equation*}
# x_0 = \textrm{initial guess} \\
# x_1 = g(x_0) \\
# x_2 = g(x_1) \\
# x_3 = g(x_2) \\
# \vdots
# \end{equation*}
# 
# The algorithm is repeated until either a testing condition $e_{i+1} = \left|\:x_{i+1} - x_{i}\:\right| < \alpha$, where $\alpha$ is some tolerance limit is met, or until a fixed number of iterations $N$ is reached. Note that the method may or may not converge. Note that the choice of $g(x)$ is in general not unique.

# ## Example 1
# 
# As a simple introductory example we use the fixed-point iteration to solve the equation
# 
# $$\frac{1}{2}\sin(x) - x + 1 = 0.$$
# 
# The most natural is to use $g(x)= \sin(x)/2 + 1$ and then use the fixed-point iteration to solve $g(x)=x$.

# In[5]:


import math

x = 1     # initial guess
N = 10    # iterations

for i in range(1, N):
    x = math.sin(x) + 1
    print("x%i:\t%.5f"%(i, x))


# ## Convergence
# 
# For the iteration scheme to return a fixed point $r$, it needs to converge. A criterion for convergence is that the error, $e_i = \left|\:r - x_i \:\right|$, decreases for each iteration step. This means that the change from $x_i$ to $x_{i+1}$ also has to decrease for each step. The convergence criterion is explained in the following theorem [1]:
# 
# > **Theorem 1** Assume that $g$ is continuously differentiable, that  $g(r) = r$, and that  $S = \left|\:g'(r)\:\right| < 1$. Then fixed-point iteration converges linearly with rate $S$ to the fixed point $r$ for initial guesses sufficiently close to $r$. 
# 
# Note that in the example above we have $|g'(x)|=\cos(x)/2 < 1$ which means that the method will converge for all initial guesses.
# 
# The convergence can in fact be of higher order. This is explained in the following theorem [2]:
# 
# > **Theorem 2** Assume that $g$ is $p$ times continuously differentiable and that the fixed-iteration converges to $r$ for some initial guess $x_0$. If $g'(r)=g''(r)=\cdots=g^{(p-1)}(r)=0$ and $g^{(p-1)}(r)=0$, then the order of convergence is $p$.
# 
# The convergence is said to be of order $p$ when $\lim_{i\to\infty}e_{i+1}/e_{i}^p = \text{constant}$. 
# 
# **Excercise:** Prove theorem 1 using the mean value theorem and the theorem 2 using Taylor's theorem.

# ## Example 2: Babylonian Method for Finding Square Roots
# 
# As mentioned in the introduction, the choice of $g(x)$ is far from unique. We will now use fixed-point iteration to estimate the square root of a real and positive number $a$. That is, we want to solve $f(x) = x^2-a^2=0$. Two natural choices for $g(x)$ are $g_1(x)=a/x$ and $g_2(x)=x$. However, none of these converges since $|g_1'(\sqrt a)| = |g_2'(\sqrt a)| =1$. We can choose the mean as $g(x)$:
# 
# $$g(x) = \frac{1}{2}\left(\frac{a}{x} + x\right).$$
# 
# In this case we have 
# 
# $$\left|\:g'(x)\:\right|_{x=\sqrt a}=\left|\:\frac{1}{2}(1-\frac{a}{x^2})\:\right|_{x=\sqrt a} = 0 < 1$$
# 
# The theorem above thus implies that method converges for some initial guess $x_0$. The method is in fact globally convergent due to the convexity of $f(x)$.

# In[1]:


a = 0.07    # square of root
x = 1       # initial guess
N = 10    # iterations

for i in range(1, N):
    x = (a/x + x)/2
    print("x%i:\t%.5f"%(i, x))


# This method can also be derived from Newton's method (see the section on Newton's method below). The method was however first used by the Babylonian people long before Newton [2].

# ## Example 3: Adding Stopping Conditions
# 
# Underneath follows an implementation of the Fixed-Point Iteration with a testing condition $\left|\:x_{i+1} - x_i\:\right| < \alpha$. There are other stopping criteria that may be relevant such as the backward error $ \left|\:f(x_a)- f(0)\right|$. We need an additional stopping criteria in case convergence fails. We therefore stop the computation if the a given number of iterations $N$ is reached.

# In[14]:


def FPI(g, x, maxit=100, alpha=1e-5):
    """ A function using fixed-point iteration to compute 
    an approximate solution to the fixed point problem g(r) = r.
    
    Arguments:
    g           callable function, g(x)
    x           float. The initial guess
    alpha       float. Error tolerance
    maxit       int. Maximum number of iterations
    
    Returns:
    array of each fixed-point approximation
    """
    
    result = [x]
    x_next = g(x)
    i = 0               # Counter
    e = abs(x_next - x)
    
    while e > alpha :
        if i > maxit:
            print("No convergence.")
            break
        x = x_next
        x_next = g(x_next)
        result.append(x)
        if x - x_next > e:
            print("Divergence.")
            break
        e = abs(x_next - x)
        i += 1
    return result


# We want to solve the equation,
# 
# \begin{equation*}
# f(x) = x^4 + x - 1 = 0
# %\label{eq:example}
# \end{equation*}
# 
# using fixed point iteration. First, we rewrite the equation to the form, $x = g(x)$. This particular equation may be rewritten in three ways:
# 
# $$g_1(x) =    x = 1 - x^4,$$
# 
# $$g_2(x) =    x = \sqrt[4]{1 - x},$$
# 
# $$g_3(x) =    x = \frac{1 + x^4}{1 + 2x^3}.$$
# 
# Each of these has its own convergence rate. One of the solutions is clearly somewhere between 0 and 1, so we will use $x_0=0.5$ as initial guess.

# In[13]:


x_0 = 0.5
N = 100
alpha = 1e-5

def g1(x): return 1 - x**4
def g2(x): return (1 - x)**(1/4)
def g3(x): return (1 + 3*x**4)/(1 + 4*x**3)

def print_result(result):
    for i in range(len(result)):
        print("x%i:\t%.5f"%(i, result[i]))


# In[10]:


result_g1 = FPI(g1, x_0, N, alpha)
print_result(result_g1)


# In[11]:


result_g2 = FPI(g2, x_0, N, alpha)
print_result(result_g2)


# In[12]:


result_g3 = FPI(g3, x_0, N, alpha)
print_result(result_g3)


# The method does not converge for $g_1$, but both $g_2$ and $g_3$ converges towards the fixed number $0.72449$, which is a solution of the equation $f(x) = 0$. The method converges much faster for $g_3$ than for $g_2$.
# 
# This result can be verified by calculating the convergence rates $S = \left|\: g'(r) \:\right|$:
# 
# \begin{equation*}
# \left|\:g_1'(0.72449)\:\right| \approx 1.521 > 1, \\
# \left|\:g_2'(0.72449)\:\right| \approx 0.657 < 1, \\
# \left|\:g_3'(0.72449)\:\right| \approx 0.000 < 1. \\
# \end{equation*}
# 
# Theorem 1 tells us that the method will converge for $g_2$ and $g_3$ for some inital guess. Moreover, from theorem 2 we know that the method converges linearly for $g_2$ and the order of convergence for $g_3$ is two.
# 
# **Exercise:** Verify that the convergence for $g_2$ is linear and that the convergence of $g_3$ is of second order by using the numerical results above.

# ## In comparison with the Newton-Rhapson method
# 
# The Newton-Rhapson method is a special case of fixed-point iteration where the convergence rate is zero - the fastest possible. The method estimates the root of a differentiable function $f(x)$ by iteratively calculating the expression
# 
# \begin{equation*}
# x_{n+1} = x_n -\frac{f(x_n)}{f'(x_n)}.
# \end{equation*}
# 
# By using theorem 2, we can show that Newton's method is of second order (if $f'(r)\neq 0$). If $f''(r)=0$, then the method has a third order convergence.
# 
# This Newton iteration may be rewritten as a fixed point iteration
# 
# \begin{equation*}
# x_{n+1} = g(x_n),  \: n = 1, 2 ,3, ...
# \end{equation*}
# 
# where
# 
# \begin{equation*}
# r = g(r) = r - \frac{f(r)}{f'(r)}
# %\label{eq:newton}
# \end{equation*}
# 
# Provided that the iteration converges to a fixed point $r$ of $g$. From this we obtain
# 
# $$ \frac{f(r)}{f'(r)} = 0 \: \Rightarrow \: f(r) = 0$$
# 
# and thus the fixed point $r$ is a root of $f$.
# 
# The Newton's method is further discussed in the module [Root Finding - Newton-Rhapson Method](https://nbviewer.jupyter.org/urls/www.numfys.net/media/notebooks/newton_raphson_method.ipynb).

# ## Final comments
# 
# We have now used the fixed-point iteration method to solve the equation $f(x) = 0$ by rewriting it as a fixed-point problem. This is a simple method and which is easy to implement, but unlike the [Bisection Method](https://nbviewer.jupyter.org/urls/www.numfys.net/media/notebooks/bisection_method.ipynb) it only converges if the initial guess is sufficiently close to the root $r$. It is in general only locally convergent. However, the convergence rate may, or may not, be faster than that of the Bisection Method, which is 1/2.
# 

# ## References
# 
# [1] Sauer, T.: Numerical Analysis international edition, second edition, Pearson 2014  
# [2] Gautschi, W.: Numerical Analysis, second edition, Birkhäuser 2012 (https://doi.org/10.1007/978-0-8176-8259-0)