Figure 1: The Forward Euler scheme for two values of the time step.
Implementation
Implementation in Python
Making a program
Function for computing the numerical solution
Integer division
Doc strings
Formatting of numbers
Running the program
Verifying the implementation
Simplest method: run a few algorithmic steps by hand
Comparison with an exact discrete solution
Computing the numerical error
Plotting solutions
Plotting with SciTools
Creating user interfaces
Reading a sequence of command-line arguments
Working with an argument parser
Computing convergence rates
Estimating the convergence rate \( r \)
Implementation
Verification
Debugging via convergence rates
Memory-saving implementation
Software engineering
Making a module
Prefixing imported functions by the module name
Doctests
Unit testing with nose
Basic use of nose
Demonstrating nose
Classical unit testing with unittest
Basic use of unittest
Demonstration of unittest
Implementing simple problem and solver classes
The problem class
The solver class
The visualizer class
Combing the objects
Implementing more advanced problem and solver classes
A generic class for parameters
The problem class
The solver class
The visualizer class
Performing scientific experiments
Interpreting output from other programs
Making a report
Publishing a complete project
Model extensions
Extension to a variable coefficient
Extension to a source term
Schemes
Implementation of the generalized model problem
The Python code
Implementations of variable coefficients
Verification via trivial solutions
Verification via manufactured solutions
Extension to systems of ODEs
General first-order ODEs
Generic form
The Odespy software
Example: Runge-Kutta methods
Example: Adaptive Runge-Kutta methods
Other schemes
Implicit 2-step backward scheme
The Leapfrog scheme
The filtered Leapfrog scheme
2nd-order Runge-Kutta scheme
2nd-order Adams-Bashforth scheme
3rd-order Adams-Bashforth scheme
We address the differential equation problem $$ u'(t) = -au(t),\quad t\in (0,T], \quad u(0)=I, $$ where \( a>0 \) is a given constant.
The problem is solved by the numerical method
$$ u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n, $$ for \( \theta\in [0,1] \). The classical Forward Euler, Backward Euler, and Crank-Nicolson (midpoint) schemes are recovered through \( \theta=0,1,0.5 \) (resp.)
Tasks: Make a program that
Why Python?
from numpy import *
def solver(I, a, T, dt, theta):
"""Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
N = int(T/dt) # no of time intervals
T = N*dt # adjust T to fit time step dt
u = zeros(N+1) # array of u[n] values
t = linspace(0, T, N+1) # time mesh
u[0] = I # assign initial condition
for n in range(0, N): # n=0,1,...,N-1
u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
return u, t
Note about the for loop: range(0, N, s) generates all integers from 0 to N in steps of s (default 1), but not including N (!).
Sample call:
u, t = solver(I=1, a=2, T=8, dt=0.8, theta=1)
Python applies integer division: 1/2 is 0, 1./2 or 1.0/2 or 1/2. or 1/2.0 or 1.0/2.0 gives 0.5.
A safer solver function:
from numpy import *
def solver(I, a, T, dt, theta):
"""Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt."""
dt = float(dt) # avoid integer division
N = int(round(T/dt)) # no of time intervals
T = N*dt # adjust T to fit time step dt
u = zeros(N+1) # array of u[n] values
t = linspace(0, T, N+1) # time mesh
u[0] = I # assign initial condition
for n in range(0, N): # n=0,1,...,N-1
u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
return u, t
def solver(I, a, T, dt, theta):
"""
Solve
u'(t) = -a*u(t),
with initial condition u(0)=I, for t in the time interval
(0,T]. The time interval is divided into time steps of
length dt.
theta=1 corresponds to the Backward Euler scheme, theta=0
to the Forward Euler scheme, and theta=0.5 to the Crank-
Nicolson method.
"""
...
Can control formatting of reals and integers through the printf format:
print 't=%6.3f u=%g' % (t[i], u[i])
Or the alternative format string syntax:
print 't={t:6.3f} u={u:g}'.format(t=t[i], u=u[i])
Complete program file: dc_v1.py.
How to run the program:
Terminal> python dc_v1.py
Can also run it as "normal" Unix programs: ./programname:
Terminal> ./dc_v1.py
Use a calculator:
$$ A\equiv \frac{1 - (1-\theta) a\Delta t}{1 + \theta a \Delta t} = 0.298245614035$$ $$ \begin{align*} u^1 &= AI=0.0298245614035,\\ u^2 &= Au^1= 0.00889504462912,\\ u^3 &=Au^2= 0.00265290804728 \end{align*} $$
Compare hand-calculated numbers with those from solver:
def verify_three_steps():
"""Compare three steps with known manual computations."""
theta = 0.8; a = 2; I = 0.1; dt = 0.8
u_by_hand = array([I,
0.0298245614035,
0.00889504462912,
0.00265290804728])
N = 3 # number of time steps
u, t = solver(I=I, a=a, T=N*dt, dt=dt, theta=theta)
tol = 1E-15 # tolerance for comparing floats
difference = abs(u - u_by_hand).max()
success = difference <= tol
return success
def main():
u, t = solver(I=1, a=2, T=8, dt=0.8, theta=1)
# Write out a table of t and u values:
for i in range(len(t)):
print 't=%6.3f u=%g' % (t[i], u[i])
# or print 't={t:6.3f} u={u:g}'.format(t=t[i], u=u[i])
Force the verification to be done before a real computation:
if verify_three_steps():
main()
else:
print 'Bug in the implementation!'
Complete program: dc_verf1.py.
The exact discrete solution as $$ \begin{equation} u^n = IA^n \label{decay:un:exact} \thinspace . \end{equation} $$
Implementation:
def verify_exact_discrete_solution():
def exact_discrete_solution(n, I, a, theta, dt):
factor = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)
return I*factor**n
theta = 0.8; a = 2; I = 0.1; dt = 0.8
N = int(8/dt) # no of steps
u, t = solver(I=I, a=a, T=N*dt, dt=dt, theta=theta)
u_de = array([exact_discrete_solution(n, I, a, theta, dt)
for n in range(N+1)])
difference = abs(u_de - u).max() # max deviation
tol = 1E-15 # tolerance for comparing floats
success = difference <= tol
return success
The complete program: dc_verf2.py.
Task: compute exact minus numerical solution, \( \uex(t_n) - u^n \)
Exact solution: \( \uex(t)=Ie^{-at} \), implemented as
def exact_solution(t, I, a):
return I*exp(-a*t)
Compute
$$ \begin{equation} e_n = \uex(t_n) - u^n,\quad n=0,1,\ldots,N \thinspace . \end{equation} $$
by
u, t = solver(I, a, T, dt, theta) # Numerical solution
u_e = exact_solution(t, I, a)
e = u_e - u
Note: array arithmetics - we compute on entire arrays
$$ \hat E^2 = \Delta t\left(\frac{1}{2}e_0^2 + \frac{1}{2}e_N^2 + \sum_{n=1}^{N-1} e_n^2\right) $$ Common approximation (ok if used consistently in all error computations):
$$ \hat E^2 \approx E^2 = \Delta t\sum_{n=0}^{N} e_n^2 $$
\( E \) is the \( L_2 \) norm of the discrete error function \( e^n \).
$$ \begin{equation} E = \sqrt{\Delta t\sum_{n=0}^N e_n^2} \label{decay:E} \end{equation} $$
Python:
E = sqrt(dt*sum(e**2))
Similar element by element computation:
m = len(u) # length of u array (alt: u.size)
u_e = zeros(m)
t = 0
for i in range(m):
u_e[i] = exact_solution(t, a, I)
t = t + dt
e = zeros(m)
for i in range(m):
e[i] = u_e[i] - u[i]
s = 0 # summation variable
for i in range(m):
s = s + e[i]**2
error = sqrt(dt*s)
Such scalar computing
from matplotlib.pyplot import *
plot(t, u)
show()
Compare with \( \uex(t) \):
t_e = linspace(0, T, 1001) # fine mesh
u_e = exact_solution(t_e, I, a)
plot(t, u, 'r-') # red line for u
plot(t_e, u_e, 'b-') # blue line for u_e
Figure 1: The Forward Euler scheme for two values of the time step.
How to include such decorations:
from matplotlib.pyplot import *
figure() # create new plot
t_e = linspace(0, T, 1001) # fine mesh for u_e
u_e = exact_solution(t_e, I, a)
plot(t, u, 'r--o') # red dashes w/circles
plot(t_e, u_e, 'b-') # blue line for exact sol.
legend(['numerical', 'exact'])
xlabel('t')
ylabel('u')
title('theta=%g, dt=%g' % (theta, dt))
savefig('%s_%g.png' % (theta, dt))
show()
Computation and plotting:
def explore(I, a, T, dt, theta=0.5, makeplot=True):
"""
Run a case with the solver, compute error measure,
and plot the numerical and exact solutions (if makeplot=True).
"""
u, t = solver(I, a, T, dt, theta) # Numerical solution
u_e = exact_solution(t, I, a)
e = u_e - u
E = sqrt(dt*sum(e**2))
if makeplot:
figure() # create new plot
t_e = linspace(0, T, 1001) # fine mesh for u_e
u_e = exact_solution(t_e, I, a)
plot(t, u, 'r--o') # red dashes w/circles
plot(t_e, u_e, 'b-') # blue line for exact sol.
legend(['numerical', 'exact'])
xlabel('t')
ylabel('u')
title('theta=%g, dt=%g' % (theta, dt))
theta2name = {0: 'FE', 1: 'BE', 0.5: 'CN'}
savefig('%s_%g.png' % (theta2name[theta], dt))
savefig('%s_%g.pdf' % (theta2name[theta], dt))
savefig('%s_%g.eps' % (theta2name[theta], dt))
show()
return E
Complete code: dc_plot_mpl.py.
Figure 2: The Backward Euler scheme for two values of the time step.
Figure 3: The Crank-Nicolson scheme for two values of the time step.
SciTools provides a unified plotting interface (Easyviz) to many different plotting packages: Matplotlib, Gnuplot, Grace, VTK, OpenDX, ...
Can use Matplotlib syntax, or a more compact plot function syntax:
from scitools.std import *
plot(t, u, 'r--o', # red dashes w/circles
t_e, u_e, 'b-', # blue line for exact sol.
legend=['numerical', 'exact'],
xlabel='t',
ylabel='u',
title='theta=%g, dt=%g' % (theta, dt),
savefig='%s_%g.png' % (theta2name[theta], dt),
show=True)
Complete code in dc_plot_st.py.
Change backend (plotting engine, Matplotlib by default):
Terminal> python dc_plot_st.py --SCITOOLS_easyviz_backend gnuplot
Terminal> python dc_plot_st.py --SCITOOLS_easyviz_backend grace
The program dc_plot_mpl.py needs this input: \( I \), \( a \), \( T \), an option to turn the plot on or off (makeplot), and a list of \( \Delta t \) values. Give these on the command line:
import sys
def read_command_line():
if len(sys.argv) < 6:
print 'Usage: %s I a T on/off dt1 dt2 dt3 ...' % \
sys.argv[0]; sys.exit(1) # abort
I = float(sys.argv[1])
a = float(sys.argv[2])
T = float(sys.argv[3])
makeplot = sys.argv[4] in ('on', 'True')
dt_values = [float(arg) for arg in sys.argv[5:]]
return I, a, T, makeplot, dt_values
Note:
Example: give \( a \) as --a 2.5 on the command line if the default value is not suitable.
def define_command_line_options():
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('--I', '--initial_condition', type=float,
default=1.0, help='initial condition, u(0)',
metavar='I')
parser.add_argument('--a', type=float,
default=1.0, help='coefficient in ODE',
metavar='a')
parser.add_argument('--T', '--stop_time', type=float,
default=1.0, help='end time of simulation',
metavar='T')
parser.add_argument('--makeplot', action='store_true',
help='display plot or not')
parser.add_argument('--dt', '--time_step_values', type=float,
default=[1.0], help='time step values',
metavar='dt', nargs='+', dest='dt_values')
return parser
Automatically generated help functionality:
Terminal> python dc_argparse.py -h
...
--I I, --initial_condition I
initial condition, u(0)
...
Structure of this help output:
--I metavar, --initial_condition metavar
help-string
How to read data from the command line:
def read_command_line():
parser = define_command_line_options()
args = parser.parse_args()
print 'I={}, a={}, T={}, makeplot={}, dt_values={}'.format(
args.I, args.a, args.T, args.makeplot, args.dt_values)
return args.I, args.a, args.T, args.makeplot, args.dt_values
Complete program: dc_argparse.py.
Frequent assumption on the relation between the numerical error \( E \) and some discretization parameter \( \Delta t \):
$$ \begin{equation} E = C\Delta t^r, \label{decay:E:dt} \end{equation} $$ Unknown: \( C \) and \( r \).
Perform numerical experiments: \( (\Delta t_i, E_i) \), \( i=0,\ldots,m-1 \).
Method 2 is best.
Compute \( r_0, r_1, \ldots, r_{m-2} \):
from math import log
def main():
I, a, T, makeplot, dt_values = read_command_line()
r = {} # estimated convergence rates
for theta in 0, 0.5, 1:
E_values = []
for dt in dt_values:
E = explore(I, a, T, dt, theta, makeplot=False)
E_values.append(E)
# Compute convergence rates
m = len(dt_values)
r[theta] = [log(E_values[i-1]/E_values[i])/
log(dt_values[i-1]/dt_values[i])
for i in range(1, m, 1)]
for theta in r:
print '\nPairwise convergence rates for theta=%g:' % theta
print ' '.join(['%.2f' % r_ for r_ in r[theta]])
return r
Complete program: dc_convrate.py.
Example on execution:
Terminal> python dc_convrate.py --dt 0.5 0.25 0.1 0.05 0.025 0.01
...
Pairwise convergence rates for theta=0:
1.33 1.15 1.07 1.03 1.02
Pairwise convergence rates for theta=0.5:
2.14 2.07 2.03 2.01 2.01
Pairwise convergence rates for theta=1:
0.98 0.99 0.99 1.00 1.00
Strong verification method: verify that \( r \) has the expected value!
Potential bug: missing a in the denominator,
u[n+1] = (1 - (1-theta)*a*dt)/(1 + theta*dt)*u[n]
Running dc_convrate.py gives same rates.
Why? The value of \( a \)... (\( a=1 \))
0 and 1 are bad values in tests!
Better:
Terminal> python dc_convrate.py --a 2.1 --I 0.1 \
--dt 0.5 0.25 0.1 0.05 0.025 0.01
...
Pairwise convergence rates for theta=0:
1.49 1.18 1.07 1.04 1.02
Pairwise convergence rates for theta=0.5:
-1.42 -0.22 -0.07 -0.03 -0.01
Pairwise convergence rates for theta=1:
0.21 0.12 0.06 0.03 0.01
Forward Euler works...because \( \theta=0 \) hides the bug.
This bug gives \( r\approx 0 \):
u[n+1] = ((1-theta)*a*dt)/(1 + theta*dt*a)*u[n]
def solver_memsave(I, a, T, dt, theta, filename='sol.dat'):
"""
Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt.
Minimum use of memory. The solution is store on file
(with name filename) for later plotting.
"""
dt = float(dt) # avoid integer division
N = int(round(T/dt)) # no of intervals
outfile = open(filename, 'w')
# u: time level n+1, u_1: time level n
t = 0
u_1 = I
outfile.write('%.16E %.16E\n' % (t, u_1))
for n in range(1, N+1):
u = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)*u_1
u_1 = u
t += dt
outfile.write('%.16E %.16E\n' % (t, u))
outfile.close()
return u, t
Reading the computed data:
def read_file(filename='sol.dat'):
infile = open(filename, 'r')
u = []; t = []
for line in infile:
words = line.split()
if len(words) != 2:
print 'Found more than two numbers on a line!', words
sys.exit(1) # abort
t.append(float(words[0]))
u.append(float(words[1]))
return np.array(t), np.array(u)
Simpler code with numpy functionality for reading/writing tabular data:
def read_file_numpy(filename='sol.dat'):
data = np.loadtxt(filename)
t = data[:,0]
u = data[:,1]
return t, u
Similar function np.savetxt, but then we need all \( u^n \) and \( t^n \) values in a two-dimensional array (which we try to prevent now!).
Usage:
def explore(I, a, T, dt, theta=0.5, makeplot=True):
filename = 'u.dat'
u, t = solver_minmem(I, a, T, dt, theta, filename)
t, u = read_file(filename)
u_e = exact_solution(t, I, a)
e = u_e - u
E = np.sqrt(dt*np.sum(e**2))
if makeplot:
plt.figure()
...
Complete program: dc_memsave.py.
Goal: make more professional numerical software.
Topics:
Sketch:
from numpy import *
from matplotlib.pyplot import *
import sys
def solver(I, a, T, dt, theta):
...
def verify_three_steps():
...
def verify_exact_discrete_solution():
...
def exact_solution(t, I, a):
...
def explore(I, a, T, dt, theta=0.5, makeplot=True):
...
def define_command_line_options():
...
def read_command_line(use_argparse=True):
...
def main():
...
That is! It's a module dc_mod in file dc_mod.py.
Usage in some other program:
from dc_mod import solver
u, t = solver(I=1.0, a=3.0, T=3, dt=0.01, theta=0.5)
Test block:
if __name__ == '__main__':
main()
If dc_mod is imported, __name__ is dc_mod.
If dc_mod.py is run, __name__ is __main__.
Use test block for testing, demo, user interface, ...
Extended test block:
if __name__ == '__main__':
if 'verify' in sys.argv:
if verify_three_steps() and verify_discrete_solution():
pass # ok
else:
print 'Bug in the implementation!'
elif 'verify_rates' in sys.argv:
sys.argv.remove('verify_rates')
if not '--dt' in sys.argv:
print 'Must assign several dt values'
sys.exit(1) # abort
if verify_convergence_rate():
pass
else:
print 'Bug in the implementation!'
else:
# Perform simulations
main()
from numpy import *
from matplotlib.pyplot import *
This imports a large number of names (sin, exp, linspace, plot, ...).
Confusion: is a function from`numpy`? Or matplotlib.pyplot? Or is it our own function?
Alternative (recommended) import:
import numpy
import matplotlib.pyplot
Now we need to prefix functions with module name:
t = numpy.linspace(0, T, N+1)
u_e = I*numpy.exp(-a*t)
matplotlib.pyplot.plot(t, u_e)
Common standard:
import numpy as np
import matplotlib.pyplot as plt
t = np.linspace(0, T, N+1)
u_e = I*np.exp(-a*t)
plt.plot(t, u_e)
Downside: math line \( e^{-at}\sin(2\pi t) \) gets cluttered with module names,
numpy.exp(-a*t)*numpy.sin(2(numpy.pi*t)
# or
np.exp(-a*t)*np.sin(2*np.pi*t)
Doc strings can be equipped with interactive Python sessions for demonstrating usage and automatic testing of functions.
def solver(I, a, T, dt, theta):
"""
Solve u'=-a*u, u(0)=I, for t in (0,T] with steps of dt.
>>> u, t = solver(I=0.8, a=1.2, T=4, dt=0.5, theta=0.5)
>>> for t_n, u_n in zip(t, u):
... print 't=%.1f, u=%.14f' % (t_n, u_n)
t=0.0, u=0.80000000000000
t=0.5, u=0.43076923076923
t=1.0, u=0.23195266272189
t=1.5, u=0.12489758761948
t=2.0, u=0.06725254717972
t=2.5, u=0.03621291001985
t=3.0, u=0.01949925924146
t=3.5, u=0.01049960113002
t=4.0, u=0.00565363137770
"""
...
Automatic check that the code gives the same output as above:
Terminal> python -m doctest dc_mod_doctest.py
Report in case of failure:
Terminal> python -m doctest dc_mod_doctest.py
********************************************************
File "dc_mod_doctest.py", line 12, in dc_mod_doctest....
Failed example:
for t_n, u_n in zip(t, u):
print 't=%.1f, u=%.14f' % (t_n, u_n)
Expected:
t=0.0, u=0.80000000000000
t=0.5, u=0.43076923076923
t=1.0, u=0.23195266272189
t=1.5, u=0.12489758761948
t=2.0, u=0.06725254717972
t=2.5, u=0.03621291001985
t=3.0, u=0.01949925924146
t=3.5, u=0.01049960113002
t=4.0, u=0.00565363137770
Got:
t=0.0, u=0.80000000000000
t=0.5, u=0.43076923076923
t=1.0, u=0.23195266272189
t=1.5, u=0.12489758761948
t=2.0, u=0.06725254717972
********************************************************
1 items had failures:
1 of 2 in dc_mod_doctest.solver
***Test Failed*** 1 failures.
Limit the number of digits in the output in doctests! Otherwise, round-off errors on a different machine may ruin the test.
Another example on using doctests:
def explore(I, a, T, dt, theta=0.5, makeplot=True):
"""
Run a case with the solver, compute error measure,
and plot the numerical and exact solutions (if makeplot=True).
>>> for theta in 0, 0.5, 1:
... E = explore(I=1.9, a=2.1, T=5, dt=0.1, theta=theta,
... makeplot=False)
... print '%.10E' % E
...
7.3565079236E-02
2.4183893110E-03
6.5013039886E-02
"""
...
Complete program: dc_mod_doctest.py.
Avoid doctests in functions using sys.argv and print (possible, but needs careful coding).
def double(n):
return 2*n
Either in this file or in a separate file test_mymod.py, we implement a test function for checking that double works:
import nose.tools as nt
def test_double():
result = mymod.double(4)
nt.assert_equal(result, 8)
(Need import mymod if the test is in test_mymod.py.)
Running
Terminal> nosetests
makes the nose tool look for Python files with test_*() functions and run these functions.
Main point with a test function: raise AssertionError exception if the test fails.
Alternative:
if result != 8:
raise AssertionError()
Advantage of nose:
Back to the exponential decay problem and the solver function.
We design three unit tests:
import nose.tools as nt
import sys, os
sys.path.insert(0, os.pardir)
import dc_mod_unittest as dc_mod
import numpy as np
def exact_discrete_solution(n, I, a, theta, dt):
"""Return exact discrete solution of the theta scheme."""
dt = float(dt) # avoid integer division
factor = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)
return I*factor**n
def test_against_discrete_solution():
"""
Compare result from solver against
formula for the discrete solution.
"""
theta = 0.8; a = 2; I = 0.1; dt = 0.8
N = int(8/dt) # no of steps
u, t = dc_mod.solver(I=I, a=a, T=N*dt, dt=dt, theta=theta)
u_de = np.array([exact_discrete_solution(n, I, a, theta, dt)
for n in range(N+1)])
diff = np.abs(u_de - u).max()
nt.assert_almost_equal(diff, 0, delta=1E-14)
nt.assert_almost_equal: compare two floats and avoid differences due to round-off errors.
def test_solver():
"""
Compare result from solver against
precomputed arrays for theta=0, 0.5, 1.
"""
I=0.8; a=1.2; T=4; dt=0.5 # fixed parameters
precomputed = {
't': np.array([ 0. , 0.5, 1. , 1.5, 2. , 2.5,
3. , 3.5, 4. ]),
0.5: np.array(
[ 0.8 , 0.43076923, 0.23195266, 0.12489759,
0.06725255, 0.03621291, 0.01949926, 0.0104996 ,
0.00565363]),
0: np.array(
[ 8.00000000e-01, 3.20000000e-01,
1.28000000e-01, 5.12000000e-02,
2.04800000e-02, 8.19200000e-03,
3.27680000e-03, 1.31072000e-03,
5.24288000e-04]),
1: np.array(
[ 0.8 , 0.5 , 0.3125 , 0.1953125 ,
0.12207031, 0.07629395, 0.04768372, 0.02980232,
0.01862645]),
}
for theta in 0, 0.5, 1:
u, t = dc_mod.solver(I, a, T, dt, theta=theta)
diff = np.abs(u - precomputed[theta]).max()
# Precomputed numbers are known to 8 decimal places
nt.assert_almost_equal(diff, 0, places=8,
msg='theta=%s' % theta)
theta = 1; a = 1; I = 1; dt = 2
lead to integer division:
(1 - (1-theta)*a*dt) # becomes 1
(1 + theta*dt*a) # becomes 2
(1 - (1-theta)*a*dt)/(1 + theta*dt*a) # becomes 0 (!)
Unit test for this issue:
def test_potential_integer_division():
"""Choose variables that can trigger integer division."""
theta = 1; a = 1; I = 1; dt = 2
N = 4
u, t = dc_mod.solver(I=I, a=a, T=N*dt, dt=dt, theta=theta)
u_de = np.array([exact_discrete_solution(n, I, a, theta, dt)
for n in range(N+1)])
diff = np.abs(u_de - u).max()
nt.assert_almost_equal(diff, 0, delta=1E-14)
Test convergence rates:
def test_convergence_rates():
"""Compare empirical convergence rates to exact ones."""
# Set command-line arguments directly in sys.argv
sys.argv[1:] = '--I 0.8 --a 2.1 --T 5 '\
'--dt 0.4 0.2 0.1 0.05 0.025'.split()
# Suppress output from dc_mod.main()
stdout = sys.stdout # save standard output for later use
scratchfile = open('.tmp', 'w') # fake standard output
sys.stdout = scratchfile
r = dc_mod.main()
for theta in r:
nt.assert_true(r[theta]) # check for non-empty list
scratchfile.close()
sys.stdout = stdout # restore standard output
expected_rates = {0: 1, 1: 1, 0.5: 2}
for theta in r:
r_final = r[theta][-1]
# Compare to 1 decimal place
nt.assert_almost_equal(expected_rates[theta], r_final,
places=1, msg='theta=%s' % theta)
# no need for any main
Complete program: test_dc_nose.py.
Write file test_mymod.py:
import unittest
import mymod
class TestMyCode(unittest.TestCase):
def test_double(self):
result = mymod.double(4)
self.assertEqual(result, 8)
if __name__ == '__main__':
unittest.main()
import unittest
import dc_mod_unittest as decay
import numpy as np
def exact_discrete_solution(n, I, a, theta, dt):
factor = (1 - (1-theta)*a*dt)/(1 + theta*dt*a)
return I*factor**n
class TestDecay(unittest.TestCase):
def test_against_discrete_solution(self):
...
diff = np.abs(u_de - u).max()
self.assertAlmostEqual(diff, 0, delta=1E-14)
def test_solver(self):
...
for theta in 0, 0.5, 1:
...
self.assertAlmostEqual(diff, 0, places=8,
msg='theta=%s' % theta)
def test_potential_integer_division():
...
self.assertAlmostEqual(diff, 0, delta=1E-14)
def test_convergence_rates(self):
...
for theta in r:
...
self.assertAlmostEqual(...)
if __name__ == '__main__':
unittest.main()
Complete program: test_dc_unittest.py.
from numpy import exp
class Problem:
def __init__(self, I=1, a=1, T=10):
self.T, self.I, self.a = I, float(a), T
def exact_solution(self, t):
I, a = self.I, self.a # extract local variables
return I*exp(-a*t)
Basic usage:
problem = Problem(T=5)
problem.T = 8
problem.dt = 1.5
Improvement: read data from the command line.
class Problem:
def __init__(self, I=1, a=1, T=10):
self.T, self.I, self.a = I, float(a), T
def define_command_line_options(self, parser=None):
if parser is None:
import argparse
parser = argparse.ArgumentParser()
parser.add_argument(
'--I', '--initial_condition', type=float,
default=self.I, help='initial condition, u(0)',
metavar='I')
parser.add_argument(
'--a', type=float, default=self.a,
help='coefficient in ODE', metavar='a')
parser.add_argument(
'--T', '--stop_time', type=float, default=self.T,
help='end time of simulation', metavar='T')
return parser
def init_from_command_line(self, args):
self.I, self.a, self.T = args.I, args.a, args.T
def exact_solution(self, t):
I, a = self.I, self.a
return I*exp(-a*t)
Observe
class Solver:
def __init__(self, problem, dt=0.1, theta=0.5):
self.problem = problem
self.dt, self.theta = float(dt), theta
def define_command_line_options(self, parser):
parser.add_argument(
'--dt', '--time_step_value', type=float,
default=0.5, help='time step value', metavar='dt')
parser.add_argument(
'--theta', type=float, default=0.5,
help='time discretization parameter', metavar='dt')
return parser
def init_from_command_line(self, args):
self.dt, self.theta = args.dt, args.theta
def solve(self):
from dc_mod import solver
self.u, self.t = solver(
self.problem.I, self.problem.a, self.problem.T,
self.dt, self.theta)
def error(self):
u_e = self.problem.exact_solution(self.t)
e = u_e - self.u
E = sqrt(self.dt*sum(e**2))
return E
Note: reuse of the numerical algorithm from the dc_mod module.
class Visualizer:
def __init__(self, problem, solver):
self.problem, self.solver = problem, solver
def plot(self, include_exact=True, plt=None):
"""
Add solver.u curve to scitools plotting object plt,
and include the exact solution if include_exact is True.
This plot function can be called several times (if
the solver object has computed new solutions).
"""
if plt is None:
import scitools.std
plt = scitools.std
plt.plot(self.solver.t, self.solver.u, '--o')
plt.hold('on')
theta2name = {0: 'FE', 1: 'BE', 0.5: 'CN'}
name = self.solver.theta2name.get(self.solver.theta, '')
plt.legend('numerical %s' % name)
if include_exact:
t_e = linspace(0, self.problem.T, 1001)
u_e = self.problem.exact_solution(t_e)
plt.plot(t_e, u_e, 'b-')
plt.legend('exact')
plt.xlabel('t')
plt.ylabel('u')
plt.title('theta=%g, dt=%g' %
(self.solver.theta, self.solver.dt))
plt.savefig('%s_%g.png' % (name, self.solver.dt))
return plt
Note:
Let Problem, Solver, and Visualizer play together:
def main():
problem = Problem()
solver = Solver(problem)
viz = Visualizer(problem, solver)
# Read input from the command line
parser = problem.define_command_line_options()
parser = solver. define_command_line_options(parser)
args = parser.parse_args()
problem.init_from_command_line(args)
solver. init_from_command_line(args)
# Solve and plot
solver.solve()
plt = viz.plot()
E = solver.error()
if E is not None:
print 'Error: %.4E' % E
plt.show()
Complete program: dc_class.py.
class Parameters:
def set(self, **parameters):
for name in parameters:
self.prms[name] = parameters[name]
def get(self, name):
return self.prms[name]
def define_command_line_options(self, parser=None):
if parser is None:
import argparse
parser = argparse.ArgumentParser()
for name in self.prms:
tp = self.types[name] if name in self.types else str
help = self.help[name] if name in self.help else None
parser.add_argument(
'--' + name, default=self.get(name), metavar=name,
type=tp, help=help)
return parser
def init_from_command_line(self, args):
for name in self.prms:
self.prms[name] = getattr(args, name)
Slightly more advanced version in class_dc_verf1.py.
class Problem(Parameters):
"""
Physical parameters for the problem u'=-a*u, u(0)=I,
with t in [0,T].
"""
def __init__(self):
self.prms = dict(I=1, a=1, T=10)
self.types = dict(I=float, a=float, T=float)
self.help = dict(I='initial condition, u(0)',
a='coefficient in ODE',
T='end time of simulation')
def exact_solution(self, t):
I, a = self.get('I'), self.get('a')
return I*np.exp(-a*t)
class Solver(Parameters):
def __init__(self, problem):
self.problem = problem
self.prms = dict(dt=0.5, theta=0.5)
self.types = dict(dt=float, theta=float)
self.help = dict(dt='time step value',
theta='time discretization parameter')
def solve(self):
#from dc_mod import solver
from decay_theta import solver
self.u, self.t = solver(
self.problem.get('I'),
self.problem.get('a'),
self.problem.get('T'),
self.get('dt'),
self.get('theta'))
def error(self):
try:
u_e = self.problem.exact_solution(self.t)
e = u_e - self.u
E = np.sqrt(self.get('dt')*np.sum(e**2))
except AttributeError:
E = None
return E
Goal: explore the behavior of a numerical method for a differential equation and show how scientific experiments can be set up and reported.
Tasks:
$$ \begin{equation} u'(t) = -au(t),\quad u(0)=I,\ 0< t \leq T, \label{decay:experiments:model} \end{equation} $$
Solution method (\( \theta \)-rule):
$$ u^{n+1} = \frac{1 - (1-\theta) a\Delta t}{1 + \theta a\Delta t}u^n, \quad u^0=I\thinspace . $$
Tasks:
Figure 4: Illustration of the Backward Euler method for four time step values.
import os, sys
# The command line must contain dt values
if len(sys.argv) > 1:
dt_values = [float(arg) for arg in sys.argv[1:]]
else:
print 'Usage: %s dt1 dt2 dt3 ...'; sys.exit(1) # abort
# Fixed physical parameters
I = 1
a = 2
T = 5
# Run module file as a stand-alone application
cmd = 'python dc_mod.py --I %g --a %g --makeplot --T %g' % \
(I, a, T)
dt_values_str = ' '.join([str(v) for v in dt_values])
cmd += ' --dt %s' % dt_values_str
print cmd
failure = os.system(cmd)
if failure:
print 'Command failed:', cmd; sys.exit(1)
# Combine images into rows with 2 plots in each row
combine_image_commands = []
for method in 'BE', 'CN', 'FE':
imagefiles = ' '.join(['%s_%s.pdf' % (method, dt)
for dt in dt_values])
combine_image_commands.append(
'montage -background white -geometry 100%' + \
' -tile 2x %s %s.png' % (imagefiles, method))
combine_image_commands.append(
'convert -trim %s.png %s.png' % (method, method))
combine_image_commands.append(
'pdftk %s output tmp.pdf' % imagefiles)
num_rows = int(round(len(dt_values)/2.0))
combine_image_commands.append(
'pdfnup --nup 2x%d tmp.pdf' % num_rows)
combine_image_commands.append(
'pdfcrop tmp-nup.pdf')
combine_image_commands.append(
'mv -f tmp-nup-crop.pdf %s.pdf;'
'rm -f tmp.pdf tmp-nup.pdf' % method)
imagefiles = ' '.join(['%s_%s.png' % (method, dt)
for dt in dt_values])
for cmd in combine_image_commands:
print cmd
failure = os.system(cmd)
if failure:
print 'Command failed:', cmd; sys.exit(1)
# Remove the files generated by dc_mod.py
from glob import glob
filenames = glob('*_*.png') + glob('*_*.pdf') + glob('*_*.eps') +\
glob('tmp*.pdf')
for filename in filenames:
os.remove(filename)
Complete program: experiments/dc_exper0.py.
Many useful constructs in the program above:
Programs that run other programs, like dc_exper0.py does, will often need to interpret output from the other programs. For example,
Terminal> python dc_plot_mpl.py
0.0 0.40: 2.105E-01
0.0 0.04: 1.449E-02
0.5 0.40: 3.362E-02
0.5 0.04: 1.887E-04
1.0 0.40: 1.030E-01
1.0 0.04: 1.382E-02
Tasks:
from subprocess import Popen, PIPE, STDOUT
p = Popen(cmd, shell=True, stdout=PIPE, stderr=STDOUT)
output, dummy = p.communicate()
failure = p.returncode
if failure:
print 'Command failed:', cmd; sys.exit(1)
Note: The command stored in cmd is run and all text that is written to the standard output and the standard error is available in the string output. The text in output is what appeared in the terminal window while running cmd.
Next tasks:
errors = {'dt': dt_values, 1: [], 0: [], 0.5: []}
for line in output.splitlines():
words = line.split()
if words[0] in ('0.0', '0.5', '1.0'): # line with E?
# typical line: 0.0 1.25: 7.463E+00
theta = float(words[0])
E = float(words[2])
errors[theta].append(E)
import matplotlib.pyplot as plt
#import scitools.std as plt
plt.loglog(errors['dt'], errors[0], 'ro-')
plt.hold('on') # MATLAB style...
plt.loglog(errors['dt'], errors[0.5], 'b+-')
plt.loglog(errors['dt'], errors[1], 'gx-')
plt.legend(['FE', 'CN', 'BE'], loc='upper left')
plt.xlabel('log(time step)')
plt.ylabel('log(error)')
plt.title('Error vs time step')
plt.savefig('error_BE_CN_FE.png')
This is a log-log plot because we expect \( E\sim\Delta t^r \).
Figure 5: Default plot (left) and manually adjusted plot (right).
Complete program: experiments/dc_exper1.py. Fine recipe for
$$ \begin{equation} u'(t) = -a(t)u(t),\quad t\in (0,T],\quad u(0)=I \thinspace . \label{decay:problem:a} \end{equation} $$
The Forward Euler scheme:
$$ \begin{equation} \frac{u^{n+1} - u^n}{\Delta t} = -a(t_n)u^n \thinspace . \end{equation} $$
The Backward Euler scheme: $$ \begin{equation} \frac{u^{n} - u^{n-1}}{\Delta t} = -a(t_n)u^n \thinspace . \end{equation} $$
The Crank-Nicolson scheme (evaluting \( a(t_{n+\frac{1}{2}}) \) and using an average for \( u \)): $$ \begin{equation} \frac{u^{n+1} - u^{n}}{\Delta t} = -a(t_{n+\frac{1}{2}})\frac{1}{2}(u^n + u^{n+1}) \thinspace . \end{equation} $$
The Crank-Nicolson scheme (using using an average for \( a \) and \( u \)): $$ \begin{equation} \frac{u^{n+1} - u^{n}}{\Delta t} = -\frac{1}{2}(a(t_n)u^n + a(t_{n+1})u^{n+1}) \thinspace . \end{equation} $$
The \( \theta \)-rule unifies the three mentioned schemes,
$$ \begin{equation} \frac{u^{n+1} - u^{n}}{\Delta t} = -a((1-\theta)t_n + \theta t_{n+1})((1-\theta) u^n + \theta u^{n+1}) \thinspace . \end{equation} $$ or, $$ \begin{equation} \frac{u^{n+1} - u^{n}}{\Delta t} = -(1-\theta) a(t_n)u^n - \theta a(t_{n+1})u^{n+1} \thinspace . \end{equation} $$
Operator notation:
$$ \begin{align*} \lbrack D^+_t u &= -au\rbrack^n,\\ \lbrack D^-_t u &= -au\rbrack^n,\\ \lbrack D_t u &= -a\overline{u}^t\rbrack^{n+\frac{1}{2}},\\ \lbrack D_t u &= -\overline{au}^t\rbrack^{n+\frac{1}{2}},\\ \lbrack D_t u &= -a\overline{u}^{t,\theta}\rbrack^{n+\theta},\\ \lbrack D_t u &= -\overline{au}^{t,\theta}\rbrack^{n+\theta} \thinspace . \end{align*} $$
$$ \begin{equation} u'(t) = -a(t)u(t) + b(t),\quad t\in (0,T],\quad u(0)=I \thinspace . \label{decay:problem:ab} \end{equation} $$
$$ \begin{align*} \lbrack D^+_t u &= -au + b\rbrack^n,\\ \lbrack D^-_t u &= -au + b\rbrack^n,\\ \lbrack D_t u &= -a\overline{u}^t + b\rbrack^{n+\frac{1}{2}},\\ \lbrack D_t u &= \overline{-au+b}^t\rbrack^{n+\frac{1}{2}},\\ \lbrack D_t u &= -a\overline{u}^{t,\theta} + b\rbrack^{n+\theta},\\ \lbrack D_t u &= \overline{-au+b}^{t,\theta}\rbrack^{n+\theta} \thinspace . \end{align*} $$
$$ \begin{equation} u^{n+1} = ((1 - \Delta t(1-\theta)a^n)u^n + \Delta t(\theta b^{n+1} + (1-\theta)b^n))(1 + \Delta t\theta a^{n+1})^{-1} \thinspace . \end{equation} $$
Implementation where \( a(t) \) and \( b(t) \) are given as Python functions (see file dc_vc.py):
def solver(I, a, b, T, dt, theta):
"""
Solve u'=-a(t)*u + b(t), u(0)=I,
for t in (0,T] with steps of dt.
a and b are Python functions of t.
"""
dt = float(dt) # avoid integer division
N = int(round(T/dt)) # no of time intervals
T = N*dt # adjust T to fit time step dt
u = zeros(N+1) # array of u[n] values
t = linspace(0, T, N+1) # time mesh
u[0] = I # assign initial condition
for n in range(0, N): # n=0,1,...,N-1
u[n+1] = ((1 - dt*(1-theta)*a(t[n]))*u[n] + \
dt*(theta*b(t[n+1]) + (1-theta)*b(t[n])))/\
(1 + dt*theta*a(t[n+1]))
return u, t
Plain functions:
def a(t):
return a_0 if t < tp else k*a_0
def b(t):
return 1
Better implementation: class with the parameters a0, tp, and k as attributes and a special method __call__ for evaluating \( a(t) \):
class A:
def __init__(self, a0=1, k=2):
self.a0, self.k = a0, k
def __call__(self, t):
return self.a0 if t < self.tp else self.k*self.a0
a = A(a0=2, k=1) # a behaves as a function a(t)
Use one-liner lambda function:
a = lambda t: a_0 if t < tp else k*a_0
In general,
f = lambda arg1, arg2, ...: expressin
is equivalent to
def f(arg1, arg2, ...):
return expression
One can use lambda functions directly in calls:
u, t = solver(1, lambda t: 1, lambda t: 1, T, dt, theta)
for a problem \( u'=-u+1 \), \( u(0)=1 \).
A lambda function can appear anywhere where a variable can appear.
import nose.tools as nt
def test_constant_solution():
"""
Test problem where u=u_const is the exact solution, to be
reproduced (to machine precision) by any relevant method.
"""
def exact_solution(t):
return u_const
def a(t):
return 2.5*(1+t**3) # can be arbitrary
def b(t):
return a(t)*u_const
u_const = 2.15
theta = 0.4; I = u_const; dt = 4
N = 4 # enough with a few steps
u, t = solver(I=I, a=a, b=b, T=N*dt, dt=dt, theta=theta)
print u
u_e = exact_solution(t)
difference = abs(u_e - u).max() # max deviation
nt.assert_almost_equal(difference, 0, places=14)
$$ \begin{align} \lbrack D_t^+ t\rbrack^n = \frac{t_{n+1}-t_n}{\Delta t}=\frac{(n+1)\Delta t - n\Delta t}{\Delta t}=1,\label{decay:fd2:Dop:tn:fw}\\ \lbrack D_t^- t\rbrack^n = \frac{t_{n}-t_{n-1}}{\Delta t}=\frac{n\Delta t - (n-1)\Delta t}{\Delta t}=1,\label{decay:fd2:Dop:tn:bw}\\ \lbrack D_t t\rbrack^n = \frac{t_{n+\frac{1}{2}}-t_{n-\frac{1}{2}}{\Delta t}=\frac{(n+\frac{1}{2})\Delta t - (n-\frac{1}{2})\Delta t}{\Delta t}=1\label{decay:fd2:Dop:tn:cn} \thinspace . \end{align} $$ The difference equation for the Forward Euler scheme
$$ [D_^+ u = -au + b]^n, $$ with \( a^n=a(t_n) \), \( b^n=c + a(t_n)(ct_n + I) \), and \( u^n=ct_n + I \) then results in
$$ c = -a(t_n)(ct_n+I) + c + a(t_n)(ct_n + I) = c $$
Verification function as a nose test:
def test_linear_solution():
"""
Test problem where u=c*t+I is the exact solution, to be
reproduced (to machine precision) by any relevant method.
"""
def exact_solution(t):
return c*t + I
def a(t):
return t**0.5 # can be arbitrary
def b(t):
return c + a(t)*exact_solution(t)
theta = 0.4; I = 0.1; dt = 0.1; c = -0.5
T = 4
N = int(T/dt) # no of steps
u, t = solver(I=I, a=a, b=b, T=N*dt, dt=dt, theta=theta)
u_e = exact_solution(t)
difference = abs(u_e - u).max() # max deviation
print difference
# No of decimal places for comparison depend on size of c
nt.assert_almost_equal(difference, 0, places=14)
Sample system:
$$ \begin{align} u' &= -a_u u + a_vv,\\ v' &= -a_vv + a_uu, \end{align} $$ for constants \( a_u, a_v>0 \).
The Forward Euler method:
$$ \begin{align} u^{n+1} &= u^n + \Delta t (-a_u u^n + a_vv^n),\\ v^{n+1} &= u^n + \Delta t (-a_vv^n + a_uu^n) \thinspace . \end{align} $$
The Backward Euler scheme: $$ \begin{align} u^{n+1} &= u^n + \Delta t (-a_u u^{n+1} + a_vv^{n+1}),\\ v^{n+1} &= v^n + \Delta t (-a_vv^{n+1} + a_uu^{n+1}), \thinspace . \end{align} $$ which is a \( 2\times 2 \) linear system: in $$ \begin{align} (1 + \Delta t a_u)u^{n+1} + a_vv^{n+1}) &= u^n ,\\ a_uu^{n+1} + (1 + \Delta t a_v) v^{n+1} &= v^n \thinspace . \end{align} $$
The standard form for ODEs: $$ \begin{equation} u' = f(u,t),\quad u(0)=I, \label{decay:ode:general} \end{equation} $$
\( u \) and \( f \): scalar or vector.
Vectors in case of ODE systems: $$ u(t) = (u^{(0)}(t),u^{(1)}(t),\ldots,u^{(m-1)}(t)) \thinspace . $$
$$ \begin{align*} f(u, t) = ( & f^{(0)}(u^{(0)}(t),\ldots,u^{(m-1)}(t)),\\ & f^{(1)}(u^{(0)}(t),\ldots,u^{(m-1)}(t)),\\ & \vdots\\ & f^{(m-1)}(u^{(0)}(t),\ldots,u^{(m-1)}(t))) \thinspace . \end{align*} $$
Higher-order ODEs are most often expressed as first-order systems.
Odespy features simple Python implementations of the most fundamental schemes as well as Python interfaces to several famous packages for solving ODEs: ODEPACK, Vode, rkc.f, rkf45.f, Radau5, as well as the ODE solvers in SciPy, SymPy, and odelab.
Typical usage:
def f(u, t):
return -a*u
import odespy
import numpy as np
I = 1; a = 2; T = 6; dt = 1
solver = odespy.RK4(f)
solver.set_initial_condition(I)
t_mesh = np.linspace(0, T, N+1)
u, t = solver.solve(t_mesh)
solvers = [odespy.RK2(f),
odespy.RK3(f),
odespy.RK4(f),
odespy.BackwardEuler(f, nonlinear_solver='Newton')]
for solver in solvers:
solver.set_initial_condition(I)
u, t = solver.solve(t)
# + lots of plot code...
Figure 6: Behavior of different schemes for the decay equation.
The 4-th order Runge-Kutta method (RK4) is the method of choice!
Figure 7: Choice of adaptive time mesh by the Dormand-Prince method for difference tolerances.
$$ u'(t_{n+1}) \approx \frac{3u^{n+1} - 4u^{n} + u^{n-1}}{2\Delta t},$$
Scheme: $$ u^{n+1} = \frac{4}{3}u^n - \frac{1}{3}u^{n-1} + \frac{2}{3}\Delta t f(u^{n+1}, t_{n+1}) \thinspace . \label{decay:fd2:bw:2step} $$
$$ \begin{equation} u'(t_n)\approx \frac{u^{n+1}-u^{n-1}}{2\Delta t} = [D_{2t} u]^n, \end{equation} $$
Scheme:
$$ [D_{2t} u = -a + b]^n,$$ or, if we write out and solve for the unknown \( u^{n+1} \), $$ \begin{equation} u^{n+1} = u^{n-1} + \Delta t f(u^n, t_n) \thinspace . \label{decay:fd2:leapfrog} \end{equation} $$
After computing \( u^{n+1} \), stabilize Leapfrog by $$ \begin{equation} u^n\ \leftarrow\ u^n + \gamma (u^{n-1} - 2u^n + u^{n+1}) \label{decay:fd2:leapfrog:filtered} \thinspace . \end{equation} $$
Forward-Euler + approximate Crank-Nicolson: $$ \begin{align} u^* &= u^n + \Delta t f(u^n, t_n), \label{decay:fd2:RK2:s1}\\ u^{n+1} &= u^n + \Delta t \frac{1}{2} \left( f(u^n, t_n) + f(u^*, t_{n+1}) \right), \label{decay:fd2:RK2:s2} \end{align} $$
$$ \begin{equation} u^{n+1} = u^n + \frac{1}{2}\Delta t\left( 3f(u^n, t_n) - f(u^{n-1}, t_{n-1}) \right) \thinspace . \label{decay:fd2:AB2} \end{equation} $$
$$ \begin{equation} u^{n+1} = u^n + \frac{1}{12}\left( 23f(u^n, t_n) - 16 f(u^{n-1},t_{n-1}) + 5f(u^{n-2}, t_{n-2})\right) \thinspace . \label{decay:fd2:AB3} \end{equation} $$