Solving ODE numerically by GSL


  • Joshua Kantor (2004-2006)

  • Robert Marik (2010 - fixed docstrings)

class sage.calculus.ode.PyFunctionWrapper

Bases: object

class sage.calculus.ode.ode_solver(function=None, jacobian=None, h=0.01, error_abs=1e-10, error_rel=1e-10, a=False, a_dydt=False, scale_abs=False, algorithm='rkf45', y_0=None, t_span=None, params=[])

Bases: object

ode_solver() is a class that wraps the GSL libraries ode solver routines To use it instantiate a class,:

sage: T=ode_solver()

To solve a system of the form dy_i/dt=f_i(t,y), you must supply a vector or tuple/list valued function f representing f_i. The functions f and the jacobian should have the form foo(t,y) or foo(t,y,params). params which is optional allows for your function to depend on one or a tuple of parameters. Note if you use it, params must be a tuple even if it only has one component. For example if you wanted to solve \(y''+y=0\). You need to write it as a first order system:

y_0' = y_1
y_1' = -y_0

In code:

sage: f = lambda t,y:[y[1],-y[0]]
sage: T.function=f

For some algorithms the jacobian must be supplied as well, the form of this should be a function return a list of lists of the form [ [df_1/dy_1,...,df_1/dy_n], ..., [df_n/dy_1,...,df_n,dy_n], [df_1/dt,...,df_n/dt] ].

There are examples below, if your jacobian was the function my_jacobian you would do:

sage: T.jacobian = my_jacobian     # not tested, since it doesn't make sense to test this

There are a variety of algorithms available for different types of systems. Possible algorithms are

  • rkf45 - runga-kutta-felhberg (4,5)

  • rk2 - embedded runga-kutta (2,3)

  • rk4 - 4th order classical runga-kutta

  • rk8pd - runga-kutta prince-dormand (8,9)

  • rk2imp - implicit 2nd order runga-kutta at gaussian points

  • rk4imp - implicit 4th order runga-kutta at gaussian points

  • bsimp - implicit burlisch-stoer (requires jacobian)

  • gear1 - M=1 implicit gear

  • gear2 - M=2 implicit gear

The default algorithm is rkf45. If you instead wanted to use bsimp you would do:

sage: T.algorithm="bsimp"

The user should supply initial conditions in y_0. For example if your initial conditions are y_0=1,y_1=1, do:

sage: T.y_0=[1,1]

The actual solver is invoked by the method ode_solve(). It has arguments t_span, y_0, num_points, params. y_0 must be supplied either as an argument or above by assignment. Params which are optional and only necessary if your system uses params can be supplied to ode_solve or by assignment.

t_span is the time interval on which to solve the ode. There are two ways to specify t_span:

  • If num_points is not specified then the sequence t_span is used as the time points for the solution. Note that the first element t_span[0] is the initial time, where the initial condition y_0 is the specified solution, and subsequent elements are the ones where the solution is computed.

  • If num_points is specified and t_span is a sequence with just 2 elements, then these are the starting and ending times, and the solution will be computed at num_points equally spaced points between t_span[0] and t_span[1]. The initial condition is also included in the output so that num_points+1 total points are returned. E.g. if t_span = [0.0, 1.0] and num_points = 10, then solution is returned at the 11 time points [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0].

(Note that if num_points is specified and t_span is not length 2 then t_span are used as the time points and num_points is ignored.)

Error is estimated via the expression D_i = error_abs*s_i+error_rel*(a|y_i|+a_dydt*h*|y_i'|). The user can specify error_abs (1e-10 by default), error_rel (1e-10 by default) a (1 by default), a_(dydt) (0 by default) and s_i (as scaling_abs which should be a tuple and is 1 in all components by default). If you specify one of a or a_dydt you must specify the other. You may specify a and a_dydt without scaling_abs (which will be taken =1 be default). h is the initial step size which is (1e-2) by default.

ode_solve solves the solution as a list of tuples of the form, [ (t_0,[y_1,...,y_n]),(t_1,[y_1,...,y_n]),...,(t_n,[y_1,...,y_n])].

This data is stored in the variable solutions:

sage: T.solution               # not tested


Consider solving the Van der Pol oscillator \(x''(t) + ux'(t)(x(t)^2-1)+x(t)=0\) between \(t=0\) and \(t= 100\). As a first order system it is \(x'=y\), \(y'=-x+uy(1-x^2)\). Let us take \(u=10\) and use initial conditions \((x,y)=(1,0)\) and use the runga-kutta prince-dormand algorithm.

sage: def f_1(t,y,params):
....:    return[y[1],-y[0]-params[0]*y[1]*(y[0]**2-1.0)]

sage: def j_1(t,y,params):
....:    return [ [0.0, 1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-params[0]*(y[0]*y[0]-1.0)], [0.0, 0.0] ]

sage: T=ode_solver()
sage: T.algorithm="rk8pd"
sage: T.function=f_1
sage: T.jacobian=j_1
sage: T.ode_solve(y_0=[1,0],t_span=[0,100],params=[10.0],num_points=1000)
sage: outfile = os.path.join(SAGE_TMP, 'sage.png')
sage: T.plot_solution(filename=outfile)

The solver line is equivalent to:

sage: T.ode_solve(y_0=[1,0],t_span=[x/10.0 for x in range(1000)],params = [10.0])

Let’s try a system:


We will not use the jacobian this time and will change the error tolerances.

sage: g_1= lambda t,y: [y[1]*y[2],-y[0]*y[2],-0.51*y[0]*y[1]]
sage: T.function=g_1
sage: T.y_0=[0,1,1]
sage: T.scale_abs=[1e-4,1e-4,1e-5]
sage: T.error_rel=1e-4
sage: T.ode_solve(t_span=[0,12],num_points=100)

By default T.plot_solution() plots the y_0, to plot general y_i use:

sage: T.plot_solution(i=0, filename=outfile)
sage: T.plot_solution(i=1, filename=outfile)
sage: T.plot_solution(i=2, filename=outfile)

The method interpolate_solution will return a spline interpolation through the points found by the solver. By default y_0 is interpolated. You can interpolate y_i through the keyword argument i.

sage: f = T.interpolate_solution()
sage: plot(f,0,12).show()
sage: f = T.interpolate_solution(i=1)
sage: plot(f,0,12).show()
sage: f = T.interpolate_solution(i=2)
sage: plot(f,0,12).show()
sage: f = T.interpolate_solution()
sage: f(pi)

The solver attributes may also be set up using arguments to ode_solver. The previous example can be rewritten as:

sage: T = ode_solver(g_1,y_0=[0,1,1],scale_abs=[1e-4,1e-4,1e-5],error_rel=1e-4, algorithm="rk8pd")
sage: T.ode_solve(t_span=[0,12],num_points=100)
sage: f = T.interpolate_solution()
sage: f(pi)

Unfortunately because Python functions are used, this solver is slow on systems that require many function evaluations. It is possible to pass a compiled function by deriving from the class ode_sysem and overloading c_f and c_j with C functions that specify the system. The following will work in the notebook:

cimport sage.calculus.ode
import sage.calculus.ode
from sage.libs.gsl.all cimport *

cdef class van_der_pol(sage.calculus.ode.ode_system):
    cdef int c_f(self,double t, double *y,double *dydt):
        return GSL_SUCCESS
    cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt):
        return GSL_SUCCESS

After executing the above block of code you can do the following (WARNING: the following is not automatically doctested):

sage: T = ode_solver()                     # not tested
sage: T.algorithm = "bsimp"                # not tested
sage: vander = van_der_pol()               # not tested
sage: T.function=vander                    # not tested
sage: T.ode_solve(y_0 = [1,0], t_span=[0,2000], num_points=1000)   # not tested
sage: T.plot_solution(i=0, filename=os.path.join(SAGE_TMP, 'test.png'))        # not tested
ode_solve(t_span=False, y_0=False, num_points=False, params=[])
plot_solution(i=0, filename=None, interpolate=False, **kwds)

Plot a one dimensional projection of the solution.


  • i – (non-negative integer) composant of the projection

  • filename – (string or None) whether to plot the picture or save it in a file

  • interpolate – whether to interpolate between the points of the discretized solution

  • additional keywords are passed to the graphics primitive


sage: T = ode_solver()
sage: T.function = lambda t,y: [cos(y[0]) * sin(t)]
sage: T.jacobian = lambda t,y: [[-sin(y[0]) * sin(t)]]
sage: T.ode_solve(y_0=[1],t_span=[0,20],num_points=1000)
sage: T.plot_solution()

And with some options:

sage: T.plot_solution(color='red', axes_labels=["t", "x(t)"])
class sage.calculus.ode.ode_system

Bases: object