Pendulum motion with Forward Euler¶

AMath 586, Spring Quarter 2019 at the University of Washington. For other notebooks, see Index.ipynb or the Index of all notebooks on Github.

This notebook illustrates how Forward Euler behaves when modeling a simple pendulum on a coarse grid.

%matplotlib inline

from pylab import *
from matplotlib import animation
from IPython.display import HTML
from scipy.integrate import solve_ivp

Compute a fine grid solution for comparison:¶

def f(t,u):
    f0 = u[1]
    f1 = -sin(u[0])
    return array([f0,f1])

theta0 = pi/2.
u0 = array([theta0, 0.])
tf = 30.

t_span = (0., tf)
t_eval = linspace(0, tf, 1000)
solution = solve_ivp(f, t_span, u0, method='RK45', t_eval=t_eval, 
                      atol=1e-9, rtol=1e-9)
theta_eval = solution.y[0,:]

Forward Euler implementation:¶

def ForwardEuler(nsteps):
    t0 = 0.;  tfinal = 30.
    t = linspace(t0, tfinal, nsteps+1)
    dt = t[1] - t[0]
    
    # Array for computed solution
    # give it two rows so each column is solution at one time,
    
    U = empty((2,nsteps+1))
    U.fill(nan)
    U[0,0] = theta0  # initial angle theta
    U[1,0] = 0.  # initial angular velocity
    for n in range(0,nsteps):
        U[0,n+1] = U[0,n] + dt * U[1,n]
        U[1,n+1] = U[1,n] - dt * sin(U[0,n])
    
    return t,U

t,U = ForwardEuler(300)

plot(t,U[0,:],'b',label='computed')
plot(t_eval,theta_eval,'k',label='from solve_ivp')
legend()
print('Error at the final time = %g' % abs(U[0,-1]-theta_eval[-1]))

Error at the final time = 18.1484

Note that the Forward Euler solution for $\theta(t)$ grows in amplitude with time in a nonphysical manner and eventually goes over the top ($\theta > \pi$) and continues swinging around and around rather than osciallating about $\theta = 0$ as it should.

This is much easier to visualize in an animation...

Animation¶

The rest of this notebook shows how to make an animation of the solution.

This uses the animation.FuncAnimation function from matplotlib to animate the pendulum motion. See e.g. this page for a nice simple introductory example.

def pend_animation(nsteps, plot_interval=1):
    """
    Solve the pendulum problem with the method defined above using nsteps time steps.
    Create an animation showing the swinging pendulum, using the solution every
    plot_interval time steps.
    """
   
    # compute the solution with the method define above:
    t,U = ForwardEuler(nsteps)
    
    fig = plt.figure(figsize=(11,5))
    ax1 = plt.subplot(121)
    plot(t,U[0,:],'r')
    pi_ticks = pi*array([-1,0,1,2,3,4,5,6])
    pi_labels = ['$-\pi$','0','$\pi$','$2\pi$','$3\pi$','$4\pi$','$5\pi$','$6\pi$']
    ax1.set_yticks(pi_ticks)
    ax1.set_yticklabels(pi_labels)
    ax1.set_xlabel('time')
    ax1.set_ylabel('theta')
    ax1.grid(True)
    
    theta = U[0,0]
    pendlength = 2.
    xpend = pendlength*sin(theta);
    ypend = -pendlength*cos(theta);
    
    ax2 = plt.subplot(122,xlim=(-2.3, 2.3), ylim=(-2.3, 2.3))
    ax2.plot([-2.3,2.3],[0,0],'k',linewidth=0.5)
    ax2.plot([0,0],[-2.3,2.3],'k',linewidth=0.5)
    axis('scaled')
    ax2.set_xlim(-2.3,2.3)
    ax2.set_ylim(-2.3,2.3)
    
    def init():
        shaft, = ax2.plot([0,xpend], [0,ypend], linestyle='-', color='lightblue', lw=2)
        bob, = ax2.plot([xpend],[ypend],'o',color='lightblue', markersize=8)
        theta_pt, = ax1.plot([t[0]], [U[0,0]], 'o',color='pink')
        return (shaft,bob,theta_pt)

    shaft,bob,theta_pt = init()
    
    def fplot(n):
        theta = U[0,n]
        xpend = pendlength*sin(theta);
        ypend = -pendlength*cos(theta);
        shaft.set_data([0,xpend], [0,ypend])
        shaft.set_color('b')
        bob.set_data([xpend], [ypend])
        bob.set_color('b')
        theta_pt.set_data([t[n]], [U[0,n]])
        theta_pt.set_color('k')
        return (shaft,bob,theta_pt)

    frames_to_plot = range(0, len(t), plot_interval)
    close(fig)
    return animation.FuncAnimation(fig, fplot, frames=frames_to_plot, interval=100,
                                   blit=True,init_func=init,repeat=False)

anim = pend_animation(300,3);

HTML(anim.to_jshtml())