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 [ ]: