MTH 337: Week 9

Week 8 Notebook

Week 10 Notebook

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

meshgrid

The NumPy function "meshgrid" can be used to:

  • Apply a function of two variables to every point in a grid in the xy-plane.
  • Create a grid of points in the complex plane.

A function f(x, y) of two scalars

In [2]:
def distance(x, y):
    return sqrt(x**2 + y**2)

distance(3, 4)
Out[2]:
5.0

"meshgrid" creates two 2D arrays from two 1D arrays.

  • One array contains just the x-coordinates of points in the xy-plane
  • The other array contains just the y-coordinates of points in the xy-plane
In [3]:
x = arange(4)
y = arange(3)
print x
print y
X, Y = meshgrid(x, y)
print X
print Y

[0 1 2 3]
[0 1 2]
[[0 1 2 3]
 [0 1 2 3]
 [0 1 2 3]]
[[0 0 0 0]
 [1 1 1 1]
 [2 2 2 2]]

Functions defined on two scalars can then be applied without modification to two arrays

In [4]:
print distance(X, Y)

[[ 0.          1.          2.          3.        ]
 [ 1.          1.41421356  2.23606798  3.16227766]
 [ 2.          2.23606798  2.82842712  3.60555128]]

Using meshgrid to create a grid of points in the complex plane

First create the real and imaginary axes

In [5]:
x = linspace(-1, 1, 5)
y = linspace(-1j, 1j, 5)
print x
print y

[-1.  -0.5  0.   0.5  1. ]
[ 0.-1.j   0.-0.5j  0.+0.j   0.+0.5j  0.+1.j ]

Then use meshgrid to create 2d arrays of the real and imaginary parts of each point

In [6]:
X, Y = meshgrid(x, y)
print X
print Y

[[-1.  -0.5  0.   0.5  1. ]
 [-1.  -0.5  0.   0.5  1. ]
 [-1.  -0.5  0.   0.5  1. ]
 [-1.  -0.5  0.   0.5  1. ]
 [-1.  -0.5  0.   0.5  1. ]]
[[ 0.-1.j   0.-1.j   0.-1.j   0.-1.j   0.-1.j ]
 [ 0.-0.5j  0.-0.5j  0.-0.5j  0.-0.5j  0.-0.5j]
 [ 0.+0.j   0.+0.j   0.+0.j   0.+0.j   0.+0.j ]
 [ 0.+0.5j  0.+0.5j  0.+0.5j  0.+0.5j  0.+0.5j]
 [ 0.+1.j   0.+1.j   0.+1.j   0.+1.j   0.+1.j ]]

Finally, add the real and imaginary parts to create a 2D array of complex numbers

In [7]:
Z = X + Y
print Z

[[-1.0-1.j  -0.5-1.j   0.0-1.j   0.5-1.j   1.0-1.j ]
 [-1.0-0.5j -0.5-0.5j  0.0-0.5j  0.5-0.5j  1.0-0.5j]
 [-1.0+0.j  -0.5+0.j   0.0+0.j   0.5+0.j   1.0+0.j ]
 [-1.0+0.5j -0.5+0.5j  0.0+0.5j  0.5+0.5j  1.0+0.5j]
 [-1.0+1.j  -0.5+1.j   0.0+1.j   0.5+1.j   1.0+1.j ]]

Newton in the Complex Plane

First create a grid of points in the complex plane

In [8]:
npts = 100
x = linspace(-1, 1, npts)
y = linspace(-1j, 1j, npts)
X, Y = meshgrid(x, y)
Z = X + Y

Iterate using Newton's Method on the function $f(z) = z^3 - 1$

$\Rightarrow$ $z_{i + 1} = z_i - \dfrac{f(z_i)}{f'(z_i)} = z_i - \dfrac{z_i^3 - 1}{3z_i^2}$

In [9]:
niters = 20
for i in range(niters):
    Z -= (Z**3 - 1)/(3*Z**2)

Create variables for the roots of $f(z)$

In [10]:
r1 = complex(1, 0)
r2 = complex(-.5, sqrt(3)/2.)
r3 = complex(-.5, -sqrt(3)/2.)

Create a boolean array testing convergence to each root.

In [11]:
tol = 1e-6
root1s = abs(Z - r1) < tol
root2s = abs(Z - r2) < tol
root3s = abs(Z - r3) < tol

Finally, create an image array and set the red, green and blue components depending on which root each initial $z$ value has converged to.

In [12]:
img = zeros((npts, npts, 3))
img[:,:,0] = root1s
img[:,:,1] = root2s
img[:,:,2] = root3s
imshow(img[::-1])
Out[12]:
<matplotlib.image.AxesImage at 0x7fbc5c092a10>

Alternatively, all three components can be set at the same time using Boolean indexing

In [13]:
img = zeros((npts, npts, 3))
img[root1s] = [1, 1, 0]
img[root2s] = [1, 0, 1]
img[root3s] = [0, 1, 1]
imshow(img[::-1])
Out[13]:
<matplotlib.image.AxesImage at 0x7fbc57f468d0>

Random Numbers and Their Applications

Say we want to numerically estimate the mean value of a function by $\bar{f} = \dfrac{1}{n} \sum_{i=1}^n f(x_i)$

A "uniform grid" of equally spaced points may not be the best choice for $\{x_1, x_2, \ldots, x_n\}$

In [14]:
def testf(x):
    return 2 + 2*sin(50*x)
In [15]:
x = linspace(0, 1, 400)
sample = linspace(0, 1, 9)
plot(x, testf(x))
plot(sample, testf(sample), 'ro')
Out[15]:
[<matplotlib.lines.Line2D at 0x7fbc57f7d890>]

Choosing "random" numbers on the interval can help minimize errors and avoid bias.

Properties of "Random" Numbers

For these kinds of applications, we want such numbers to have the properties of being:

1) Uniformly distributed. Statistically, the numbers are evenly spread over the entire interval.

2) Independent of each other. There should be a lack of "pairwise" and higher-order correlation between numbers in the sequence.

Does the "Mayfly Model" generate random numbers with b = 4?

We can test for these properties experimentally by:

1) Plotting a "histogram" of the sequence of numbers - "A graphical representation of the distribution of numerical data".

2) Plotting $x_{i + 1}$ against $x_i$ for a large nunber of points. If there is no pairwise correlation, then the points will be scattered throughout the square.

The histogram shows a non-uniform distribution, with distinctly higher densities of points at the ends of the interval [0, 1].

In [16]:
x = .4
b = 4.
npts = 100000
results = empty(npts)
for i in xrange(npts):
    results[i] = x
    x = b*(1-x)*x

h = hist(results, bins=50)

Even worse, the "Mayfly Model" fails the "no pairwise correlation" test badly, since every pair ($x_i$, $x_{i+1}$) lies on the graph of $g(x) = 4(1 - x)x$.

In [17]:
npts = 1000
results = empty(npts)
for i in xrange(npts):
    results[i] = x
    x = b*(1-x)*x

plot(results[:-1], results[1:], 'b.', ms=1, alpha=.7)
Out[17]:
[<matplotlib.lines.Line2D at 0x7fbc57dca6d0>]

We want a function that effectively "fills in" the square.

Consider the graph of $g(x) = ax \mod 1$

With a = 2, we have the following graph.

In [18]:
a = 2
x = linspace(0, 1, 400)
plot(x, a*x % 1, 'b.')
Out[18]:
[<matplotlib.lines.Line2D at 0x7fbc57ba67d0>]

Increasing a, the graph gets steeper and steeper.

In [19]:
a = 20
x = linspace(0, 1, 5000)
plot(x, a*x % 1, 'b.', ms=2)
Out[19]:
[<matplotlib.lines.Line2D at 0x7fbc57ad1550>]