I have a problem in physics formulated via an ODE. Now I like to solve it numerically using Pythons scipy.integrate and the therein complex_ode. I figured out how and it works but now I like to optimize my code and I ran into a question:

I need a massive amount of grid points so that the solution converges. So I tried setting another integrator method with stepsize control (dop853). I am not really familiar with that feature but it works as well and is much faster than my approach. Here's a simple MWE which does not reflect my complicated ODE but the structure of the setup (complex valued etc.) is the same:

from scipy import *
from scipy.integrate import *
from pylab import *

grid = linspace(0, 8, 10e3)
dgrid = (grid[1]-grid[0])*ones(100e3)

def RHS(t, x):
    return -x

y = zeros(len(grid), dtype=complex)
y[0] = 1.0

solver = complex_ode(RHS)
solver.set_initial_value(y[0], grid[0]).set_integrator('dop853')

for t in range(len(grid)):
    y[t] = solver.y[0]

plot(grid, y.real)

But now the question: I now of course want to relate the time grid to the solution. What grid is now connected to the solution? Is it still the grid I was introducing (grid) or does the stepsize control calculate at different time steps and I would need another time array?


In this situation, you need an integrator with dense output, which means it has a capability to calculate the solution to the ODE at any point in time between the time steps taken by the integrator. dop853 (from Hairer and Wanner's eighth-order Dormand-Prince integrator) has this capability, in addition to adaptive time-stepping. Likely what happens is that the integrator keeps track of enough time steps internally to be able to interpolate solutions at intermediate times -- that is, the integrator has an "internal discretization". When you call the integrator using solver.integrate(solver.t + dt), it adjusts the internal discretization (if necessary) so that the solution at solver.t + dt can be interpolated from solution values calculated in the internal discretization. If you can read the original Fortran source, the part around "SAVE THE FIRST FUNCTION EVALUATION" is probably the relevant code illustrating this strategy.

So, to answer your question, the stepsize control generally does calculate at different time steps, but you don't need another time array for that. Just call the SciPy interface to dop853 directly to calculate the solution at the time points you want, and it will take care of everything.

  • $\begingroup$ That's how I interpreted the documentation as well; it's a bit opaque. I think if one wanted to keep a record of every internal timestep they'd need to use the callback facility in the integrator. $\endgroup$ – Aurelius Jan 14 '14 at 19:28
  • $\begingroup$ I think I am not free in choosing the time array I want. If you try solver.integrate(grid[t]) and the default integrator vode then you will not a solution. Instead, you will see: DVODE-- ISTATE (=I1) .gt. 1 but DVODE not initialized In above message, I1 = 2 $\endgroup$ – DaPhil Jan 15 '14 at 8:55
  • $\begingroup$ @DaPhil Could you clarify what you modified in the code sample above? Did you just substitute vode for dop853? $\endgroup$ – Geoff Oxberry Jan 15 '14 at 19:38

This site is temporarily in read only mode and not accepting new answers.

Not the answer you're looking for? Browse other questions tagged .