MTH 337: Week 5

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

The NumPy function ``loadtxt'' loads data from a file into a 2D array

The data file is expected to consist of header information, followed by a series of rows containing equal numbers of data values.

In [2]:
# The balloon data has been downloaded from the web and saved as "balloon.dat"
# There are 138 rows of header information, so we skip these using the "skiprows" keyword argument
data = loadtxt("balloon.dat", skiprows=138)
In [3]:
# The balloon data is now loaded into a 2D numpy array
data
Out[3]:
array([[  0. ,   0. ,   0. , ...,  27.5,  70. ,  21.6],
       [  0. ,   1. ,   4.1, ...,  27.4,  70. ,  21.5],
       [  0. ,   2. ,   4.1, ...,  27.3,  69. ,  21.2],
       ..., 
       [ 34. ,   8. ,   4.5, ..., -30. ,  56. , -36. ],
       [ 34. ,   9. ,   4.5, ..., -30. ,  56. , -36. ],
       [ 34. ,  10. ,   4.5, ..., -30.1,  56. , -36.1]])
In [4]:
# Use slicing to select the minutes
# Select from every row, 
mins = data[:, 0]
In [5]:
# Slice the seconds
secs = data[:, 1]
In [6]:
# Slice the temperature
temp = data[:, 5]
In [7]:
print mins
print secs
print temp

[  0.   0.   0. ...,  34.  34.  34.]
[  0.   1.   2. ...,   8.   9.  10.]
[ 27.5  27.4  27.3 ..., -30.  -30.  -30.1]
In [8]:
# Use numpy array operations to combine integer minutes and seconds into a single floating-point minutes
t = mins + secs/60.

Floating Point Arithmetic - Hazards and Fixes

Consider the function:

$g(x) = \left\{ \begin{array} \log(1 + x)/x & x \neq 1 \\ 1 & x = 1 \end{array}\right.$

Coding this function directly in Python leads to the graph below.

In [9]:
end = 10**-7
x = linspace(-end, end, 641)
y = where(x==0.0, 1.0, log(1 +  x)/x)
xlim(-end, end)
plot(x, y)

/home/adam/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: invalid value encountered in divide
  app.launch_new_instance()
Out[9]:
[<matplotlib.lines.Line2D at 0x7f8213aa57d0>]

Zooming in, we find increasing large oscillations close to zero

In [10]:
end = 10**-15
x = linspace(-end, end, 641)
y = where(x==0.0, 1.0, log(1 + x)/x)
xlim(-end, end)
plot(x, y)

/home/adam/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: invalid value encountered in divide
  app.launch_new_instance()
Out[10]:
[<matplotlib.lines.Line2D at 0x7f82137515d0>]

The problem is that:

  • When small x values are added to 1, values that differ by less than machine epsilon can end up being rounded to the same number. The numerator for such number is therefore the same.
  • These small values of x can still be represented in the denominator to full precision.
  • In the range in which $\log(1 + x)$ is returning a fixed number, the numerator is fixed but the denominator is increasing.
  • Eventully the numerator catches up, and we have a "jump" in the graph.

The Taylor's Theorem approximation to $\log(1 + x)$ is

$x - \dfrac{x^2}{2} + \dfrac{x^3}{3} - R_4$, where $|R_4| < \dfrac{x^4}{4}$

which gives

$\dfrac{\log(1 + x)}{x} = 1 - \dfrac{x}{2} + \dfrac{x^2}{3} - Q$, where $|Q| < \dfrac{x^3}{4}$

This approximation is as accurate as we need if $|Q| < $ machine epsilon ($\sim 10^{-16}) \Rightarrow x < 7 \times 10^{-6}$

In [11]:
end = 10**-7
x = linspace(-end, end, 641)
y = where(x==0, 1.0, log(1+x)/x)
xlim(-end, end)
plot(x, y)
plot(x, 1 - x/2. + (x**2)/3., 'r')

/home/adam/anaconda/lib/python2.7/site-packages/IPython/kernel/__main__.py:3: RuntimeWarning: invalid value encountered in divide
  app.launch_new_instance()
Out[11]:
[<matplotlib.lines.Line2D at 0x7f8213593fd0>]
In [ ]: