7
$\begingroup$

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 with lsoda

Result for integrator = 'dopri5':

dopri5: mean error: 3.7564128626655536e-11

result with dopri5

$\endgroup$
8
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you, that makes sense. It would be interesting to understand why DOPRI5 manages this without problem though! $\endgroup$ – joh Oct 12 '16 at 14:38
  • 1
    $\begingroup$ @joh My guess is they use different error control schemes. $\endgroup$ – Kirill Oct 12 '16 at 23:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.