On 4/1/07, Joshua Kantor <[EMAIL PROTECTED]> wrote:
> I just installed a new version of sage and
> ode_solver? fails with the same errors as in his message. Did something
> change
> which would obviously cause this.
This is probably bug in the misc/sageinspect.py, which Nick Alexander
recently rewrote
(so the bug is either mine or his). In any case, a temporary work
around is to type
ode_solver?? which will show the documentation (and source code). I've reported
this as trac #340:
http://www.sagemath.org:9002/sage_trac/ticket/340
>
>
> Josh
>
>
> ---------- Forwarded message ----------
> From: Hamptonio <[EMAIL PROTECTED]>
> Date: Mar 31, 2007 6:37 AM
> Subject: Re: ode_solver
> To: Joshua Kantor <[EMAIL PROTECTED]>
>
> Since I hope to give a demo on sage to other faculty in my department,
> I've been taking a look at the ode_solver again. (If anyone reading
> this has suggestions on impressive sage demos let me know!).
> Anyway, I hit a funny bug: if I try to do: ode_solver? in the
> notebook, I get the following:
>
> Traceback (most recent call last):
> File "<stdin>", line 1, in <module>
> File
> "/Users/mh/sage-2.0/ODE1/worksheets/ode1/code/10.py", line
> 4,
> in <module>
> print _support_.docstring("ode_solver", globals())
> File
> "/Users/mh/sage-2.0/local/lib/python2.5/site-packages/sage/
> server/support.py", line 131, in docstring
> s += 'Definition: %s\n'%sageinspect.sage_getdef(obj,
> obj_name)
> File
> "/Users/mh/sage-2.0/local/lib/python2.5/site-packages/sage/misc/
> sageinspect.py ", line 264, in sage_getdef
> spec = sage_getargspec(obj)
> File
> "/Users/mh/sage-2.0/local/lib/python2.5/site-packages/sage/misc/
> sageinspect.py", line 246, in sage_getargspec
> return _sage_getargspec_sagex(source)
> File
> "/Users/mh/sage-2.0/local/lib/python2.5/site-packages/sage/misc/
> sageinspect.py", line 200, in _sage_getargspec_sagex
> raise ValueError, "Could not parse sagex argspec"
> ValueError: Could not parse sagex argspec
>
> This is with sage 2.4.1.2, upgraded from 2.0. The ode_solver? command
> works fine from the command line interpreter, and the actual functions
> are fine - i.e. I can use ode_solver in the notebook and it works
> fine, its just the documentation that is messed up.
>
> Is this a more general notebook problem, or is it specific to
> ode_solver?
>
> Cheers,
> Marshall Hampton
>
>
> On Jan 4, 1:32 am, "Joshua Kantor" < [EMAIL PROTECTED]> 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
> > 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)
>
>
--
William Stein
Associate Professor of Mathematics
University of Washington
http://www.williamstein.org
--~--~---------~--~----~------------~-------~--~----~
To post to this group, send email to [email protected]
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/
-~----------~----~----~----~------~----~------~--~---
