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.

Basic Usage

Overview

The typical usage of the tools consists of six steps. These are outlined in generic form below.

Step 1

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.

Step 2

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.

Step 3

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)

Step 4

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.

Step 5

Solve the ODE problem, which means to compute u(t) at some discrete user-specified time points t_i, for i=0,1,\ldots,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 $j$-th unknown function at the $i$-th time point 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.)

Step 6

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()

Example: Exponential Growth

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.