scipy.integrate.ode gives up on integration

Question

I am encountering a strange problem with scipy.integrate.ode. Here is a minimal working example:

import sys
import time
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from scipy.integrate import ode, complex_ode

def fun(t):
    return np.exp( -t**2 / 2. )

def ode_fun(t, y):
    a, b = y
    f = fun(t)
    c = np.conjugate(b)
    dt_a = -2j*f*c + 2j*f*b
    dt_b = 1j*f*a
    return [dt_a, dt_b]

t_range = np.linspace(-10., 10., 10000)

init_cond = [-1, 0]
trajectory = np.empty((len(t_range), len(init_cond)), dtype=np.complex128)

### setup ###
r = ode(ode_fun).set_integrator('zvode', method='adams', with_jacobian=False)
r.set_initial_value(init_cond, t_range[0])
dt = t_range[1] - t_range[0]

### integration ###
for i, t_i in enumerate(t_range):
    trajectory[i,:] = r.integrate(r.t+dt)

a_traj    = trajectory[:,0]
b_traj    = trajectory[:,1]
fun_traj  = fun(t_range)

### plot ###
plt.figure(figsize=(10,5))
plt.subplot(121, title='ODE solution')
plt.plot(t_range, np.real(a_traj))
plt.subplot(122, title='Input')
plt.plot(t_range, fun_traj)
plt.show()

This code works correctly and the output figure is (the ODE is explicitly dependent on the input variable, right panel shows the input, left panel the ode solution for the first variable).

So in principle my code is working. What is strange is that if I simply replace the integration range

t_range = np.linspace(-10., 10., 10000)

by

t_range = np.linspace(-20., 20., 10000)

I get the output

So somehow the integrator just gave up on the integration and left my solution as a constant. Why does this happen? How can I fix it?

Some things I've tested: It is clearly not a resolution problem, the integration steps are really small already. Instead, it seems that the integrator does not even bother calling the ode function anymore after a few steps. I've tested that by including a print statement in the ode_fun().

My current suspicion is that the integrator decided that my function is constant after it did not change significantly during the first few integration steps. Do I maybe have to set some tolerance levels somewhere?

Any help appreciated!

Answer 1

"My current suspicion is that the integrator decided that my function is constant after it did not change significantly during the first few integration steps." Your suspicion is correct. ODE solvers typically have an internal step size that is adaptively adjusted based on error estimates computed by the solver. These step sizes are independent of the times for which the output is requested; the output at the requested times is computed using interpolation of the solution at the points computed at the internal steps.

When you start your solver at t = -20, apparently the input changes so slowly that the solver's internal step size becomes large enough that by the time the solver gets near t = 0, the solver skips right over the input pulse.

You can limit the internal step size with the option max_step of the set_integrator method. If I set max_step to 2.0 (for example),

r = ode(ode_fun).set_integrator('zvode', method='adams', with_jacobian=False,
                                max_step=2.0)

I get the output that you expect.