Solving ODE numerically by GSL#
AUTHORS:
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=[])[source]#
Bases:
object
ode_solver()
is a class that wraps the GSL library’s ode solver routines.To use it, instantiate the class:
sage: T = ode_solver()
>>> from sage.all import * >>> 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 formfoo(t,y)
orfoo(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 would 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
>>> from sage.all import * >>> f = lambda t, y: [y[Integer(1)], -y[Integer(0)]] >>> T.function = f
For some algorithms, the jacobian must be supplied as well, the form of this should be a function returning 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
>>> from sage.all import * >>> 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'
– Runge-Kutta-Fehlberg (4,5)'rk2'
– embedded Runge-Kutta (2,3)'rk4'
– 4th order classical Runge-Kutta'rk8pd'
– Runge-Kutta Prince-Dormand (8,9)'rk2imp'
– implicit 2nd order Runge-Kutta at gaussian points'rk4imp'
– implicit 4th order Runge-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"
>>> from sage.all import * >>> 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]
>>> from sage.all import * >>> T.y_0 = [Integer(1),Integer(1)]
The actual solver is invoked by the method
ode_solve()
. It has argumentst_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 usesparams
can be supplied toode_solve
or by assignment.t_span
is the time interval on which to solve the ode. There are two ways to specifyt_span
:If
num_points
is not specified, then the sequencet_span
is used as the time points for the solution. Note that the first elementt_span[0]
is the initial time, where the initial conditiony_0
is the specified solution, and subsequent elements are the ones where the solution is computed.If
num_points
is specified andt_span
is a sequence with just 2 elements, then these are the starting and ending times, and the solution will be computed atnum_points
equally spaced points betweent_span[0]
andt_span[1]
. The initial condition is also included in the output so thatnum_points + 1
total points are returned. E.g. ift_span = [0.0, 1.0]
andnum_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 andt_span
is not length 2 thent_span
are used as the time points andnum_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 specifyerror_abs
(1e-10 by default),error_rel
(1e-10 by default),a
(1 by default),a_dydt
(0 by default) ands_i
(asscaling_abs
which should be a tuple and is 1 in all components by default).
If you specify one of
a
ora_dydt
you must specify the other. You may specifya
anda_dydt
withoutscaling_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
>>> from sage.all import * >>> T.solution # not tested
EXAMPLES:
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 Runge-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: import tempfile sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot ....: T.plot_solution(filename=f.name)
>>> from sage.all import * >>> def f_1(t, y, params): ... return [y[Integer(1)], -y[Integer(0)] - params[Integer(0)]*y[Integer(1)]*(y[Integer(0)]**Integer(2)-RealNumber('1.0'))] >>> def j_1(t, y, params): ... return [[RealNumber('0.0'), RealNumber('1.0')], ... [-RealNumber('2.0')*params[Integer(0)]*y[Integer(0)]*y[Integer(1)] - RealNumber('1.0'), -params[Integer(0)]*(y[Integer(0)]*y[Integer(0)]-RealNumber('1.0'))], ... [RealNumber('0.0'), RealNumber('0.0')]] >>> T = ode_solver() >>> T.algorithm = "rk8pd" >>> T.function = f_1 >>> T.jacobian = j_1 >>> T.ode_solve(y_0=[Integer(1),Integer(0)], t_span=[Integer(0),Integer(100)], params=[RealNumber('10.0')], num_points=Integer(1000)) >>> import tempfile >>> with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot ... T.plot_solution(filename=f.name)
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])
>>> from sage.all import * >>> T.ode_solve(y_0=[Integer(1),Integer(0)], t_span=[x/RealNumber('10.0') for x in range(Integer(1000))], params=[RealNumber('10.0')])
Let’s try a system:
y_0'=y_1*y_2 y_1'=-y_0*y_2 y_2'=-.51*y_0*y_1
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)
>>> from sage.all import * >>> g_1 = lambda t,y: [y[Integer(1)]*y[Integer(2)], -y[Integer(0)]*y[Integer(2)], -RealNumber('0.51')*y[Integer(0)]*y[Integer(1)]] >>> T.function = g_1 >>> T.y_0 = [Integer(0),Integer(1),Integer(1)] >>> T.scale_abs = [RealNumber('1e-4'), RealNumber('1e-4'), RealNumber('1e-5')] >>> T.error_rel = RealNumber('1e-4') >>> T.ode_solve(t_span=[Integer(0),Integer(12)], num_points=Integer(100))
By default
T.plot_solution()
plots the \(y_0\); to plot general \(y_i\), use:sage: with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot ....: T.plot_solution(i=0, filename=f.name) ....: T.plot_solution(i=1, filename=f.name) ....: T.plot_solution(i=2, filename=f.name)
>>> from sage.all import * >>> with tempfile.NamedTemporaryFile(suffix=".png") as f: # needs sage.plot ... T.plot_solution(i=Integer(0), filename=f.name) ... T.plot_solution(i=Integer(1), filename=f.name) ... T.plot_solution(i=Integer(2), filename=f.name)
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() # needs sage.plot sage: f = T.interpolate_solution(i=1) sage: plot(f,0,12).show() # needs sage.plot sage: f = T.interpolate_solution(i=2) sage: plot(f,0,12).show() # needs sage.plot sage: f = T.interpolate_solution() sage: from math import pi sage: f(pi) 0.5379...
>>> from sage.all import * >>> f = T.interpolate_solution() >>> plot(f,Integer(0),Integer(12)).show() # needs sage.plot >>> f = T.interpolate_solution(i=Integer(1)) >>> plot(f,Integer(0),Integer(12)).show() # needs sage.plot >>> f = T.interpolate_solution(i=Integer(2)) >>> plot(f,Integer(0),Integer(12)).show() # needs sage.plot >>> f = T.interpolate_solution() >>> from math import pi >>> f(pi) 0.5379...
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) 0.5379...
>>> from sage.all import * >>> T = ode_solver(g_1, y_0=[Integer(0),Integer(1),Integer(1)], scale_abs=[RealNumber('1e-4'),RealNumber('1e-4'),RealNumber('1e-5')], ... error_rel=RealNumber('1e-4'), algorithm="rk8pd") >>> T.ode_solve(t_span=[Integer(0),Integer(12)], num_points=Integer(100)) >>> f = T.interpolate_solution() >>> f(pi) 0.5379...
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_system
and overloadingc_f
andc_j
with C functions that specify the system. The following will work in the notebook:%cython 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): dydt[0]=y[1] dydt[1]=-y[0]-1000*y[1]*(y[0]*y[0]-1) return GSL_SUCCESS cdef int c_j(self, double t,double *y,double *dfdy,double *dfdt): dfdy[0]=0 dfdy[1]=1.0 dfdy[2]=-2.0*1000*y[0]*y[1]-1.0 dfdy[3]=-1000*(y[0]*y[0]-1.0) dfdt[0]=0 dfdt[1]=0 return GSL_SUCCESS
After executing the above block of code you can do the following (WARNING: the following is not automatically doctested):
sage: # not tested sage: T = ode_solver() sage: T.algorithm = "bsimp" sage: vander = van_der_pol() sage: T.function = vander sage: T.ode_solve(y_0=[1, 0], t_span=[0, 2000], ....: num_points=1000) sage: from tempfile import NamedTemporaryFile sage: with NamedTemporaryFile(suffix=".png") as f: ....: T.plot_solution(i=0, filename=f.name)
>>> from sage.all import * >>> # not tested >>> T = ode_solver() >>> T.algorithm = "bsimp" >>> vander = van_der_pol() >>> T.function = vander >>> T.ode_solve(y_0=[Integer(1), Integer(0)], t_span=[Integer(0), Integer(2000)], ... num_points=Integer(1000)) >>> from tempfile import NamedTemporaryFile >>> with NamedTemporaryFile(suffix=".png") as f: ... T.plot_solution(i=Integer(0), filename=f.name)
- plot_solution(i=0, filename=None, interpolate=False, **kwds)[source]#
Plot a one dimensional projection of the solution.
INPUT:
i
– (non-negative integer) composant of the projectionfilename
– (string orNone
) whether to plot the picture or save it in a fileinterpolate
– whether to interpolate between the points of the discretized solutionadditional keywords are passed to the graphics primitive
EXAMPLES:
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() # needs sage.plot
>>> from sage.all import * >>> T = ode_solver() >>> T.function = lambda t,y: [cos(y[Integer(0)]) * sin(t)] >>> T.jacobian = lambda t,y: [[-sin(y[Integer(0)]) * sin(t)]] >>> T.ode_solve(y_0=[Integer(1)],t_span=[Integer(0),Integer(20)],num_points=Integer(1000)) >>> T.plot_solution() # needs sage.plot
And with some options:
sage: T.plot_solution(color='red', axes_labels=["t", "x(t)"]) # needs sage.plot
>>> from sage.all import * >>> T.plot_solution(color='red', axes_labels=["t", "x(t)"]) # needs sage.plot
- class sage.calculus.ode.ode_system#
Bases:
object