[SciPy-User] SciPy ODE integrator

Warren Weckesser warren.weckesser@enthought....
Tue Oct 26 13:35:56 CDT 2010


On Tue, Oct 26, 2010 at 12:38 PM, Sebastian Castillo <castillohair@gmail.com
> wrote:

> David Pine <djpine <at> gmail.com> writes:
>
> >
> > Anne,
> >
> > Thanks.  Actually I finally figured this (the VODE option) out but I
> agree
> that scipy's ODE solvers need a
> > makeover.  The routines under the hood seem to be quite nice but the
> interface
> to Python is clumsy at best and
> > the documentation on how to use it is pretty awful.  I'll take a look at
> pydstool.  Thanks.
> >
> > David
> >
> > On Jul 28, 2010, at 10:45 AM, Anne Archibald wrote:
> >
> > > On 26 July 2010 12:46, David Pine <djpine <at> gmail.com> wrote:
> > >> Is there a SciPy ODE integrator that does adaptive stepsize
> integration AND
> produces output with the
> > adaptive time steps intact?
> > >
> > > It is not obvious, but the object-oriented integrator, based on VODE,
> > > can be run in this mode. You normally tell it how much to advance on
> > > each call and it does as many adaptive steps as it takes to get there,
> > > but there is an optional argument you can pass it that will make it
> > > take just one step of the underlying integrator. You can then write a
> > > python loop to produce the solution you want.
> > >
> > > If this seems messy, I have to agree. scipy's ODE integrators are in
> > > desperate need of an API redesign (they've had one already, which is
> > > why there are two completely different interfaces, but they need
> > > another). You could try pydstool, which is designed for the study of
> > > dynamical systems and has many more tools for working with ODEs and
> > > their solutions.
> > >
> > > Anne
> > >
> > >> The standard SciPy ODE integrator seems to be scipy.integrate.odeint
> and
> its simpler cousin
> > scipy.integrate.ode.  These work just fine but  both take a
> user-specified
> time series and returns the
> > solution at those points only.  Often, I prefer to have a more classic
> adaptive stepsize integrator that
> > returns the solution at time steps determined by the integrator (and the
> degree of desired precision
> > input by the user).  This is often the most useful kind of solution
> because it
> tends to produce more points
> > where the solution is varying rapidly and fewer where it is not varying
> much.
>  A classic Runge-Kugga
> > adaptive stepsize ODE solver does this as to many others, but I can't
> find a
> nice implementation in SciPy or
> > NumPy.  Please advise.  Thanks.
> > >>
> > >> David
> > >> _______________________________________________
> > >> SciPy-User mailing list
> > >> SciPy-User <at> scipy.org
> > >> http://mail.scipy.org/mailman/listinfo/scipy-user
> > >>
> > > _______________________________________________
> > > SciPy-User mailing list
> > > SciPy-User <at> scipy.org
> > > http://mail.scipy.org/mailman/listinfo/scipy-user
> >
>
> Hello. I am trying to use scipy.integrate.ode with time steps determined by
> the
> integrator. I can see that you have acomplished this already. I understand
> that
> there is an option in the VODE object, and I have found a "step" method in
> it,
> but I still don't get how to use it. Can you please post some example code
> of
> this? Thank you very much!
>
>
Sebastian,

Here's an example that uses the 'step=True' option of the integrate()
method:

-----
from numpy import array
from scipy.integrate import ode
from pylab import figure, show, xlabel, ylabel
from mpl_toolkits.mplot3d import Axes3D

def lorenz_sys(t, q, sigma, rho, beta):
    x = q[0]
    y = q[1]
    z = q[2]
    f = [sigma * (y - x),
         rho*x - y - x*z,
         x*y - beta*z]
    return f


ic = [1.0, 2.0, 1.0]
t0 = 0.0
t1 = 100.0
#dt = 0.1

sigma = 10.0
rho = 28.0
beta = 10.0/3

solver = ode(lorenz_sys)

t = []
sol = []
solver.set_initial_value(ic, t0)
solver.set_integrator('vode')
solver.set_f_params(sigma, rho, beta)
while solver.successful() and solver.t < t1:
    solver.integrate(t1, step=True)
    t.append(solver.t)
    sol.append(solver.y)

t = array(t)
sol = array(sol)

fig = figure()
ax = Axes3D(fig)
ax.plot(sol[:,0], sol[:,1], sol[:,2])
xlabel('x')
ylabel('y')
show()
-----

Warren



>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20101026/cbadf5db/attachment.html 


More information about the SciPy-User mailing list