%matplotlib inline
import numpy
import matplotlib.pyplot as plt

"""
Explanation: <table>
 <tr align=left><td><img align=left src="../images/CC-BY.png">
 <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Marc Spiegelman,  Based on ipython notebook by Kyle Mandli from his course [Introduction to numerical methods](https://github.com/mandli/intro-numerical-methods)</td>
</table>
End of explanation
"""


f = lambda x: x - numpy.cos(x)

"""
Explanation: Root Finding
GOAL: Find where $f(x) = 0$.
Example:  A simple non-linear problem
let $f(x) = x - \cos(x)$,  find  values of $x$ where $f(x)=0$.
Because $f$ is non-linear it is possible that it has no roots, a finite number of roots or even an infinite number of roots.  Thus the first thing one should do is try to visualize $f$ over a range of $x$ to see, qualitatively, whether there are any zero crossings and to identify brackets where $f$ changes sign.
Here we will introduce our function using an "inlined-function" or "lambda function" in python
End of explanation
"""


print 'f(0.)=',f(0.)

"""
Explanation: which simply replaces any use of f(x) with x - cos(x), for example
End of explanation
"""


x = numpy.linspace(-10.,10., 100)

plt.figure()
plt.plot(x,f(x),'b')
plt.hold(True)
plt.plot(x,numpy.zeros(x.shape),'r--')
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("A plot")
plt.show()

"""
Explanation: or let's plot this for $x\in [-10,10]$
End of explanation
"""


from scipy.optimize import brentq

# give a bracket [a,b] such that f(a)*f(b) <= 0
a = 0.
b = 5.
x0 = brentq(f,a,b)
print
print "root x0 = {0}, in bracket [{1},{2}]".format(x0,a,b)
print "residual f(x0) = {0}".format(f(x0))

plt.figure()
plt.plot(x,f(x),'b')
plt.hold(True)
plt.plot(x,numpy.zeros(x.shape),'r--')
plt.plot(x0,f(x0),'go')
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("A root at $x_0={0}$".format(x0))
plt.show()

"""
Explanation: Which, in this range, has a single root somewhere in the bracket $0<x<5$.  The question is how to find it?
The "right answer"
The correct approach is to use (and understand) a good algorithm for bracketed root finding of functions of a single variable of which scipy.optimize provides several.  Here we will use the brentq algorithm which is a workhorse for rootfinding as it is guaranteed to find at least one root given a proper bracket (where $f$ changes sign).  The method is a generalized secant method so doesn't require derivatives of $f$ and has super-linear convergence (a simple bisection scheme also guarantees a root but has only linear convergence).
End of explanation
"""


g = lambda x: numpy.cos(x) 

xn = numpy.zeros(21)
for i in xrange(len(xn)-1):
    print "step {0}: x = {1}, residual f(x) = {2}".format(i,xn[i], f(xn[i]))
    xn[i+1] = g(xn[i])

"""
Explanation: Successive Substition (Fixed point iteration)
A more naive approach is to rewrite $f(x) = 0$ as a fixed point iteration
$$ x = g(x) $$
where $g(x)$ is another function such that when $x$ satisfies this equation,  it is a root of $f$.  For example here we could choose $g(x) = \cos(x)$. 
We can turn this equation into an iterative method by setting $x_0 =0.$ and forming a sequence of numbers
$$ x_n = g(x_{n-1})$$
and hope it converges.  Algorithmically we could do something like
End of explanation
"""


plt.figure()
plt.plot(range(len(xn)),f(xn),'b-o')
plt.xlabel('Iterations')
plt.ylabel('Residual $f(x)$')
plt.title('Convergence of fixed point iteration)')
plt.show()

"""
Explanation: and plot out the residual.
End of explanation
"""


plt.figure()
plt.semilogy(range(len(xn)),numpy.abs(f(xn)),'b-o')
plt.xlabel('Iterations')
plt.ylabel('Residual $|f(x)|$')
plt.title('Convergence of fixed point iteration')
plt.show()

"""
Explanation: which oscillates towards the value of the root.  If we look at the absolute value of the error, we see that it converges linearly i.e.
$$ |e_{n+1}| = K|e_{n}| $$
Where $K$ is a value that analysis shows should be $K=|g'(x^)|$ where $x^$ is the root.  For our problem $g'(x) = -\sin(x)$ and $K=0.673612029183$.  Because $K<1$ the fixed point iteration is a "contraction" and the error eventually $\rightarrow 0$ as $n\rightarrow\infty$.  We demonstrated that this works for this problem graphically
End of explanation
"""


print
for i in range(len(xn)-1):
    print 'Step = {0}, K={1}'.format((i+1),numpy.abs(f(xn[i+1]))/numpy.abs(f(xn[i])))
    
gprime = lambda x: -numpy.sin(x)
print
print "|g'(x0)| = {0}".format(numpy.abs(gprime(x0)))

"""
Explanation: and numerically by comparing the ratio of $|f(x_{n+1})|/|f(x_n)|$
End of explanation
"""


fprime = lambda x: 1. + numpy.sin(x)

xnn = numpy.zeros(10)
print "\nNewton's Method\n"
i = 0
tol = 1.e-16
while numpy.abs(f(xnn[i])) > tol:
    print "step {0}: x = {1}, residual f(x) = {2}".format(i,xnn[i], f(xnn[i]))
    xnn[i+1] = xnn[i] - f(xnn[i])/fprime(xnn[i])
    i += 1
    
imax = i
xnn = xnn[:imax]

"""
Explanation: Newton's Method
A  potentially more efficient method for Non-linear problems is Newton's method which can be considered another fixed-point iteration but promises much better convergence (near the fixed point, for simple roots).  
The basic idea is that given some initial guess $x_n$ such that $f(x_n) \neq 0$, there is some correction $\delta_n$ such that $f(x_n + \delta_n) = 0$.  Expanding $f$ in a Taylor series around $x_n$ we get the linear approximation
$$ f(x_n + \delta_n) \approx f(x_n) + f'(x_n)\delta_n + O(\delta_n^2) = 0$$
neglecting terms of order $\delta_n^2$ we can solve for the correction that would be exact if the problem were linear, i.e.
$$ \delta_n = -f(x_n)/f'(x_n) $$
then the next iterate is given by 
$$ x_{n+1} = x_{n} + \delta_n $$
and iterate until the residual $|f(x)| < \mathrm{tol}$ for some tolerance.
Algorithmically...
End of explanation
"""


plt.figure()
plt.semilogy(range(len(xn)),numpy.abs(f(xn)),'b-o',label='fixed point')
plt.hold(True)
plt.semilogy(range(len(xnn)),numpy.abs(f(xnn)),'r-o',label='newton')
plt.xlabel('Iterations')
plt.ylabel('Residual $|f(x)|$')
plt.legend(loc='best')
plt.title('Comparison of Fixed point iteration to Newtons Method')
plt.show()

"""
Explanation: Analysis shows that near a simple root,  Newton's method converges quadratically i.e.
$$|e_{n+1}| = C |e_n|^2$$
thus doubling the number of significant digits per iteration. This analysis is only valid near the root and in general, Newton's method can be highly unstable (for example if it finds a region where $f'(x)$ is close to zero), and in general requires some additional controls to maintain a bracket.
Comparing the two methods for this problem, however, shows that Newton's method converges quadratically,  while the fixed point iteration converges linearly
End of explanation
"""


