In response to Williams sage-2.0 plan I wanted to describe what I had done
with using gsl to implement a numerical ode solver. I believe that the
patch containing this will be applied after
doing a recent pull or upgrade but I'm not sure(is this true?). If not I
can send patches for people to play with. I have included the
documentation at the end so people can look at the syntax.
To see the documentation and examples
sage: ode_solver?
To be done
1. Testing/feedback: Is it easy to use compared to matlab,mathematica.
Any bugs or ways to crash it (unexpected input or user specified
functions that do weird things)?
2. More examples, doctests, improved documentation (I noticed a bunch of
typos just copying the documentation below)
3. Obvious additions people think would be useful. Currently it has a
plotting routine and a routine to produce an interpolated function from
the data points.
Ideas for extension:
1. It would be nice if there was some facility for automatically
converting a nth order ODE
to a system of first order ones.
2. It would be nice if there was some facility for automatically computing
the jacobian when the functions involved are rational functions and
elementary functions.
3. Numerically computing the jacobian: For the algorithms that require the
jacobian It would be possible to numerically compute the jacobian,
however I was wary of doing this by default. Does anyone have any knowledge
about
the benefits of this, can it cause instability (using the numerical
jacobian
instead of the exact one).
Accuracy testing:
1. I believe there are standard batches of tests for ode solvers.
it would be interesting if some of these could be applied to compare the
accuracy with matlab for example. This ode solver will be slower than
matlab's on systems that require many function evaluations because we are
forced to use python functions and there is significant overhead in
calling these. Still comparisons would be interesting. There is a way to
use compiled functions describe in the documentation but requires
writing the functions in C and so isn't really suitable for
general users.
____________________________________________________________________
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
y1_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
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 if 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, would do
sage: T.y_0=[1,0]
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 is 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.
If t_span is a tuple with just two time
values then the user must specify num_points,
and the system will be evaluated at num_points equally spaced
points between t_span[0]
and t_span[1]. If t_span is a tuple with more than two values
than the
values of y_i at points in time listed in t_span will be
returned.
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
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)]
sage: def j_1(t,y,params):
... return [ [0,
1.0],[-2.0*params[0]*y[0]*y[1]-1.0,-params[0]*(y[0]*y[0]-1.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],num_points=1000)
sage: T.plot_solution()
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])
Lets try a system
y_0'=y_2*y_3
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)
sage: T.plot_solution(i=1)
sage: T.plot_solution(i=2)
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)
Unfortunately because python functions are used, this solver is
slow on system 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 notebook
%pyrex
cimport sage.gsl.ode
import sage.gsl.ode
include 'gsl.pxi'
cdef class van_der_pol(sage.gsl.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
sage: T=ode_solver()
sage: T.algorithm="bsimp"
sage: vander=van_der_pol()
sage: T.function=vander
sage: vander=van_der_pol()
sage: T.ode_solve(y_0=[1,0],t_span=[0,2000],num_points=1000)
sage: T.plot_solution(i=0)
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to sage-devel@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at http://groups.google.com/group/sage-devel
URLs: http://sage.scipy.org/sage/ and http://modular.math.washington.edu/sage/
-~----------~----~----~----~------~----~------~--~---
