Tutorial
The odesolvers package contains tools for solving ordinary
differential equations (ODEs). The user specifies the problem through
high-level Python code. Both scalar ODEs and systems of ODEs are
supported. A wide range of numerical methods for ODEs are offered.
Tutorial: html/index.html
The typical usage of the tools consists of six steps. These are outlined in generic form below.
Write the ODE problem on the generic form u' = f(u, t), where u(t) is the unknown function to be solved for, or a vector of unknown functions of time in case of a system of ODEs.
Implement the right-hand side function f(u, t), which
defines the ODEs to be solved, as a Python function f(u, t). The
argument u is either a float object (in case of a scalar ODE) or a
numpy array object (in case of a system of ODEs). Alternatively, a
class with a __call__ method can be used instead of a plain
function. Some tools in this package also allow implementation of f
in Fortran or C for increased efficiency.
Create a method object:
method = classname(f)
where classname is the name of a class in this package implementing
the desired numerical method.
Many solver classes has a range of parameters that the user can set to
control various parts of the solution process. The parameters are
documented in the doc string of the class (pydoc classname will list
the documentation in a terminal window). One can either specify parameters
at construction time, via extra keyword arguments to the constructor:
method = classname(f, prm1=value1, prm2=value2, ...)
or at any time using the set method:
method.set(prm1=value1, prm2=value2, prm3=value3) ... method.set(prm4=value4)
Set the initial condition, u(0)=U_0:
method.set_initial_condition(U0)
where U0 is either a number, for a scalar ODE, or a sequence (list, tuple,
numpy array), for a system of ODEs.
Solve the ODE problem, which means to compute u(t) at some discrete user-specified time points t_i, for i=0,1,...,n:
T = ... # end time time_points = numpy.linspace(0, T, n+1) u, t = method.solve(time_points)
In case of a scalar ODE, the returned solution u is a one-dimensional
numpy array where u[i] holds the solution at time point t[i].
For a system of ODEs, the returned u is a two-dimensional numpy
array where u[i,j] holds the solution of the t[i] (u_j(t_i) in mathematics
notation).
The time_points array specifies the time points where we want the
solution to be computed. The returned array t is the same as
time_points. The simplest numerical methods in the odesolvers
package apply the time_points array directly in the solution. That
is, the time steps used are given by time_points[i] -
time_points[i-1], for i=0,1,...,len(time_points)-1. Some more
advanced adaptive methods compute the time steps internally. Then the
time_points array is just a specification of the time points where
we want to know the solution.
(hpl: Need to document how the tp array is used.)
(hpl: Need to comment on storage.)
The solve method in solver classes also allows a second argument,
terminate, which is a user-implemented Python function specifying
when the solution process is to be terminated. For example,
terminating when the solution reaches an asymptotic (known) value
a can be done by:
def terminate(u, t, step_no):
# u and t are arrays. Most recent solution is u[step_no].
tolerance = 1E-6
return True if abs(u[step_no] - a) < tolerance else False
u, t = method.solve(time_points, terminate)
The arguments transferred to the terminate function are the
solution array u, the corresponding time points t, and
an integer step_no reflecting the most recently computed u
value. That is, u[step_no] is most recently computed value of u.
(The array data u[step_no+1:] will typically be zero as these
are uncomputed future values.)
Extract solution components for plotting and further analysis.
Since the u array returned from method.solve stores all unknown
functions at all discrete time levels, one usually wants to extract
individual unknowns as one-dimensional arrays. Here is an example
where unknown 0 and k are extracted in individual arrays and plotted:
u_0 = u[:,0] u_k = u[:,k] from matplotlib.pyplot import plot, show plot(t, u_0, t, u_k) show()
Our first example concerns the simple scalar ODE problem u'=cu, u(0)=A, where A>0 and c>0 are known constants. Using a standard Runge-Kutta method of order four, the code for solving the problem in the time interval [0,10], with N=30 time steps, looks like this:
c = 0.1
A = 1.5
def f(u, t):
return c*u
import odesolvers
method = odesolvers.RK4(f)
method.set_initial_condition(A)
import numpy
N = 30 # no of time steps
time_points = numpy.linspace(0, 10, N+1)
u, t = method.solve(time_points)
from matplotlib.pyplot import *
plot(t, u)
show()
With the RK4 method and other non-adaptive methods
the time steps are dictated by the time_points array.
A constant time step of size time_points[1] - time_points[0] = 1/3
is implied in the present example.
The right-hand side function and all physical parameters are often lumped together in a class, for instance:
class ExponentialGrowth:
def __init__(self, c=1, A=1):
self.c, self.A = c, A
def __call__(self, u, t):
return self.c*u
f = ExponentialGrowth(c=0.1, A=1.5)
import odesolvers
method = odesolvers.RK4(f)
method.set_initial_condition(f.A)
We may also compare the numerical and exact solution:
u_exact = f.A*numpy.exp(f.c*t) error = numpy.abs(u_exact - u).max() print 'Max deviation of numerical solution:', error from matplotlib.pyplot import * plot(t, u, 'r-', t, u_exact, 'bo') legend(['RK4, N=%d' % N, 'exact']) show()
Instead of having the c variable as a global variable or
in a class, we may include it as an extra argument to f, either
as a positional argument or as a keyword argument. Positional arguments
can be sent to f via the constructor argument f_args (a list/tuple of
variables), while a dictionary f_kwargs is used to transfer
keyword arguments to f via the constructor:
# f has extra positional argument
def f(u, t, c):
return c*u
method = odesolvers.RK4(f, f_args=[c])
# Alternative: f has extra keyword argument
def f(u, t, c=1):
return c*u
method = odesolvers.RK4(f, f_kwargs={'c': c})
The right-hand side function f may feature both extra positional
arguments and extra keyword arguments:
def f(u, t, arg1, arg2, arg3, ..., kwarg1=val1, kwarg2=val2, ...):
...
method = odesolvers.classname(f,
f_args=[arg1, arg2, arg3, ...],
f_kwargs=dict(kwarg1=val1, kwarg2=val2, ...))
# Alternative setting of f_args and f_kwargs
method.set(f_args=[arg1, arg2, arg3, ...],
f_kwargs=dict(kwarg1=val1, kwarg2=val2, ...))
Solvers will call f as f(u, t, *f_args, **f_kwargs).
It is easy to simulate for some time interval [0, T_1], the continue with u(T_1) as new initial condition and simulate for t in [T_1, T_2] and so on. Let us divide the time domain [0,10] into subdomains and compute the solution for each subdomain in sequence. The following program performs the steps:
c = 0.1
A = 1.5
def f(u, t):
return c*u
import odesolvers, numpy
from matplotlib.pyplot import plot, hold, show
method = odesolvers.RK4(f)
# Split time domain into subdomains and
# integrate the ODE in each subdomain
T = [0, 1, 3, 6, 10, 20]
N_tot = 30 # no of time intervals in total
dt = float(T[-1])/N_tot # time step, kept fixed
u = []; t = [] # collectors for u and t
for i in range(len(T)-1):
T_interval = T[i+1] - T[i]
N = int(round(T_interval/dt))
time_points = numpy.linspace(T[i], T[i+1], N+1)
method.set_initial_condition(A) # at time_points[0]
print 'Solving in [%s, %s] with %d points' % \
(T[i], T[i+1], N)
ui, ti = method.solve(time_points)
A = ui[-1] # newest u is next initial condition
plot(ti, ui)
hold('on')
u.append(ui); t.append(ti)
# Can concatenate all the elements of u and t, if desired
u = numpy.concatenate(u); t = numpy.concatenate(t)
#plot(t, u, 'bo') # same curve
show()
The second example also concerns the ODE u'=cu, but this time with
c<0. We want to integrate until the asymptotic value u=0 is
reached with an accuracy of 0.001. For this purpose we need a
terminate function, which returns True when the termination
criterion is met:
tol = 0.001 # tolerance for terminating the simulation
def terminate(u, t, step_no):
# Most recent solution is in u[step_no] at time t[step_no]
if abs(u[step_no]) <= tol:
return True
else:
return False
The complete program may take the form:
c = -0.1
A = 1.5
def f(u, t):
return c*u
import odesolvers
method = odesolvers.RK4(f)
method.set_initial_condition(A)
import numpy
# Make sure integration interval [0, T] is large enough
N = 50
T = 150
time_points = numpy.linspace(0, T, N+1)
tol = 0.001 # tolerance for terminating the simulation
def terminate(u, t, step_no):
# Most recent solution is in u[step_no] at time t[step_no]
if abs(u[step_no]) <= tol:
return True
else:
return False
u, t = method.solve(time_points, terminate)
print "Solve u'=%g*u, u(0)=%g, for t in [%g, %g] and u>%g" % \
(c, A, time_points[0], time_points[-1], tol)
print 'Final u(t=%g)=%g after %d steps' % (t[-1], u[-1], len(u)-1)
from matplotlib.pyplot import *
print plot
plot(t, u)
show()
Note that the size of the returned arrays u and t fits
the time interval up to the point of termination by terminate.
The odesolvers package applies u for the unknown function or
vector of unknown functions and t as the name of the independent
variable. Many problems involve other symbols for functions and
independent variables. These symbols should be reflected in the code.
For example, here is a coding example involving the differential
equation y'(x)=-y(x), y(0)=1:
def myrhs(y, x):
return -y
import odesolvers
method = odesolvers.RK4(myrhs)
method.set_initial_condition(1)
import numpy
# Make sure integration interval [0, L] is large enough
N = 50
L = 10
x_points = numpy.linspace(0, L, N+1)
def terminate(y, x, stepnumber):
tol = 0.001
return True if abs(y[stepnumber]) < tol else False
y, x = method.solve(x_points, terminate)
print 'Final y(x=%g)=%g' % (x[-1], y[-1])
from matplotlib.pyplot import *
plot(x, y)
show()
As shown, we use y for u, x for t, and x_points instead
of time_points.
The angle theta of a pendulum with mass m and length L is governed by the equation (neglecting air resistance for simplicity) .. pure ASCII
theta'' + c*sin(theta) = 0, with theta(0) = Theta and theta'(0) = 0.
This problem must be expressed as a first-order ODE system if it is
going to be solved by the tools in the odesolvers package.
Introducing omega as theta', we get the first equation as
theta' = omega and the second as omega' = -c*sin(theta).
The initial conditions are theta(0) = Theta and omega'(0) = 0.
Now the f function must return a list or array with the two
right-hand side functions:
def f(u, t):
theta, omega = u
return [omega, -c*sin(theta)]
Note that when we have a system of ODEs with n components, the u
object sent to the f function is an array of length n,
representing the value of all components in the ODE system at time t.
Here we extract the two components of u in separate local variables
with names equal to what is used in the mathematical description of
the current problem.
The initial conditions must be specified as a list:
method = odesolvers.Heun(f) method.set_initial_condition([Theta, 0])
To specify the time points we here first decide on a number of periods (oscillations back and forth) to simulate and then on the time resolution of each period. (When Theta is small we can replace sin(theta) by theta and find an analytical solution theta(t) = Theta*cos(sqrt(c)*t) whose period is 2*pi/sqrt(c).):
freq = sqrt(c) # frequency of oscillations when Theta is small period = 2*pi/freq # the period of the oscillations T = 8*period # final time N_per_period = 20 # resolution of one period N = N_per_period*period time_points = numpy.linspace(0, T, N+1) u, t = method.solve(time_points)
The u returned from method.solve is a two-dimensional array, where the
columns hold the various solution functions of the ODE system. That is,
the first column holds theta and the second column holds
omega. For convenience we extract the individual solution
components in individual arrays:
theta = u[:,0] omega = u[:,1] from matplotlib.pyplot import * theta_small = Theta*numpy.cos(sqrt(c)*t) plot(t, theta, 'r-', t, theta_small, 'y-') legend(['theta', 'small angle approximation']) show()
We shall now make a more advanced solver by extending the previous example. More specifically, we shall
- represent the right-hand side function as class,
- compare several different solvers,
- compute error of numerical solutions.
Since we want to compare numerical errors in the various solvers we need a test problem where the exact solution is known. Approximating sin(theta) by theta (valid for small theta), gives the ODE system theta' = omega and omega' = -c*theta with initial conditions theta(0) = Theta and omega'(0) = 0.
Right-hand side functions with parameters can be handled by
including extra arguments via the f_args and f_kwargs functionality,
or by using a class where the parameters are attributes and
a __call__ method makes the class instance callable as a function. The section Parameters in the Right-Hand Side Function exemplifies the details.
A minimal class representation of the right-hand side
function in the present case looks like this:
class Problem:
def __init__(self, c, Theta):
self.c, self.Theta = float(c), float(Theta)
def __call__(self, u, t):
theta, omega = u; c = self.c
return [omega, -c*theta]
problem = Problem(c=1, Theta=pi/4)
It would be convenient to add an attribute period which holds
an estimate of the period of oscillations as we need this for
deciding on the complete time interval for solving the differential
equations. An appropriate extension of class Problem is therefore:
class Problem:
def __init__(self, c, Theta):
self.c, self.Theta = float(c), float(Theta)
self.freq = sqrt(c)
self.period = 2*pi/self.freq
def __call__(self, u, t):
theta, omega = u; c = self.c
return [omega, -c*theta]
problem = Problem(c=1, Theta=pi/4)
The second extension is to loop over many solvers. All solvers can be listed by:
>>> import odesolvers, pprint >>> solvers = list_all_solvers() >>> for solver in solvers: ... print solver ... AdamsBashMoulton2 AdamsBashMoulton3 AdamsBashforth2 AdamsBashforth3 AdamsBashforth4 AdaptiveResidual Backward2Step BackwardEuler BogackiShampine CashKarp Dop853 Dopri5 DormandPrince Fehlberg Euler ForwardEuler Heun Leapfrog LeapfrogFiltered Lsoda Lsodar Lsode Lsodes Lsodi Lsodis Lsoibt MidpointImplicit MidpointIter MyRungeKutta MySolver RKC RKF45 RK2 RungeKutta2 RK3 RungeKutta3 RK4 RungeKutta4 RKFehlberg SymPy_odefun ThetaRule Trapezoidal Vode
A similar function, list_available_solvers, returns a list of the
names of the solvers that are available in the current installation
(e.g., the Vode solver is only available if the comprehensive
scipy package is installed).
This is the list that is usually most relevant.
For now we explicitly choose a subset of the commonly available solvers:
import odesolvers
solvers = [
odesolvers.ThetaRule(problem, theta=0), # Forward Euler
odesolvers.ThetaRule(problem, theta=0.5), # Midpoint
odesolvers.ThetaRule(problem, theta=1), # Backward Euler
odesolvers.RK4(problem),
odesolvers.MidpointIter(problem, max_iter=2, eps_iter=0.01),
odesolvers.LeapfrogFiltered(problem),
]
It will be evident that the ThetaRule solver with theta=0 and
theta=1 (Forward and Backward Euler methods) gives growing and
decaying amplitudes, respectively, while the other solvers are
capable of reproducing the constant amplitude of the oscillations of
in the current mathematical model.
The loop over the chosen solvers may look like:
N_per_period = 20
T = 3*problem.period # final time
import numpy
import matplotlib.pyplot as mpl
legends = []
for method in solvers:
method_name = str(method)
print method_name
method.set_initial_condition([problem.Theta, 0])
N = N_per_period*problem.period
time_points = numpy.linspace(0, T, N+1)
u, t = method.solve(time_points)
theta = u[:,0]
legends.append(method_name)
mpl.plot(t, theta)
mpl.hold('on')
mpl.legend(legends)
mpl.savefig(__file__[:-3] + '.png')
mpl.show()
Comparison of solutions
We can extend this program to compute the error in each numerical
solution for different time step sizes.
Let results be a dictionary with the method name as
key, containing two sub dictionaries dt and error, which hold
a sequence of time steps and a sequence of corresponding
errors, respectively. The errors are computed by subtracting
the numerical solution from the exact solution:
theta_exact = lambda t: problem.Theta*numpy.cos(sqrt(problem.c)*t) u, t = method.solve(time_points) theta = u[:,0] error = numpy.abs(theta_exact(t) - theta)
The so-called L2 norm of the error array is a suitable
scalar error measure (square root of total error squared and integrated,
here numerically):
error_L2 = sqrt(numpy.sum(error**2)/dt)
where dt is the time step size.
Typical loops over solvers and resolutions look as follows:
T = num_periods*problem.period # final time
results = {}
resolutions = [10, 20, 40, 80, 160] # intervals per period
import numpy
for method in solvers:
method_name = str(method)
results[method_name] = {'dt': [], 'error': []}
method.set_initial_condition([problem.Theta, 0])
for N_per_period in resolutions:
N = N_per_period*problem.period
time_points = numpy.linspace(0, T, N+1)
u, t = method.solve(time_points)
theta = u[:,0]
error = numpy.abs(theta_exact(t) - theta)
error_L2 = sqrt(numpy.sum(error**2)/N)
if not numpy.isnan(error_L2): # drop nan
results[method_name]['dt'].append(t[1] - t[0])
results[method_name]['error'].append(error_L2)
Assuming the error to be of the form CDelta t^r, we can estimate
C and r from two consequtive experiments to obtain a sequence
of r values which (hopefully) convergences to a value that we can
view as the empirical convergence rate of a method.
Given the sequence of time steps and errors, a function in
the scitools package automatically computes the sequence of
r values:
# Analyze convergence
from scitools.convergencerate import OneDiscretizationPrm
pairwise_rates = OneDiscretizationPrm.pairwise_rates # short form
print '\n\nConvergence results for %d periods' % num_periods
for method_name in results:
rates, C = pairwise_rates(results[method_name]['dt'],
results[method_name]['error'])
rates = ', '.join(['%.1f' % rate for rate in rates])
print '%-20s r: %s E_min=%.1E' % \
(method_name, rates, min(results[method_name]['error']))
With 4 periods we get:
ThetaRule(theta=0) r: 2.9, 1.9, 1.4, 1.2 E_min=1.1E-01 RK2 r: 2.1, 2.0, 2.0, 2.0 E_min=8.4E-04 ThetaRule(theta=1) r: 0.3, 0.5, 0.7, 0.8 E_min=9.0E-02 RK4 r: 4.0, 4.0, 4.0, 4.0 E_min=2.6E-08 ThetaRule r: 2.0, 2.0, 2.0, 2.0 E_min=4.2E-04 Leapfrog r: 2.1, 2.0, 2.0, 2.0 E_min=8.5E-04 LeapfrogFiltered r: 0.2, 0.4, 0.6, 0.8 E_min=1.3E-01
The rates of the Forward and Backward Euler methods (1st and 3rd line) have
not yet converged to unity, as expected, while the 2nd-order
Runge-Kutta method, Leapfrog, and the ThetaRule with default value of theta) shows the expected
r=2 value. The 4th-order Runge-Kutta holds the promise of being of 4th
order, while the filtered Leapfrog method has slow convergence and
a fairly large error, which is also evident in the previous figure.
Extending the time domain to 20 periods makes many of the simplest methods inaccurate:
ThetaRule(theta=0) r: 10.4, 22.0, 18.2, 10.4 E_min=3.3E+02 RK2 r: 51.4, 16.1, 2.3, 2.1 E_min=1.1E-01 ThetaRule(theta=1) r: 11.2, 0.0, 0.1 E_min=5.0E-01 RK4 r: 1.0, 3.6, 4.0, 4.0 E_min=8.2E-05 ThetaRule r: 87.8, 0.2, 1.7, 2.0 E_min=5.3E-02 Leapfrog r: 95.9, 18.0, 1.0, 2.0 E_min=1.1E-01 LeapfrogFiltered r: -0.0, 121.2, 0.3, 0.1 E_min=5.2E-01