In [7]:
from matplotlib import rcdefaults
rcdefaults()  # restore default matplotlib rc parameters
%config InlineBackend.figure_format='retina'   
import seaborn as sns  # wrapper for matplotlib that provides prettier styles and more
import matplotlib.pyplot as plt  # use matplotlib functionality directly
%matplotlib inline   
sns.set()
In [8]:
import numpy as np

Euler and Midpoint methods compared

In [84]:
def f(t,y):
    u,v = y
    return np.array([-v,u])

def exact(t):
    return 3* np.array(  [np.cos(t),np.sin(t)]  )

T = 2*np.pi  # final value of t

fig,axes = plt.subplots(3,1,figsize=(5,10))

def euler(m):
    t = 0
    y = np.array([3.,0.])

    h = (T-t)/m

    ta = np.empty(m+1)
    ta[0] = t

    ya = np.empty((2,m+1))  # a 2d array
    ya[:,0] = y
    for k in range(m): # 0, 1, 2, ...
        y = y + h*f(t,y)  # scalar mult and addition done element-wise
        t += h
        ta[k+1] = t
        ya[:,k+1] = y
    if doplot:
        axes[0].set_aspect(1)
        axes[0].plot(ya[0,:],ya[1,:],'.')
        axes[1].plot(ta,ya[0,:],'.')
        axes[2].plot(ta,ya[1,:],'.')
    return y

def midpoint(m):
    t = 0
    y = np.array([3.,0.])

    h = (T-t)/m

    ta = np.empty(m+1)
    ta[0] = t

    ya = np.empty((2,m+1))  # a 2d array
    ya[:,0] = y
    for k in range(m): # 0, 1, 2, ...
        ymid = y + h/2*f(t,y)  # scalar mult and addition done element-wise
        y = y + h*f(t+h/2,ymid)
        t += h
        ta[k+1] = t
        ya[:,k+1] = y

    if doplot:
        axes[0].set_aspect(1)
        axes[0].plot(ya[0,:],ya[1,:],'.')
        axes[1].plot(ta,ya[0,:],'.')
        axes[2].plot(ta,ya[1,:],'.')
    return y

doplot = True
for method in [euler,midpoint]:
    method(40)
t = np.linspace(0,2*np.pi,200)
y = exact(t)
axes[0].plot( y[0,:],y[1,:],alpha=0.5)
axes[1].plot( t,y[0,:],alpha=0.5)
axes[2].plot( t,y[1,:],alpha=0.5)
axes[0].set_xlabel('u'); axes[0].set_ylabel('v')
axes[1].set_xlabel('t'); axes[1].set_ylabel('u')
axes[2].set_xlabel('t'); axes[2].set_ylabel('v');

Midpoint (orange) is dramatically better than Euler (blue). Green is exact solution.

Error vs stepsize $h$

Goal is to plot the error at t=2 pi versus h

In [85]:
ytrue = [3,0]
plt.figure(figsize=(15,8))
plt.subplot(111,aspect=1)
doplot = False
for m in [40*2**i for i in range(8)]:
    h = 2*np.pi/m  # step size

    error = np.abs((euler(m) - ytrue)).max()
    print(m,error)
    plt.loglog(h,error,'.')

    error = np.abs((midpoint(m) - ytrue)).max()
    print(m,error)
    plt.loglog(h/2,error,'.')

    
    plt.xlabel('-work')
plt.ylabel('error infinity norm');
40 1.8783425143840855
40 0.07716297600176875
80 0.8363126695876031
80 0.019350209808586427
160 0.3935687933801377
160 0.0048427163683405355
320 0.19084286414227902
320 0.0012110497460797662
640 0.09396455425374928
640 0.00030278713962006085
1280 0.04662176425960629
1280 7.569837749178167e-05
2560 0.02322122520496661
2560 1.8924695438952976e-05
5120 0.01158825775073602
5120 4.7311802333070176e-06

We observe that for Euler the error is $O(h)$ (slope 1 on log-log plot), and for midpoint it is $O(h^2)$ (slope 2).

List comprehensions in python

In [18]:
[ i**2 for i in [1,2,3] ]
Out[18]:
[1, 4, 9]
In [19]:
[ 'hello' for i in [1,2,3] ]
Out[19]:
['hello', 'hello', 'hello']

A non-autonomous system

for a more stringent test of our methods

Same system as before but with a distorted t

In [81]:
gamma = 1.1

def f(t,y):
    u,v = y
    return np.array([-v*(1-np.cos(t)/gamma),u*(1-np.cos(t)/gamma)])

def exact(t):
    return np.array([ 3*np.cos( t - np.sin(t)/gamma), 3*np.sin( t - np.sin(t)/gamma ) ])

T = 2*np.pi  # final value of t

fig,axes = plt.subplots(3,1,figsize=(10,10))

def euler(m):
    t = 0
    y = np.array([3.,0.])

    h = (T-t)/m

    ta = np.empty(m+1)
    ta[0] = t

    ya = np.empty((2,m+1))  # a 2d array
    ya[:,0] = y
    for k in range(m): # 0, 1, 2, ...
        y = y + h*f(t,y)  # scalar mult and addition done element-wise
        t += h
        ta[k+1] = t
        ya[:,k+1] = y
    if doplot:
        axes[0].set_aspect(1)
        axes[0].plot(ya[0,:],ya[1,:],'.')
        axes[1].plot(ta,ya[0,:],'.')
        axes[2].plot(ta,ya[1,:],'.')
    return y

def midpoint(m):
    t = 0
    y = np.array([3.,0.])

    h = (T-t)/m

    ta = np.empty(m+1)
    ta[0] = t

    ya = np.empty((2,m+1))  # a 2d array
    ya[:,0] = y
    for k in range(m): # 0, 1, 2, ...
        ymid = y + h/2*f(t,y)  # scalar mult and addition done element-wise
        y = y + h*f(t+h/2,ymid)
        t += h
        ta[k+1] = t
        ya[:,k+1] = y

    if doplot:
        axes[0].set_aspect(1)
        axes[0].plot(ya[0,:],ya[1,:],'.')
        axes[1].plot(ta,ya[0,:],'.')
        axes[2].plot(ta,ya[1,:],'.')
    return y

doplot = True
for method in [euler,midpoint]:
    method(40)
t = np.linspace(0,2*np.pi,200)
y = exact(t)
axes[0].plot( y[0,:],y[1,:],alpha=0.5)
axes[1].plot( t,y[0,:],alpha=0.5)
axes[2].plot( t,y[1,:],alpha=0.5)
axes[0].set_xlabel('u'); axes[0].set_ylabel('v')
axes[1].set_xlabel('t'); axes[1].set_ylabel('u')
axes[2].set_xlabel('t'); axes[2].set_ylabel('v');

Orange is midpoint. Blue is Euler. Green is exact solution.

In [86]:
def f(t,y):
    u,v = y
    gamma = 1.2
    return np.array([-v*(1-np.cos(t)/gamma),u*(1-np.cos(t)/gamma)])

ytrue = [3,0]
plt.figure(figsize=(15,8))
plt.subplot(111,aspect=1)
doplot = False
for method in [euler,midpoint]:
    for m in [40*2**i for i in range(8)]:
        h = 2*np.pi/m  # step size

        error = np.abs((method(m) - ytrue)).max()
        print(m,error)
        plt.loglog(h,error,'.')
    
    plt.xlabel('h')
plt.ylabel('error infinity norm');
40 2.693105146223414
80 1.1713036233098988
160 0.5412787938318488
320 0.25983175547827386
640 0.1272657887041846
1280 0.06297766785210435
2560 0.031325984705819376
5120 0.015622402389822376
40 0.15616441940768194
80 0.03939063570217323
160 0.009878909565573545
320 0.002472006534789503
640 0.0006181547261827807
1280 0.0001545485896212828
2560 3.863777780412341e-05
5120 9.659484626687274e-06

For this non-autonomous system, we again confirm that for Euler the error is $O(h)$ (slope 1 on log-log plot), and for midpoint it is $O(h^2)$ (slope 2).

In [ ]: