| 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) |
|
In [ ]:
%matplotlib inline
import numpy
import matplotlib.pyplot as plt
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
In [ ]:
f = lambda x: x - numpy.cos(x)
which simply replaces any use of f(x) with x - cos(x), for example
In [ ]:
print 'f(0.)=',f(0.)
or let's plot this for $x\in [-10,10]$
In [ ]:
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()
Which, in this range, has a single root somewhere in the bracket $0<x<5$. The question is how to find it?
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).
In [ ]:
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()
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
In [ ]:
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])
and plot out the residual.
In [ ]:
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()
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
In [ ]:
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()
and numerically by comparing the ratio of $|f(x_{n+1})|/|f(x_n)|$
In [ ]:
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)))
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...
In [ ]:
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]
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
In [ ]:
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()
Neither Fixed Point iteration or Newton's method are guaranteed to converge even given a bracket and more sophisticated methods are needed for convergence and robustness (and thus the recommendation of a proper root finder like brentq from scipy). For further analysis of the convergence of fixed-point iterations and Newton's iteration as well as a host of other root finding algorithm's see Kyle Mandli's notes Introduction to numerical methods, in particular lecture 05_root_finding_optimization.ipynb
In [ ]: