Python: Grid with step control ODE solver

Question

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)):
    solver.integrate(solver.t+dgrid[t])
    y[t] = solver.y[0]

plot(grid, y.real)
show()

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?

Answer 1

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.