Numerical Analysis

In class lab

The main point of today's lab is to try to make a groovy animation. We'll need ImageMagick so do go ahead and download it, if you haven't yet. And let's load some modules.

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import numpy as np

The dance of Earth and Venus

If we're going to make a movie, we might as well generate a cool one. Here's a cool picture:

In [2]:
earth_period = 365.256
venus_period = 224.7
venus_radius = 0.721925
def earth_position(t): 
    return (np.cos(2*np.pi*t/earth_period), np.sin(2*np.pi*t/earth_period))
def venus_position(t):
    return (
        venus_radius*np.cos(2*np.pi*t/venus_period), 
        venus_radius*np.sin(2*np.pi*t/venus_period))

fig,ax = plt.subplots()
earth, = ax.plot(0,0)
venus, = ax.plot(0,0)
ax.plot(0,0, 'yo', markersize=33)
ax.add_patch(Circle((0,0),1, fill=False))
ax.add_patch(Circle((0,0),venus_radius, fill=False))
ax.set_xlim(-1.1,1.1)
ax.set_ylim(-1.1,1.1)
ax.set_aspect(1)
ax.set_frame_on(False)
ax.set_xticks([])
ax.set_yticks([])
for i in range(0,int(8*365.25),6):
    ex,ey = earth_position(i)
    vx,vy = venus_position(i)
    plt.plot([ex,vx],[ey,vy], 'k-', alpha=0.1)
    earth.remove()
    venus.remove()
    earth, = plt.plot(ex,ey, 'bo')
    venus, = plt.plot(vx,vy, 'ro')
#    if 0<=i<10:
#        ext = '000'+str(i)
#    elif 10<=i<100:
#        ext = '00'+str(i)
#    elif 100<=i<1000:
#        ext = '0'+str(i)
#    else:
#        ext = str(i)
#    plt.savefig('pic'+ext+'.png')

The image that you see is just the last in a sequence of images that generate quite a lovely animation. If you uncomment the last nine lines, they should generate a sequence of images in your working directory. Once you've got that, you should be able to generate a movie with an ImageMagick command (executed from your terminal command line) like the following:

convert -delay 1 *.png -loop 0 anim.gif

After that, you can open anim.gif in your web browser to see it. You can play with the parameters that generate the image and the command line arguments of convert.

This example is called the dance of Earth and Venus. It assumes circular orbits and uniform speed and does not use differential equations at all. Perhaps, we should have a look at an example that does involve a differential equation?

The harmonic oscillator

Let's try to generate that animation from our Discourse problem: I'll get you started!

In [3]:
from scipy.integrate import odeint
def f(yyp,t): 
    y = yyp[0]
    yp = yyp[1]
    return [yp,-y]
ts = np.linspace(0,12,100)
sol = odeint(f,[0,1],ts)
y,yp = sol.T
def pic(i):
    plt.plot([-1.2,1.2],[0,0], 'k', linewidth=1)
    plt.plot([0,0],[-0.1,0.1], 'k', linewidth=0.5)
    plt.plot(y[i], 0, 'o')
    ax = plt.gca()
    ax.set_frame_on(False)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_aspect(1)
pic(50)

That much shows you how to set up the IVP in Python, solve it with odeint, and plot one of the resulting images. You've still got to interpolate the points, find the period, solve it again over the period, and animate the result. Sounds fun!