SciPy ODE integrator

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

SciPy ODE integrator

David Pine
Is there a SciPy ODE integrator that does adaptive stepsize integration AND produces output with the adaptive time steps intact?

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
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy ODE integrator

Anne Archibald-2
On 26 July 2010 12:46, David Pine <[hidden email]> 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
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy ODE integrator

David Pine
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 <[hidden email]> 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
>> [hidden email]
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy ODE integrator

Sebastian Castillo-3
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!


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy ODE integrator

Warren Weckesser-3


On Tue, Oct 26, 2010 at 12:38 PM, Sebastian Castillo <[hidden email]> 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
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
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy ODE integrator

Sebastian Castillo-3
Warren Weckesser <warren.weckesser <at> enthought.com> writes:

> Sebastian,Here's an example that uses the 'step=True' option of the
integrate() method:-----from numpy import array
>
> from scipy.integrate import odefrom pylab import figure, show, xlabel,
ylabelfrom mpl_toolkits.mplot3d import Axes3Ddef 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 fic = [1.0, 2.0, 1.0]t0 = 0.0t1 = 100.0#dt = 0.1sigma = 10.0rho =
28.0beta = 10.0/3solver = 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 listSciPy-User <at>
scipy.orghttp://mail.scipy.org/mailman/listinfo/scipy-user

Warren:

Thank you for your quick answer! Now it works with regular nonstiff-problems,
but I'm still having problems with the stiff ones. I am trying to solve a stiff
van-der-pol oscillator as described in matlab's ode15s reference (finish time
tf=3000). The problem is that it takes too long, so long that I have to kill
the process (in fact, if I change the finish time to tf=10, it takes
approximately 1 minute to solve it whereas matlab takes 0.2s to solve it for
tf=3000, and odeint took even less but with specified time points, which I don't
want). Here is my code:

def f_test(t,y):
    dy=numpy.zeros(2)
    dy[0]=y[1]
    dy[1]=1000*(1 - y[0]**2)*y[1] - y[0];
    return dy

import numpy
import scipy
import scipy.integrate
import matplotlib.pyplot as plt

# Define initial conditions
t0=numpy.array(0.0)
tf=numpy.array(3000.0)
y0=numpy.array([2,0])
y=y0
t=t0

# Solve
r=scipy.integrate.ode(f_test).set_integrator('vode', method='bdf', order=15)
r.set_initial_value(y0,t0)
while r.successful() and r.t<tf:
    r.integrate(tf,step=True)
    t=numpy.hstack((t,r.t))
    y=numpy.column_stack((y,r.y))

# Plot
plt.plot(t,y[0,:],'-o')

Hope you or anyone can help. Thank you!

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: SciPy ODE integrator

Rob Clewley
Sebastian,

> Thank you for your quick answer! Now it works with regular nonstiff-problems,
> but I'm still having problems with the stiff ones. I am trying to solve a stiff
> van-der-pol oscillator as described in matlab's ode15s reference (finish time
> tf=3000). The problem is that it takes too long, so long that I have to kill
> the process (in fact, if I change the finish time to tf=10, it takes
> approximately 1 minute to solve it whereas matlab takes 0.2s to solve it for
> tf=3000, and odeint took even less but with specified time points, which I don't
> want). Here is my code:


In case it helps, PyDSTool's Radau stiff integrator will solve this
system much faster than Matlab (with adaptive time step, of course),
and there is a van der pol example (actually involving a bifurcation
analysis) provided in the /tests directory, as PyCont_vanDerPol.py.

Let me know if you need help setting up. The full wiki documentation
pages are offline right now (the old server crashed) but some of the
core pages are being mirrored at my web site, linked from
pydstool.sourceforge.net. These do include installation information
for the Dopri and Radau integrators.

-Rob
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user