[SciPy-User] SciPy ODE integrator

Sebastian Castillo castillohair@gmail....
Tue Oct 26 17:29:59 CDT 2010


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!



More information about the SciPy-User mailing list