Python SciPy ODE solver not converging

Question

I'm trying to use scipy's ode solver to plot the interaction between a 2D system of equations. I'm attempting to alter the parameters passed to the solver by the following block of code:

# define maximum number of iteration steps for ode solver iteration
m = 1 #power of iteration
N = 2**m #number of steps

# setup a try-catch formulation to increase the number of steps as needed for solution to converge
while True:
    try:
        z = ode(stateEq).set_integrator("vode",nsteps=N,method='bdf',max_step=5e5)
        z.set_initial_value(x0, t0)
        for i in range(1, t.size):
            if i%1e3 == 0:
                print 'still integrating...'
            x[i, :] = z.integrate(t[i]) # get one more value, add it to the array
            if not z.successful():
                raise RuntimeError("Could not integrate")
        break
    except:
        m += 1
        N = 2**m
        if m%2 == 0:
            print 'increasing nsteps...'
            print 'nsteps = ', N

Running this never breaks the while loop. It keeps increasing the nsteps forever and the system never gets solved. If I don't put it in the while loop, the system gets solved, I think, because the solution gets plotted. Is the while loop necessary? Am I formulating the solver incorrectly?

*If I don't put it in the while loop, the system gets solved, I think, because the solution gets plotted.* – I fail to make sense of this sentence. — Wrzlprmft, Aug 03 '18 at 08:34

score 2 · Accepted Answer · answered Aug 03 '18 at 08:34

The parameter nsteps regulates how many integration steps can be maximally performed during one sampling step (i.e., a call of z.integrate). Its default value is okay if your sampling step is sufficiently small to capture the dynamics. If you want to integrate over a huge time span in one large sampling step (e.g., to get rid of transient dynamics), the value can easily be too small.

The point of this parameter is to avoid problems arising from unexpectedly very long integrations. For example, if you want to perform a given integration for 100 values of a control parameter in a loop over night, you do not want to see on the next morning that the No. 14 was pathological and is still running.

If this is not relevant to you, just set nsteps to a very high value and stop worrying about it. There is certainly no point to successively increase nsteps, you are just performing the same calculations all over again.

Running this never breaks the while loop. It keeps increasing the nsteps forever and the system never gets solved.

This suggests that you have a different problem than nsteps being exceeded, most likely that the problem is not well posed. Carefully read the error message produced by the integrator. I also recommend that you check your differential equations. It may help to look at the solutions until the integration fails to see what is going wrong, i.e., plot x after running this:

z = ode(stateEq)
z.set_integrator("vode",nsteps=1e10,method='bdf',max_step=5e5)
z.set_initial_value(x0, t0)

for i,time in enumerate(t):
    x[i,:] = z.integrate(time)
    if not z.successful():
        break

Your value for max_step is very high (this should not be higher than the time scale of your dynamics). Depending on your application, this may very well be reasonable, but then it suggests that you are working with large numbers. This in turn may mean that the default values of the parameters atol and first_step are not suited for your situation and you need to adjust them.

Python SciPy ODE solver not converging

1 Answers1