This lecture by Aron Ahmadia and David Ketcheson is licensed under a Creative Commons Attribution 4.0 International License. All code examples are also licensed under the MIT license.
from IPython.core.display import HTML
css_file = './example.css'
HTML(open(css_file, "r").read())
One of the major advantages of the notebook over a traditional textbook is that plots of numerical solutions can be included with the code that produces them -- and the student can modify and execute that code, to see how the result changes! Better yet, animations and interactive widgets can be included to show the evolution of time-dependent solutions, or the dependence of results on some parameter(s). Here we'll explore some of the tools for producing such plots and animations. Some of these are very new and are still evolving rapidly.
In order to make things more interesting and relevant, we'll demonstrate these options in the context of a PDE solution.
We'll solve a system of reaction-diffusion PDEs in two dimensions:
\begin{align} u_t & = \delta D_1 \nabla^2 u + f(u,v) \\ v_t & = \delta D_2 \nabla^2 v + g(u,v) \end{align}where $\nabla^2 u = u_{xx} + u_{yy}$ denotes the Laplacian and $f,g$ represent reaction terms.
For simplicity, we'll consider the square domain $[-1,1]\times[-1,1]$ with periodic boundary conditions; i.e., the conditions on $u(x,y,t)$ are
\begin{align} u(-1,y,t) & = u(1,y,t) \\ u(x,-1,t) & = u(x,1,t) \end{align}with corresponding conditions on $v$. The reaction terms we will use are
\begin{align} f(u,v) & = \alpha u (1-\tau_1 v^2) + v(1-\tau_2 u) \\ g(u,v) & = \beta v + \alpha \tau_1 u v^2 + u (\gamma + \tau_2 v). \end{align}import numpy as np
import scipy.optimize
import scipy.sparse
def f(u,v):
return alpha*u*(1-tau1*v**2) + v*(1-tau2*u);
def g(u,v):
return beta*v*(1+alpha*tau1/beta*u*v) + u*(gamma+tau2*v);
def five_pt_laplacian_sparse_periodic(m,a,b):
"""Construct a sparse matrix that applies the 5-point laplacian discretization
with periodic BCs on all sides."""
e=np.ones(m**2)
e2=([1]*(m-1)+[0])*m
e3=([0]+[1]*(m-1))*m
h=(b-a)/(m+1)
A=scipy.sparse.spdiags([-4*e,e2,e3,e,e],[0,-1,1,-m,m],m**2,m**2)
# Top & bottom BCs:
A_periodic = scipy.sparse.spdiags([e,e],[m-m**2,m**2-m],m**2,m**2).tolil()
# Left & right BCs:
for i in range(m):
A_periodic[i*m,(i+1)*m-1] = 1.
A_periodic[(i+1)*m-1,i*m] = 1.
A = A + A_periodic
A/=h**2
A = A.todia()
return A
def one_step(u,v,k,A,delta,D1=0.5,D2=1.0):
u_new = u + k * (delta*D1*A*u + f(u,v))
v_new = v + k * (delta*D2*A*v + g(u,v))
return u_new, v_new
def step_size(h,delta):
return h**2/(5.*delta)
As a first example of a plot, let's just look at the sparsity structure of the numerical laplacian matrix produced by the provided function above. By default, matplotlib plots will open in a separate window, just as they would if we were plotting from an IPython command line. To get them to appear in the notebook, we use an IPython magic function
%matplotlib inline
A = five_pt_laplacian_sparse_periodic(4,-1.,1.)
import matplotlib.pyplot as plt
plt.spy(A)
Now let's solve the reaction-diffusion PDE and plot the solution (note that the methods we're using here are not very accurate or efficient, but they're good enough to give a qualitatively correct solution in reasonable time on a small grid).
delta=0.0021; tau1=3.5; tau2=0; alpha=0.899; beta=-0.91; gamma=-alpha;
def set_up(m, T, a=-1., b=1.):
# Set up the grid
a=-1.; b=1.
h=(b-a)/m; # Grid spacing
x = np.linspace(a,b,m) # Coordinates
y = np.linspace(a,b,m)
# Initial data
u=np.random.randn(m,m)/2.;
v=np.random.randn(m,m)/2.;
plt.clf(); plt.hold(False)
plt.pcolormesh(x,y,u); plt.colorbar(); plt.axis('image');
plt.draw()
u=u.reshape(-1)
v=v.reshape(-1)
A=five_pt_laplacian_sparse_periodic(m,-1.,1.)
k = step_size(h,delta) # Time step size
N = int(round(T/k)) # Number of steps to take
return x, y, u, v, A, k, N
def pattern_formation(m,T):
r"""Model pattern formation by solving a reaction-diffusion PDE on a periodic
square domain with an m x m grid."""
x, y, u, v, A, k, N = set_up(m,T)
t=0. # Initial time
#Now step forward in time
for j in range(N):
u,v = one_step(u,v,k,A,delta)
t = t+k;
# Plot the final solution
U=u.reshape((m,m))
plt.pcolormesh(x,y,U)
plt.colorbar()
plt.axis('image')
plt.show()
pattern_formation(m=100,T=200)
Having the plot in the notebook is nice, but there are major drawbacks:
There are straightforward ways to get multiple plots from one cell, and even to plot several snapshots of a time-dependent solution. But it would be much nicer if the solution plot actually evolved in time. Here's a crude way to accomplish that:
import time
from IPython.display import display, clear_output
def pattern_formation(m,T):
x, y, u, v, A, k, N = set_up(m,T)
t=0. # Initial time
#Now step forward in time
next_plot = 0
for j in range(N):
u,v = one_step(u,v,k,A,delta)
t = t+k;
#Plot every t=5 units
if t>next_plot:
clear_output(wait=True)
next_plot = next_plot + 5
U=u.reshape((m,m))
time.sleep(0.2)
plt.pcolormesh(x,y,U); plt.axis('image'); plt.title(str(t))
fig=plt.gcf(); display(fig)
return U
U = pattern_formation(m=100,T=50)
plt.close('all')
What if we want to be able to zoom and pan? That's a bit harder, but there are multiple solutions either just developed or in the works.
One is Plotly, a web-based service that generates interactive plots and allows them to be embedded in the notebook. To use it, you should set up an account. But for this tutorial, you can just run the following code using my account.
import plotly
import plotly.plotly as py
py.sign_in('DavidKetcheson','mgs2lgb203')
py.iplot([plotly.graph_objs.Heatmap(z=U)],width=500,height=500)
Notice that you can get actual data values by hovering over the plot. You can also adjust the plot's look interactively on the Plotly website. For many more examples of using Plotly in the IPython notebook, see this notebook and this notebook, both of which are introductions to Plotly. For a few more examples of plotting numerical simulation results with Plotly, look at my notebook on Stegotons.
Importantly, usage of Plotly requires an internet connection.
Another package that can produce interactive plots in the notebook (this time without needing an internet connection) is Bokeh. I couldn't get a pcolor
plot example working with Bokeh (I think it's possible, but it's in beta and changing rapidly). Here is a very simple line plot example.
from bokeh.plotting import output_notebook, scatter, show, line
output_notebook()
x = np.linspace(0, 4*np.pi, 100)
y = np.sin(x)
line(x,y, color="#FF00FF", tools="pan,wheel_zoom,box_zoom,reset,resize")
show()
Animated plots are also possible with Bokeh, but we won't spend more time on it because the interface will change in the near future (they plan to integrate with IPython widgets).
Some other projects that allow (or may soon allow) for interactive plots in the browser include:
The last two are still too new to go into in this tutorial, but may be useful soon.
Thanks to Jake Vanderplas' JSAnimation library, we have a much better way to include animations of time-dependent solutions. I've been using JSAnimation heavily in my teaching notebooks, but it will soon be replaced by tools based on IPython widgets.
If you're running on SageMathCloud, you already have access to JSAnimation through Clawpack, via
from clawpack.visclaw.JSAnimation import IPython_display
If you're working locally, the easiest way to get it is
git clone https://github.com/jakevdp/JSAnimation.git
cd JSAnimation
python setup.py install
after which you can
from JSAnimation import IPython_display
def pattern_formation(m=10,T=1000):
# Set up the grid
a=-1.; b=1.
h=(b-a)/m; # Grid spacing
x = np.linspace(a,b,m) # Coordinates
y = np.linspace(a,b,m)
Y,X = np.meshgrid(y,x)
# Initial data
u=np.random.randn(m,m)/2.;
v=np.random.randn(m,m)/2.;
frames = [u]
u=u.reshape(-1)
v=v.reshape(-1)
A=five_pt_laplacian_sparse_periodic(m,-1.,1.)
t=0. # Initial time
k = step_size(h,delta) # Time step size
N = int(round(T/k)) # Number of steps to take
#Now step forward in time
next_plot = 0
for j in range(N):
#Plot every t=5 units
if t>=next_plot:
next_plot = next_plot + 5
U=u.reshape((m,m))
frames.append(U)
u,v = one_step(u,v,k,A,delta)
t = t+k;
return x,y,frames
x, y, frames = pattern_formation(m=100,T=200)
from matplotlib import animation
import matplotlib.pyplot as plt
#from clawpack.visclaw.JSAnimation import IPython_display # Works on SMC
from JSAnimation import IPython_display
import numpy as np
fig = plt.figure(figsize=[4,4])
U = frames[0]
# This essentially does a pcolor plot, but it returns the appropriate object
# for use in animation. See http://matplotlib.org/examples/pylab_examples/pcolor_demo.html.
# Note that it's necessary to transpose the data array because of the way imshow works.
plot_handle = plt.imshow(U.T, vmin=U.min(), vmax=U.max(),
extent=[x.min(), x.max(), y.min(), y.max()],
interpolation='nearest', origin='lower')
def fplot(frame_number):
U = frames[frame_number]
plot_handle.set_data(U.T)
return plot_handle,
animation.FuncAnimation(fig, fplot, frames=len(frames), interval=20)
Notice that FuncAnimation
takes three arguments:
Using JSAnimation can be a little tricky, because we need to provide a function (called fplot
here) that returns a handle to a plot of one frame of the solution. It's not enough for fplot
to simply plot the solution. That's why we first get a handle to the plot and then use set_data
to modify the plot each time fplot
is called.
Also, notice that the frames
argument can be any list. Typically, it would be a list of numbers, but more generally it is just the list of argument values that need to be passed to fplot
.
See also: IPython widgets (up next)