import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
%matplotlib inline
from ipywidgets import interact
from ipywidgets import widgets
sym.init_printing()
from scipy import optimize
import pandas as pd
pd.set_option("display.colheader_justify","center")
pd.options.display.float_format = '{:.10f}'.format
from colorama import Fore, Back, Style
# https://pypi.org/project/colorama/
# conda install -c anaconda colorama
# Fore: BLACK, RED, GREEN, YELLOW, BLUE, MAGENTA, CYAN, WHITE, RESET.
# Back: BLACK, RED, GREEN, YELLOW, BLUE, MAGENTA, CYAN, WHITE, RESET.
# Style: DIM, NORMAL, BRIGHT, RESET_ALL
textBold = lambda x: Style.BRIGHT+x+Style.RESET_ALL
textBoldH = lambda x: Style.BRIGHT+Back.YELLOW+x+Style.RESET_ALL

"""
Explanation: <center>
    <img src="http://sct.inf.utfsm.cl/wp-content/uploads/2020/04/logo_di.png" style="width:60%">
    <h1> INF285 - Computación Científica </h1>
    <h2> Roots of 1D equations </h2>
    <h2> <a href="#acknowledgements"> [S]cientific [C]omputing [T]eam </a> </h2>
    <h2> Version: 1.37</h2>
</center>
<div id='toc' />

Table of Contents

Introduction
Bisection Method
Fixed Point Iteration and Cobweb diagram
FPI - example from etxtbook
Newton Method
Wilkinson Polynomial
Acknowledgements
Extra Examples
End of explanation
"""


f = lambda x: x+np.log(x)-3

"""
Explanation: <div id='intro' />

Introduction
Back to TOC
In this document we're going to study how to find roots of 1D equations using numerical methods. 
First, let's start with the definition of a root:
<b>Definition</b>: The function $f(x)$ has a <b>root</b> in $x = r$ if $f(r) = 0$.
Example: Let's say we want to solve the equation $r + \log(r) = 3$.
We can re-arrange the equation as follows: $r + \log(r) - 3 = 0$. 
Thus, solving the previous equation is equivalent to find the root of $f(x) = x + \log(x) - 3$.
This example shows how we can translate an equation into a root-finding problem!
We will study now several numerical methods to find roots.
We will start by defining a function $f(x)$ using a lambda definition.
End of explanation
"""


# Definition of symbolic variable
x = sym.Symbol('x')
# Defining 'symbolic' function
fsym = lambda x: x+sym.log(x)-3
# Finding the root 'symbolically' and obtaining the only root
r=sym.solve(sym.Eq(fsym(x), 0), x)[0]
print(textBoldH('Root obtained:'),r)
print(textBoldH('Numerical root:'),r.evalf())

"""
Explanation: Notice that we have used the NumPy implementation for the logarithmic function. Quick question: what is the base for this implementation? Is it the natural logarithm or logarithm base 10?
But before we start working on a numerical implementation, we should always consider in solving the problem algebraically.
This can be done with SymPy, i.e. using symbolic computation.
For instance, we will start by defining a symbolic variable:
End of explanation
"""


def find_root_manually(r=2.0):
    # Defining a vector to evaluate f(x) in a vectorized fashion
    x = np.linspace(1,3,1000)
    # reating the figure
    plt.figure(figsize=(8,8))
    # Plotting the function in a vectorized way
    plt.plot(x,f(x),'b-')
    # Plotting the x-axis. 
    # Quick question: Why do we have to multiply 'x' by '0'? What would happen if we only put '0' instead of 'x*0'?
    plt.plot(x,x*0,'r--')
    # Adding the background grid to the plot. We strongly recommend it!
    plt.grid(True)
    # Just adding labels.
    plt.ylabel('$f(x)$',fontsize=16)
    plt.xlabel('$x$',fontsize=16)
    plt.title('$r='+str(r)+',\, f(r)='+str(f(r))+'$',fontsize=16)
    plt.plot(r,f(r),'k.',markersize=20)
    plt.show()
interact(find_root_manually,r=(1,3,1e-3))

"""
Explanation: Lamber W function: https://en.wikipedia.org/wiki/Lambert_W_function
We will now obtain the root 'manually'.
This is not a recommended path but it is useful initially.
End of explanation
"""


def bisect(f, a, b, tol=1e-5, maxNumberIterations=100):
    # Evaluating the extreme points of the interval provided
    fa = f(a)
    fb = f(b)
    # Iteration counter.
    i = 0
    # Just checking if the sign is not negative => not root  necessarily 
    if np.sign(f(a)*f(b)) >= 0:
        print('f(a)f(b)<0 not satisfied!')
        return None
  
    # Output table to store the numerical evolution of the algorithm
    output_table = []
    
    # Main loop: it will iterate until it satisfies one of the two criterias:
    # The tolerance 'tol' is achived or the max number of iterations is reached.
    while ((b-a)/2 > tol) and i<=maxNumberIterations:
        # Obtaining the midpoint of the interval. Quick question: What could happen if a different point is used?
        c = (a+b)/2.
        # Evaluating the mid point
        fc = f(c)
        # Saving the output data
        output_table.append([i, a, c, b, fa, fc, fb, b-a])

        # Did we find the root?
        if fc == 0:
            print('f(c)==0')
            break
        elif np.sign(fa*fc) < 0:
            # This first case consider that the new inetrval is defined by [a,c]
            b = c
            fb = fc
        else:
            # This second case consider that the new interval is defined by [c,b]
            a = c
            fa = fc
        # Increasing the iteration counter
        i += 1
    
    # Showing final output table
    columns    = ['$i$', '$a_i$', '$c_i$', '$b_i$', '$f(a_i)$', '$f(c_i)$', '$f(b_i)$', '$b_i-a_i$']
    df = pd.DataFrame(data=output_table, columns=columns)
    display(df)
    
    # Computing the best approximation obtaind for the root, which is the midpoint of the final interval.
    xc = (a+b)/2.
    return xc

# Initial example
f1 = lambda x: x+np.log(x)-3
# A different function, notice that x is multiplied to the exponential now and not added, as before.
f2 = lambda x: x*np.exp(x)-3
# This is the introductory example about Fixed Point Iteration
f3 = lambda x: np.cos(x)-x
bisect(f1,1e-10,3) # Recall to change the 'tol'!

"""
Explanation: <div id='bisection' />

Bisection Method
Back to TOC
The bisection method finds the root of a function $f$. 
It requires that:
1. $f$ be a continuous function.
2. The interval $[a,b]$, such that $f(a)\cdot f(b) < 0$.
If these 2 conditions are satisfied, it means that there is a value $r$, between $a$ and $b$, for which $f(r) = 0$.
To summarize how this method works, start with the aforementioned interval (checking that there's a root in it), and split it into two smaller intervals $[a,c]$ and $[c,b]$. 
Then, check which of the two intervals contains a root. 
Keep splitting each "eligible" interval until the algorithm converges or the tolerance is achived.
End of explanation
"""


# Just plotting the Cobweb diagram: https://en.wikipedia.org/wiki/Cobweb_plot
def cobweb(x,g=None):
    min_x = np.amin(x)
    max_x = np.amax(x)
    
    plt.figure(figsize=(10,10))
    ax = plt.axes()
    plt.plot(np.array([min_x,max_x]),np.array([min_x,max_x]),'b-')
    for i in np.arange(x.size-1):
        delta_x = x[i+1]-x[i]
        head_length =  np.abs(delta_x)*0.04
        arrow_length = delta_x-np.sign(delta_x)*head_length
        ax.arrow(x[i], x[i], 0, arrow_length, head_width=1.5*head_length, head_length=head_length, fc='k', ec='k')
        ax.arrow(x[i], x[i+1], arrow_length, 0, head_width=1.5*head_length, head_length=head_length, fc='k', ec='k')
    
    if g!=None:
        y = np.linspace(min_x,max_x,1000)
        plt.plot(y,g(y),'r')
    
    plt.title('Cobweb diagram')
    plt.grid(True)
    plt.show()

# This code performs the fixed point iteration.
def fpi(g, x0, k, flag_cobweb=False):
    # This is where we store all the approximation, 
    # this is technically not needed but we store them because we need them for the cobweb diagram at the end.
    x = np.empty(k+1)
    # Just starting the fixed point iteration from the 'initial guess'
    x[0] = x0
    # Initializing the error in NaN
    error_i = np.nan
    
    # Output table to store the numerical evolution of the algorithm
    output_table = []
    
    # Main loop
    for i in range(k):
        # Iteration
        x[i+1] = g(x[i])
        # Storing error from previous iteration
        error_iminus1 = error_i
        # Computing error for current iteration.
        # Notice that from the theory we need to compute e_i=|x_i-r|, i.e. we need the root 'r'
        # but we don't have it, so we approximate it by 'x_{i+1}'.
        error_i = abs(x[i]-x[i+1])
        output_table.append([i,x[i],x[i+1],error_i,error_i/error_iminus1,error_i/(error_iminus1**((1+np.sqrt(5))/2.)),error_i/(error_iminus1**2)])
    
    # Showing final output table
    columns    = ['$i$', '$x_i$', '$x_{i+1}$', '$e_i$', r'$\frac{e_i}{e_{i-1}}$', r'$\frac{e_i}{e_{i-1}^\alpha}$', r'$\frac{e_i}{e_{i-1}^2}$']
    df = pd.DataFrame(data=output_table, columns=columns)
    display(df)
    
    # Just showing cobweb if required
    if flag_cobweb:
        cobweb(x,g)
    return x[-1]

# First example
g = lambda x: np.cos(x)
# Examples from classnotes
g1 = lambda x: -(3/2)*x+5/2
g2 = lambda x: -(1/2)*x+3/2
g3 = lambda x: -x+2
fpi(g1, 1.1, 5, True)

# Suggestions:
# 1.- A very useful and simple 'limit cicle' example! Try it. Credit: anonymous student from class 20210922.
# fpi(lambda x: -x, 2, 10, True)
# 2.- Try the next example. Why do we see 1.0000 for e_{i+1}/e_i over 90 iterations?
# fpi(g, 1, 100, True)
# 3.- The following fixed-point iteration obtain sqrt(2). Credits: Anonymous student from class 20210921 and 20210922.
# gD = lambda x: (x+2)/(x+1)
# fpi(gD, 1, 10, True)

"""
Explanation: It's very important to define a concept called convergence rate. 
This rate shows how fast the convergence of a method is at a specified point.
The convergence rate for the bisection is always 0.5 because this method uses the half of the interval for each iteration.
In this particular case we observe $e_{i+1} \approx \dfrac{e_{i}}{2}$, why? where?
<div id='fpi' />

Fixed Point Iteration and Cobweb diagram
Back to TOC
To learn about the Fixed-Point Iteration we will first learn about the concept of a Fixed Point.
A Fixed Point of a function $g$ is a real number $r$, where $g(r) = r$
The Fixed-Point Iteration is based in the Fixed Point concept and works like this to find the root of a function:
\begin{align} 
    x_{0} &= initial_guess \ 
    x_{i+1} &= g(x_{i})
\end{align}
To find an equation's root using this method, we'll have to rearrange the equation to make it of the following form $x = g(x)$.
For example, if we want to obtain the root of $f(r)=0$, one could add a zero convenient this way, $f(r)+r=r$, i.e. we add it $r$ on both sides.
This way we have $g(r)=r+f(r)$ and the fixed point iteration could be performed.
In the following example, we'll find the root of $f(x) = x - \cos(x)$ by iterating over the funcion $g(x) = \cos(x)$.
End of explanation
"""


# Consider this funtion
g = lambda x: -(3/2)*x**2+(11/2)*x-2
# Here we compute the derivative of it.
gp = lambda x: -3*x+11/2

# We plot now the funcion itself (red), its derivative (magenta) and the function y=x (blue).
# We also plot the values -1 and 1 with green dashed curves.
# This analyis shows that the fixed point, which is the intersection between teh red and blue curves, 
# does not generate a convergent fix-point-iteration since the derivative (magenta curve) has a value
# lower then -1 about the fized point.
x=np.linspace(2,3,100)
plt.figure(figsize=(8,8))
plt.plot(x,g(x),'r-',label=r'$g(x)$')
plt.plot(x,x,'b-',label=r'$x$')
plt.plot(x,gp(x),'m-',label=r'$g\'(x)$')
plt.plot(x,gp(x)*0+1,'g--')
plt.plot(x,gp(x)*0-1,'g--')
plt.grid(True)
plt.legend(loc='best')
plt.show()

"""
Explanation: Let's quickly explain the Cobweb Diagram we have here. The <font color="blue">blue</font> line is the function $y=x$ and the <font color="red">red</font> is the function $y=g(x)$. 
The point in which they meet is $r=g(r)$, i.e. the fixed point. 
In this particular example, we start at <font color="blue">$y = x = 1.5$</font> (the top right corner) and then we "jump" vertically to <font color="red">$y = \cos(1.5) \approx 0.07$</font>. 
After this, we jump horizontally to <font color="blue">$x = \cos(1.5) \approx 0.07$</font>. 
Then, we jump again vertically to <font color="red">$y = \cos\left(\cos(1.5)\right) \approx 0.997$</font> and so on. 
See the pattern here? We're just iterating over $x = \cos(x)$, getting closer to the center of the diagram where the fixed point resides, in $x \approx 0.739$.
It's very important to mention that the algorithm will converge only if the rate of convergence $S < 1$, where $S = \left| g'(r) \right|$. 
If you want to use this method, you'll have to construct $g(x)$ starting from $f(x)$ accordingly. 
In this example, $g(x) = \cos(x) \Rightarrow g'(x) = -\sin(x)$ and $|-\sin(0.739)| \approx 0.67$.
Quick question: Do you see the value 0.67 in the previous table?
Another example. Look at this web page to undertand the context: https://divisbyzero.com/2008/12/18/sharkovskys-theorem/amp/?__twitter_impression=true
End of explanation
"""


fpi(g, 1, 12, True)

# Suggestion, try the following alternative.
#fpi(g, 2.5, 100, True)

"""
Explanation: What it is interesting about the previous example is that it generates an interesting limit cicle! In the next cell we evaluate the fixed point with initial guess equal to 1. The iteration oscilates generating the following sequence: 1, 2, 3, 1, 2, 3, .... Which is nice!
End of explanation
"""


# This is a "palta" hidden in the code! Think about it. Quick question: what is it doing?
a=-1/(1-(-1.7)) # a = -1 / (1 - g'(r)), where does "a" come from?
G = lambda x: x+a*(x-g(x))
fpi(G, 1, 14, True)

"""
Explanation: However, we prefer convergent fixed-point-iterations! Here is interesting way to make a non-convergent FPI into a convergent one.
End of explanation
"""


# These are the three functions proposed.
g1 = lambda x: 1-x**3
g2 = lambda x: (1-x)**(1/3)
g3 = lambda x: (1+2*x**3)/(1+3*x**2)
# Change the input function to evaluate different functions.
# Are the three functions convergent fixed point iterations?
fpi(g3, 0.75, 20, True)

# This is a "hack" to improve the convergence of g2!
a  = -0.6
g4 = lambda x: x+a*(x-g2(x))
fpi(g4, 0.75, 10, True)
# Why does this hack works?

"""
Explanation: <div id='fpi-textbook-example' />

FPI - example from textbook
Back to TOC
This example is from the textbook. We are trying to find a root of $f(x)=x^3+x-1$.
End of explanation
"""


g1p = lambda x: -3*x**2
g2p = lambda x: -(1/3)*(1-x)**(-2/3)
g3p = lambda x: ((1+3*x**2)*(6*x**2)-(1+2*x**3)*6*x)/((1+3*x**2)**2)
g4p = lambda x: 1+a*(1-g2p(x))
r=0.6823278038280194
print('What is the conclusion then?')
print([g1p(r), g2p(r), g3p(r), g4p(r)])

# Or it may be better to apply the absolute value.
print(np.abs([g1p(r), g2p(r), g3p(r), g4p(r)]))

"""
Explanation: Now that we have found the root, let's compute the derivative of each $g(x)$ used previously and understand what exactly was going on.
End of explanation
"""


def newton_method(f, fp, x0, rel_error=1e-8, m=1, maxNumberIterations=100):
    #Initialization of hybrid error and absolute
    hybrid_error = 100
    error_i = np.inf
    
    #print('i |     x_i      |     x_{i+1}  |     |x_{i+1}-x_i|      |  e_{i+1}/e_i | e_{i+1}/e_i^2')
    #print('----------------------------------------------------------------------------------------')
    
    # Output table to store the numerical evolution of the algorithm
    output_table = []
    
    #Iteration counter
    i = 0
    while (hybrid_error > rel_error and hybrid_error < 1e12 and i<=maxNumberIterations):
        #Newton's iteration
        x1 = x0-m*f(x0)/fp(x0)
        
        #Checking if root was found
        if f(x1) == 0.0:
            hybrid_error = 0.0
            break
        
        #Computation of hybrid error
        hybrid_error = abs(x1-x0)/np.max([abs(x1),1e-12])
        
        #Computation of absolute error
        error_iminus1 = error_i
        error_i = abs(x1-x0)
        
        # Storing output data
        output_table.append([i,x0,x1,error_i,error_i/error_iminus1,error_i/(error_iminus1**((1+np.sqrt(5))/2.)),error_i/(error_iminus1**2)])
        
        #Updating solution
        x0 = x1
        #Increasing iteration counter
        i += 1
    
    # Showing final output table
    columns    = ['$i$', '$x_i$', '$x_{i+1}$', '$e_i$', r'$\frac{e_i}{e_{i-1}}$', r'$\frac{e_i}{e_{i-1}^\alpha}$', r'$\frac{e_i}{e_{i-1}^2}$']
    df = pd.DataFrame(data=output_table, columns=columns)
    display(df)
    
    #Checking if solution was obtained
    if hybrid_error < rel_error: 
        return x1
    elif i>=maxNumberIterations:
        print('Newton''s Method did not converge. Too many iterations!!')
        return None
    else:
        print('Newton''s Method did not converge!')
        return None

"""
Explanation: <div id='nm' />

Newton's Method
Back to TOC
The Newton's method also finds a root of a function $f(x)$ but it requires its derivative, i.e. $f'(x)$.
The algorithm is as follows:
\begin{align}
    x_0 &= \text{initial guess},\
    x_{i+1} &= x_i - \dfrac{f(x_i)}{f'(x_i)}.
\end{align}
For roots with multiplicity equal to 1, Newton's method convergens quadratically. However, when the multiplicity is larger that 1, it will show linear convergence. Fortunately, we can modify Newton's method if we know the multiplicity of the root, say $m$, this is as follows:
\begin{align}
    x_0 &= \text{initial guess},\
    x_{i+1} &= x_i - m\,\dfrac{f(x_i)}{f'(x_i)}.
\end{align}
This modified version will also show quadratic convergence!
End of explanation
"""


# Example funtion
f = lambda x: np.sin(x)
# The derivative of f
fp = lambda x: np.cos(x)
newton_method(f, fp, 3.1,rel_error=1e-15)

"""
Explanation: First example, let's compute a root of $\sin(x)$, near $x_0=3.1$.
End of explanation
"""


f = lambda x: x**2
fp = lambda x: 2*x # the derivative of f
newton_method(f, fp, 3.1, rel_error=1e-1, m=1, maxNumberIterations=10)

"""
Explanation: Now, we will look at the example when Newton's method shows linear convergence.
End of explanation
"""


x = sym.symbols('x', reals=True)
W=1
for i in np.arange(1,21):
        W*=(x-i)
W # Printing W nicely

# Expanding the Wilkinson polynomial
We=sym.expand(W)
We 

# Just computiong the derivative
Wep=sym.diff(We,x)
Wep

# Lamdifying the polynomial to be used with sympy
P=sym.lambdify(x,We)
Pp=sym.lambdify(x,Wep)

# Using scipy function to compute a root
root = optimize.newton(P,16)
print(root)

newton_method(P, Pp, 16.01, rel_error=1e-10, maxNumberIterations=10)

"""
Explanation: So, in the previous example Newton's method showed linear convergence.
But, how can we uss its outcome to improve the convergence?
This can fixed by understanding the following facts:
1. Linear convergence definition: $e_{i+1}/e_i=S$
2. Linear convergence exhibit by Newton's method when the root has multiplicity greater than 1: $S=(m-1)/m$
Connecting the two previous two facts we get, 
$$e_{i+1}/e_i=(m-1)/m$$.
From the table we obtain that $e_{i+1}/e_i\approx 0.5$, this implies the following equation,
$$0.5=(m-1)/m.$$
Solving for $m$ we get $m=2$.
Knowing this is very useful because we can use it with the Newton's method and recover its quadratic convergence!
<div id='wilkinson' />

Wilkinson Polynomial
Back to TOC
https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial
Final question: Why is the root far far away from $16$?
End of explanation
"""


print('Please try to solve it before you see the solution!!!')

"""
Explanation: <div id='acknowledgements' />

Acknowledgements
Back to TOC

Material created by professor Claudio Torres (ctorres@inf.utfsm.cl) and assistants: Laura Bermeo, Alvaro Salinas, Axel Simonsen and Martín Villanueva. DI UTFSM. March 2016. v1.1.
Update April 2020 - v1.32 - C.Torres : Re-ordering the notebook.
Update April 2021 - v1.33 - C.Torres : Updating format and re-re-ordering the notebook. Adding 'maxNumberIterations' to bisection, fpi and Newton's method. Adding more explanations.
Update April 2021 - v1.33 - C.Torres : Updating description and solution of 'Proposed classwork'.
Update September 2021 - v1.35 - C.Torres : Updating and commeting code more.
Update September 2021 - v1.36 - C.Torres : Updating the way we show the output tables.
Update September 2021 - v1.37 - C.Torres : Fixing typo suggested by Nicolás Tapia 2021-2. Thanks Nicolás! And removing extra code.

<div id='extraexamples' />

Extra examples
Back to TOC
Propose Classwork

Build a FPI such that given $a$ computes $\displaystyle \frac{1}{a}$. The constraint is that you can't use a division in the 'final' FPI. Write down your solution below or go and see the solution. 

Hint: I strongly suggest to use Newton's method.
End of explanation
"""


g1 = lambda x: (4./5.)*x+1./x
g2 = lambda x: x/2.+5./(2*x)
g3 = lambda x: (x+5.)/(x+1)
fpi(g1, 3.0, 10, True)

"""
Explanation: Build an algorithm that computes $\log(x_i)$ for $x_i=0.1i+0.5$, for $i\in{0,1,2,\dots,10}$. The only special function available is $\exp(x)$, in particular use np.exp(x). You can also use $$, $÷$, $+$, and $-$. It would be nice to use the result from previous example to replace $÷$.

In class
Which function shows quadratic convergence? Why?
End of explanation
"""


# What is 'a'? Can we find another 'a'?
a = -3*(1.7**2)
print(a)

f = lambda x: x**3-7
g =  lambda x: f(x)/a+x
r=fpi(g, 1.7, 14, True)
print(f(r))

"""
Explanation: Building a FPI to compute the cubic root of 7
End of explanation
"""


f = lambda x: 8*x**4-12*x**3+6*x**2-x
fp = lambda x: 32*x**3-36*x**2+12*x-1

x = np.linspace(-1,1,10000)
plt.figure(figsize=(10,10))

plt.title('What are we seeing with the semiloigy plot? Is this function differentiable?')
plt.semilogy(x,np.abs(f(x)),'b-',label=r'$|f(x)|$')
plt.semilogy(x,np.abs(fp(x)),'m-',label=r'$|fp(x)|$')
plt.grid()
plt.legend()
plt.xlabel(r'$x$',fontsize=16)
plt.show()

r=newton_method(f, fp, 0.4, rel_error=1e-8, m=1)
print([r,f(r)])
# Is this showing quadratic convergence? If not, can you fix it?

"""
Explanation: Playing with some roots
The following example proposed a particular function $f(x)$. 
The idea here is first obtain an initial guess for applying the Newton's method from the plot in semilogy scale.
The plot of $f(x)$ (blue) shows that there seems to be 2 roots in the interval plotted.
Now, the plot of $f'(x)$ (magenta) indicates that the derivative may also have a 0 together with a root, this means that the multiplicity of that root may be higher than 1. Do you see it?
End of explanation
"""


# We are finding the 1/a
# Solution code:
a = 2.1
g = lambda x: 2*x-a*x**2
gp = lambda x: 2-2*a*x
r=fpi(g, 0.7, 7, flag_cobweb=False)
print('Reciprocal found :',r)
print('Reciprocal computed explicitly: ', 1/a)
# Are we seeing quadratic convergence?

"""
Explanation: <div id='sol1' />

Solutions
Problem: Build a FPI such that given $a$ computes $\displaystyle \frac{1}{a}$
Solution: Apply Newton's method to the following function: $f(x)=a-\frac{1}{x}$. Make sure you simplify the final expression.
End of explanation
"""


xx=np.linspace(0.2,0.8,1000)
plt.figure(figsize=(10,10))
plt.plot(xx,g(xx),'r-',label=r'$g(x)$')
plt.plot(xx,gp(xx),'m-',label=r'$gp(x)$')
plt.plot(xx,xx,'b-',label=r'$x$')
plt.plot(xx,0*xx+1,'k--')
plt.plot(xx,0*xx-1,'k--')
plt.legend(loc='best')
plt.grid()
plt.show()

"""
Explanation: What is this plot telling us?
This plots shows that, even if we don't know the exact value of $g'(r)$, we can determine if the FPI will convergan by looking at the plot.
In this plot we observe that when plotting $g'(x)$ (magenta), we can determine that the value of $|g'(r)|$ will be less than 1 since it is between the black lines, that are located at $y=-1$ and $y=1$.
End of explanation
"""


