Some simple computation

The main purpose of this notebook is to illustrate a bit of relatively simple Python code - the kind of code that I expect folks in my numerical analysis classes to be able to read, understand, and emulate. There is some mathematics behind the code that we'll deal with later.

Here's the basic problem: Solve the equation $\cos(x)=x$. Like a lot of interesting problems, this is very easy to state and understand, but quite hard to solve. It's impossible to solve symbolically, in fact! As we'll learn in a few weeks, though, there are fabulous algorithms to find numerical approximations to solutions. One of these, called Brent's algorithm, searches for a single root of a function in an interval. This is implemented in Scipy's optimize module as brentq. Here's how to invoke it to find a solution to our equation in the interval $[0,1]$.

In [1]:
# import the Numpy library, where the cosine is defined.
import numpy as np

# import brentq
from scipy.optimize import brentq

# define the function whose root is the solution we seek
def f(x): return np.cos(x)-x

# invoke brentq
brentq(f,0,1)
Out[1]:
0.7390851332151559

The comments and spaces are not necessary, so that's really just four lines of code.

Like many algorithms, Brent's algorithm is iterative. As it turns out in fact, we can solve the equation by simply iterating the cosine function! We'll learn a bit later why this works but here's how to implement the idea in Python.

In [2]:
x1 = 1
x2 = np.cos(1)
cnt = 0
while np.abs(x1-x2) > 10**(-8) and cnt < 100:
    x1 = x2
    x2 = np.cos(x2)
    cnt = cnt + 1
(x2,cnt)
Out[2]:
(0.73908513664657183, 45)

I guess this says that, after 45 iterates, the difference between $x$ and $\cos(x)$ is less than $1/10^8$. So, we're pretty close.

The code so far illustrates that there's a couple of libraries that we'll use a lot - NumPy (which implements a lot of relatively low level functionality) and SciPy (which implements higher level functionality built on top of NumPy).

Very often, it's nice to visualize your results. There are a number of very nice Python libraries for visualization but one of the most popular and widely used is called Matplotlib. In addition, Matplotlib is included in Anaconda, so that's what we'll use for the most part.

My expectation is that you can at least do a little basic plotting with Matplotlib. Here's how to $f(x)=\cos(x)$ together with the line $y=x$, as well as the point where they are equal.

In [3]:
# A cell magic to run Matplotlib inside the notebook.
# Not totally necessary, but sometimes nice.
%matplotlib inline

# Import the plotting module
import matplotlib.pyplot as plt

# Set up the x-values we wish to plot.
# 100 of them evenly distributed over [-0.2,1.2]
xs = np.linspace(-0.2,1.2,100)

# Plot the cosine function over those x-values
plt.plot(xs,np.cos(xs))

# Plot the line y=x
plt.plot(xs,xs)

# Plot the soution as a red dot.
plt.plot(x2,x2,'ro')
Out[3]:
[<matplotlib.lines.Line2D at 0x108f0ff60>]

If we keep track of our orbit as we progress, we can go a bit farther and illustrate the convergence of the orbit with a cobweb plot.

In [4]:
x0 = .1
x1 = x0
orbitx = [x1]
orbity = [x1]

x2 = np.cos(x1)
cnt = 0
while np.abs(x1-x2) > 10**(-8) and cnt < 100:
    orbitx.extend([x1,x2])
    orbity.extend([x2,x2])
    x1 = x2
    x2 = np.cos(x2)
    cnt = cnt + 1
xs = np.linspace(-0.2,1.2,100)
plt.plot(orbitx,orbity)
plt.plot(xs,np.cos(xs), 'k', linewidth=2)
plt.plot(xs,xs, 'k')
plt.plot(x0,x0,'go')
plt.plot(x2,x2,'ro')
Out[4]:
[<matplotlib.lines.Line2D at 0x109028f60>]

In my own images, I will often go a bit overboard with options to specify things exactly the way I want. I emphasize, this level of graphic detail is not your responsibility!

In [5]:
import matplotlib as mpl
plt.plot(orbitx,orbity)
plt.plot(xs,np.cos(xs), 'k', linewidth=2)
plt.plot(xs,xs, 'k')
plt.plot(x0,x0,'go')
plt.plot(x2,x2,'ro')

ax = plt.gca()
ax.set_aspect(1)
ax.set_ylim(-0.1,1.1)
ax.set_xlim(-0.1,1.1)
ax.set_xticks([1/2,1])
ax.set_yticks([1/2,1])

ax.spines['left'].set_position('zero')
ax.spines['right'].set_color('none')
ax.spines['bottom'].set_position('zero')
ax.spines['top'].set_color('none')

xaxis,yaxis = ax.findobj(mpl.axis.Axis)
xticks = xaxis.get_major_ticks()
for tick in xticks:
    tick.get_children()[1].set_color('w')
yticks = yaxis.get_major_ticks()
for tick in yticks:
    tick.get_children()[1].set_color('w')
In [ ]: