MTH 337: Week 7¶
Error Handling¶
A try-except statement can be used to handle errors more deliberately than having the code break. The syntax is:
try:
<code to try to execute>
except:
<code to execute if an error is generated>n = 1
while True:
try:
x = 2.**n
n += 1
except:
print n - 1
break
An error type can also be specified. In this case, different code blocks can be specified for handling different types of error.
try:
<code to try to execute>
except ErrorType:
<code to execute if an error is generated>n = 1
while True:
try:
x = 2.**n
n += 1
except OverflowError:
print n - 1
break
Animation¶
We will be using the matplotlib "animation" module.
To use this with IPython notebook, the "%pylab" magic needs to be used rather than "%pylab inline".
%pylab
from matplotlib import animation
Planetary orbit simulation¶
The planetary state is a 4-element array [x y vx vy].
The function "diff_P" returns the time derivative of the state of the planet.
def diff_P(P):
G = 1.0
M = 1.0
x, y, vx, vy = P
r3 = (sqrt(x**2 + y**2))**3
ax = -G*M*x/r3
ay = -G*M*y/r3
return array([vx, vy, ax, ay])
Define a function to return a 2d array "solving" the differential equations of planetary motion using Heun's Method.
The rows contain the state of the planet for successive time steps.
def calculate_orbit(steps, h, x=3., y=0., vx=0, vy=.5):
orbit = empty((steps, 4))
P = array([x, y, vx, vy])
for i in xrange(steps):
orbit[i] = P
Ptilde = P + h*diff_P(P)
P += h*(diff_P(P) + diff_P(Ptilde))/2.
return orbit
Creating an Animation¶
"FuncAnimation" is used to create an animation. We need to define:
An initialization function to be called once at the start.
An animation function to be called once for each frame.
Note that these functions should return a "tuple" of the objects which need to be redrawn.
steps = 5000
h = .1
fig, ax = subplots(figsize=(6, 6))
ax.set_aspect('equal')
# "point" repesents the planet
point, = plot([], [], 'ro')
sun = plot(0, 0, 'yo', ms=8)
# The size of the displayed area is fixed at the start, to avoid resizing during the animation
size = 4
xlim(-size, size)
ylim(-size, size)
# Calculate the complete orbit in advance.
orbit = calculate_orbit(steps, h)
# Initialization function for the start of the animation
def init():
point.set_data([], [])
return point,
# Function to update the data for each frame of the animation
def animate(i):
point.set_data(orbit[i, 0], orbit[i, 1])
return point,
# This is called to generate and run the animation
animation.FuncAnimation(fig, animate, init_func=init, frames=steps, interval=20, blit=True, repeat=True)
Plotting the orbit as well as the planet¶
A line needs to be drawn for the orbit, as well as a point for the planet.
steps = 5000
h = .1
fig, ax = subplots(figsize=(6, 6))
ax.set_aspect('equal')
ax.set_axis_bgcolor('k')
point, = plot([], [], 'ro')
line, = plot([], [], 'cyan')
sun = plot(0, 0, 'yo', ms=8)
size = 4
xlim(-size, size)
ylim(-size, size)
orbit = calculate_orbit(steps, h)
def init():
point.set_data([], [])
line.set_data([], [])
return line, point,
def animate(i):
point.set_data(orbit[i, 0], orbit[i, 1])
line.set_data(orbit[:i, 0], orbit[:i, 1])
return line, point,
animation.FuncAnimation(fig, animate, init_func=init, frames=steps, interval=20, blit=True, repeat=True)
Using the "subplot'' function¶
The function "subplot" can be used to divide a figure into a grid of rows and columns. Individual cells in the grid can be selected for plotting using
subplot(nrows, ncols, i)
where
- nrows is the number of rows in the grid
- ncols is the number of columns in the grid
- i is the cell index. This starts at 1 in the top-left corner of the grid, and increases moving left to right and top to bottom.
After a cell is selected using subplot, all plotting is done in that cell until another cell is selected.
%pylab inline
Example: Lissajous Figures¶
end = 3
t = linspace(0, 2*pi, 100)
i = 1
figure(figsize=(6,6))
for row in range(1, end + 1):
for col in range(1, end + 1):
subplot(end, end, i)
plot(cos(row*t), sin(col*t))
i += 1
Axis ticks can be turned off if they are getting too crowded¶
end = 5
t = linspace(0, 2*pi, 100)
i = 1
figure(figsize=(8, 8))
for row in range(1, end + 1):
for col in range(1, end + 1):
subplot(end, end, i)
plot(cos(row*t), sin(col*t))
xticks([])
yticks([])
if i % 5 == 1:
ylabel(str(row))
if i > 20:
xlabel(str(col))
i += 1
# A title can be added to the figure using the "suptitle" function
suptitle("Lissajous Figures")
