Animation techniques

The IPython Notebook is here.

1. Moving frame for 2d plotting

(Warning) If you run the following code in the notebook, the interrupt key Ctrl-C doesn't work. To stop the run, you should restart kernel, which is included in the "Kernel" menu. (停止には、jupyterのメニュー"Kernel"にある"Restart"を選択)

To avoid this inconvenience, you are recommended to edit a "py" file containing the below code and run it in the ipython console. To do this, an development environment such as spyder is a good choice.

Principal idea for ploting a moving frame:

• append the latest small part of the trace to the x and y values (line 44)
• delete the oldest segment of the trace by the same length (line 45)
• redraw the plot

The width of plot is defined by "width_plot" and the width of shifting the frame per a plotting step is given by "with_ode". (lines 24,25)

%matplotlib auto
# coding: UTF-8
"""
Lorenz model visualization as an exercise on animation
    continue to plot untile the run quits.

This doesn't work "%matplotlib inline". You should use "%matplotlib auto", "qt" or etc.
"""
import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def lorentz_deriv(x, t0, sigma=10., beta=8./3, rho=28.0):
    """Compute the time-derivative of a Lorentz system."""
    return [sigma * (x[1] - x[0]), x[0] * (rho - x[2]) - x[1], x[0] * x[1] - beta * x[2]]

# Create new Figure and an Axes which fills it.
fig = plt.figure()
#ax = fig.add_axes([0,0,1,1], axes=True)

# initial setting
step = 5
width_ode = 0.05
width_plot = 10
t = np.linspace(0,width_ode,step)
x0 = [1., 1., 1.]  # starting vector

t4plot = np.linspace(0, width_plot, int(step*width_plot/width_ode))
data4plot = integrate.odeint(lorentz_deriv, x0, t4plot)  # ode calc.  for the first plot
x_plot, = plt.plot(t4plot, data4plot[:,0], 'b-')
y_plot, = plt.plot(t4plot, data4plot[:,1], 'r-')
z_plot, = plt.plot(t4plot, data4plot[:,2], 'g-')

x0 = data4plot[-1]
start_t =0

try :
    while True :
        start_t = start_t + width_ode
        # integration
        x_t = integrate.odeint(lorentz_deriv, x0,t)
        # update data for plotting
        data4plot = np.vstack((data4plot, x_t))   # append new columns
        data4plot = data4plot[step:]             # delete "step" columns at the head
        # new value for the horizontal axis
        t4plot = np.linspace(start_t, start_t+width_plot, int(step*width_plot/width_ode))
                  
        x_plot.set_xdata(t4plot)
        x_plot.set_ydata(data4plot[:,0])
        y_plot.set_xdata(t4plot)
        y_plot.set_ydata(data4plot[:,1])
        z_plot.set_xdata(t4plot)
        z_plot.set_ydata(data4plot[:,2])
        plt.xlim(min(t4plot), max(t4plot))
        x0 = x_t[-1]                                        # initial value for the next calculation of ode
        plt.pause(0.01)                                    # important for realtime plot: you should use "plt.pause" rather than "plt.draw"
except KeyboardInterrupt:                                   # doesn't work when this starts from the notebook (???)
      pass

2. animation functions of matplotlib

By using animation function of matplotlib, you can save the animationtion to mpeg4 and moving gif file.

(a) a sample from web page

https://jakevdp.github.io/blog/2013/02/16/animating-the-lorentz-system-in-3d/

%matplotlib auto
# coding: UTF-8
# copied from https://jakevdp.github.io/blog/2013/02/16/animating-the-lorentz-system-in-3d/

import numpy as np
from scipy import integrate

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.colors import cnames
from matplotlib import animation

N_trajectories = 20


def lorentz_deriv(x, t0, sigma=10., beta=8./3, rho=28.0):
    """Compute the time-derivative of a Lorentz system."""
    return [sigma * (x[1] - x[0]), x[0] * (rho - x[2]) - x[1], x[0] * x[1] - beta * x[2]]

# Choose random starting points, uniformly distributed from -15 to 15
np.random.seed(1)
x0 = -15. + 30. * np.random.random((N_trajectories, 3))

# Solve for the trajectories　（全部計算してしまう）
t = np.linspace(0, 10, 5000)
x_t = np.asarray([integrate.odeint(lorentz_deriv, x0i, t)
                  for x0i in x0]) # 初期値ごとに軌跡用の配列を生成

# Set up figure & 3D axis for animation
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], projection='3d')
ax.axis('off')

# choose a different color for each trajectory
colors = plt.cm.jet(np.linspace(0, 1, N_trajectories))

# set up lines and points
lines = sum([ax.plot([], [], [], '-', c=c)
             for c in colors], [])
pts = sum([ax.plot([], [], [], 'o', c=c)
           for c in colors], [])

# prepare the axes limits
ax.set_xlim((-25, 25))
ax.set_ylim((-35, 35))
ax.set_zlim((5, 55))

# set point-of-view: specified by (altitude degrees, azimuth degrees)
ax.view_init(30, 0)

# initialization function: plot the background of each frame
def init():
    for line, pt in zip(lines, pts):
        line.set_data([], [])
        line.set_3d_properties([])

        pt.set_data([], [])
        pt.set_3d_properties([])
    return lines + pts

# animation function.  This will be called sequentially with the frame number
def animate(i):
    # we'll step two time-steps per frame.  This leads to nice results.
    i = (2 * i) % x_t.shape[1] # shape[1]はx_tの列数　（繰り返し描く工夫）

    for line, pt, xi in zip(lines, pts, x_t):
        x, y, z = xi[:i].T  # T は transposeのこと
        line.set_data(x, y) # 1-dim arrayの場合は列ベクトルか行ベクトルかの区別がない
        line.set_3d_properties(z)

        pt.set_data(x[-1:], y[-1:])
        pt.set_3d_properties(z[-1:])

    ax.view_init(30, 0.3 * i)
    fig.canvas.draw()
    return lines + pts

# instantiate the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=500, interval=30, blit=True)

# Save as mp4. This requires mplayer or ffmpeg to be installed
# Set blit=False when you use the following "save" function
# anim.save('lorentz_attractor.mp4', fps=15, extra_args=['-vcodec', 'libx264'])

plt.show()

(b) two traces with initial values which have a tiny amount of difference

The following is an application where two trajectories which have a tiny different initial condition. You can see the difference rapidly grows.

The trajectories are updated infinitely until it's interrupted, while in the above program, trajectories are calculated initially 5000 steps and then plotted in a definite interval as the interval is shifted.

If you run jupyter on your local host, try to run the below code after uncommenting the bottom line "plt.show". Otherwise, you can try to view inline mpeg4 animation by using mpeg encoder "ffmpeg" installed in the host of jupyter.

%matplotlib auto
# coding: UTF-8
"""
Lorenz model visualization as an exercise on animation
    continue to plot until the run quits.
"""


import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation

def lorentz_deriv(x, t0, sigma=10., beta=8./3, rho=28.0):
    """Compute the time-derivative of a Lorentz system."""
    return [sigma * (x[1] - x[0]), x[0] * (rho - x[2]) - x[1], x[0] * x[1] - beta * x[2]]

# Create new Figure and an Axes which fills it.
fig = plt.figure()
ax = fig.add_axes([0, 0, 1, 1], projection='3d')
# prepare the axes limits
ax.set_xlim((-25, 25))
ax.set_ylim((-35, 35))
ax.set_zlim((5, 55))

# set point-of-view: specified by (altitude degrees, azimuth degrees)
# ax.view_init(60, 45)

# initial setting
t = np.linspace(0,0.1,10)
N_trajectories = 2

x_ini = [[1., 1., 1.],[1.00001, 1., 1.] ]  # starting vector:  two traces with different initial value
line1, = ax.plot([], [], [], 'b-')
line2, = ax.plot([], [], [], 'r-')
line = [line1, line2]
pt1, = ax.plot([],[],[], 'bo')
pt2, = ax.plot([],[],[], 'ro')
pt = [pt1, pt2]

data4plot = [integrate.odeint(lorentz_deriv, x_ini[0],t),
             integrate.odeint(lorentz_deriv, x_ini[1],t)]

#update
def update(i) :
    global line, pt, x_ini, t, data4plot
    for n in range(N_trajectories):
        x_t = integrate.odeint(lorentz_deriv, x_ini[n],t)
        data4plot[n] = np.vstack((data4plot[n], x_t))
        line[n].set_data(data4plot[n][:,0],data4plot[n][:,1])         # trajectory in 3D
        line[n].set_3d_properties(data4plot[n][:,2]) 
        pt[n].set_data(data4plot[n][-1:,0], data4plot[n][-1:,1])    # head (current point)
        pt[n].set_3d_properties(data4plot[n][-1:,2])
        x_ini[n] = x_t[-1]
        
    ax.view_init(30, 0.3 * i)            # rotating eye's point
    return line                  # return value is a tuple of graphic objects. So we should add comma here for a single trajectory

anim = animation.FuncAnimation(fig, update,
                               frames=500, interval=100,blit=True)

# To make movie file, use anim.save function like the below line
# anim.save('lorentz3d.mp4', fps=15, extra_args=['-vcodec', 'libx264'])

# To make movie as a moving "gif" file, use the below (this will be huge)
# anim.save('lorentz3d.gif', writer='imagemagick',  fps=30)

# plt.show()

If you would like to view this animation in line, try the below code. To run this, you have to have installed ffmpeg.

# It takes somewhat long time until the movie start because it starts after completing the mpeg4 file for a given number of frames
from IPython.display import HTML
HTML(anim.to_html5_video())