Why does LSODA fail to integrate the logistic function?

Question

I'm comparing some of the different ODE integrators in scipy.integrate.ode on solving the logistic function:

$$x(t) = \frac{1}{1+e^{-rt}}$$ $$\dot{x} = rx(1-x)$$

I've heard that LSODA should be very good, so I was a bit surprised to find that it fails completely for $r>4$. The dopri5 integrator, on the other hand, seems to have no problem for very large values of $r$ (see plots below).

Why does LSODA fail here? Do I need to tune some parameters, or is it simply a poor choice for this particular problem?

Here's my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

# Logistic function
r = 5
def func(t):
    return 1. / (1 + np.exp(-r * t))

def dfunc(t, X):
    x, = X
    return [r * x * (1 - x)]

t0 = -10
t1 = 10
dt = 0.01
X0 = func(t0)

integrator = 'lsoda'

t = [t0]
Xi = [func(t0)]
ode = integrate.ode(dfunc)
ode.set_integrator(integrator)
ode.set_initial_value(X0, t0)

while ode.successful() and ode.t < t1:
    t.append(ode.t + dt)
    Xi.append(ode.integrate(ode.t + dt))

t = np.array(t)     # Time
X = func(t)         # Solution
Xi = np.array(Xi)   # Numerical

# Plot analytical and numerical solutions
X = func(t)

plt.subplot(211)
plt.title("Solution")
plt.plot(t, X, label="analytical", color='g')
plt.plot(t, Xi, label=integrator)
plt.xlabel("t")
plt.ylabel("x")
plt.legend(loc='upper left', bbox_to_anchor=(1.0, 1.05))

# Plot numerical integration errors
err = np.abs(X - Xi.flatten())
print("{}: mean error: {}".format(integrator, np.mean(err)))

plt.subplot(212)
plt.title("Error")
plt.plot(t, err, label=integrator)
plt.xlabel("t")
plt.ylabel("error")
plt.legend(loc='upper left', bbox_to_anchor=(1.0, 1.05))

plt.tight_layout()
plt.show()

Result for integrator = 'lsoda':

lsoda: mean error: 0.500249750249742

Result for integrator = 'dopri5':

dopri5: mean error: 3.7564128626655536e-11

Answer 1

When you use $r=5$, the initial condition is $$ x(-10) \approx e^{-50} \approx 1.9\times 10^{-22}. $$ This is much smaller than the machine epsilon, $2\times 10^{-16}$, and it is very likely that LSODA just concludes that the solution is identically zero.

To check this idea, I replaced $t_0$ with $-5$ instead of $-10$, so that $x(-5)\approx 1.4\times 10^{-11}$, and then it seems to work. The solution isn't very accurate, the error rises to $0.01$, but this is likely because the numerical error places the location of the bump in slightly the wrong place, while the shape of the solution is correct.

The way LSODA seems to decide that the solution is identically zero is that it uses its absolute tolerance, and presumably $x(-10)$ is smaller than that, causing the problem. So to fix that, instead of changing $t_0$, one can also replace the line

ode.set_integrator(integrator)

with

ode.set_integrator(integrator, atol=0.0)

Then it also seems to work, but the error is $10^{-5}$. To get the error down to the same level as DOPRI5, I can set rtol=1e-10, but I'm not quite sure why DOPRI5 manages with its default parameters, or whether those defaults are different. In any case, I think it's the default tolerances at work here.