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=[])¶ 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 functionf
representingf_i
. The functionsf
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 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-kuttark8pd
- runga-kutta prince-dormand (8,9)rk2imp
- implicit 2nd order runga-kutta at gaussian pointsrk4imp
- implicit 4th order runga-kutta at gaussian pointsbsimp
- implicit burlisch-stoer (requires jacobian)gear1
- M=1 implicit geargear2
- M=2 implicit gear
The default algorithm is
rkf45
. If you instead wanted to usebsimp
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 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 uses params 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
(as scaling_abs which should be a tuple and is 1 in all components by default). If you specify one ofa
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
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 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:
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)
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) 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...
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 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: 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
-
interpolate_solution
(i=0)¶
-
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.
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()
And with some options:
sage: T.plot_solution(color='red', axes_labels=["t", "x(t)"])
-
class
sage.calculus.ode.
ode_system
¶ Bases:
object