Very nice! Things seem fine so far after a few crude tests (on a Mac
Pro, sage 1.6, no packages installed). I will try to push it harder
soon. I think there is one small typo in the documentation: for the
initial condition example shouldn't it be
y_0=1,y_1=1, would do
sage: T.y_0=[1,1]
instead of
y_0=1,y_1=1, would do
sage: T.y_0=[1,0]
Cheers,
Marshall Hampton
Joshua Kantor wrote:
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
