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 |
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 |
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 |
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 > 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 |
On Tue, Oct 26, 2010 at 12:38 PM, Sebastian Castillo <[hidden email]> wrote:
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 |
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 |
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 |
Free forum by Nabble | Edit this page |