Why does scipy's odeint function give a non-monotonic solution for a problem whose solution should be monotone?

Question

The solution to the ode below looks like it is monotonically increasing:

However on closer inspection we see that it is not:

How can I ensure that the numerical solution is monotonically increasing?

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

def f1(theta):
    return -1*np.cos(theta)/np.tan(theta)
def f2(theta):
    return -1*np.cos(theta)

def f3(y, t) :
    return 1/np.sqrt(f1(y)**2+f2(y)**2)

s = np.linspace(0,30,100)
theta = odeint(f3, -np.pi/2, s)

plt.plot(s, theta)
plt.show()

Edit 1

As Geoff Oxberry mentions in his answer the RHS of the ODE is infinite at the initial condition. So I thought of finding another point to start integrating from which should be on the same trajectory as in the original case above. I now solve $\frac{ds}{d\theta} = f_3(\theta)$ instead of $\frac{d\theta}{ds} = \frac{1}{f_3(\theta)}$. The former does not have the singularity at the initial condition $s=0$ at $\theta = -\pi/2$. I then integrate to a arbitrary $\theta = -\pi/4$, and get the value of $s$ there. I then use this point as my initial condition in the original problem.

def f4(theta):
    return np.sqrt(f1(theta)**2+f2(theta)**2)

# Get another initial condition

theta0 = -np.pi/4
s0 =integrate.quad(f4, -np.pi/2, theta0) 

However I still get the problems of oscillations for this case, so I'm not sure what to make of this method.

Answer 1

If you want a numerical method to satisfy certain properties, you have to prove (or construct the method in such a way) that these properties hold.

Examples:

• conservative finite-volume formulations for conservation laws
• Runge-Kutta methods preserve linear invariants
• symplectic integrators preserve symplecticity, which is important for Hamiltonian systems
• symplectic Runge-Kutta methods preserve quadratic invariants
• total variation diminishing discretizations (TVD) are monotonicity-preserving (typically, these are spatial discretizations for hyperbolic problems)
• strong-stability preserving discretizations in time are those that preserve the TVD property, if a spatial discretization is TVD when integrated in time using forward Euler
• you can force certain integration schemes (example: BDF implementation in SUNDIALS) to respect bounds in properties of ODE solutions (example: you're integrating an ODE tracking the mass of an object; that object's mass must be nonnegative)

Tightening the tolerances, as DavidKetcheson correctly points out, will not necessarily make the numerical solution monotonic. There's no way to guarantee that. To summarize the following discussion, in a practical sense, if you are approaching a steady-state monotonically, tightening tolerances will limit spurious oscillations (specifically, it will limit their amplitude, but not necessarily their total variation), but it will not eliminate them for certain.

Your particular ODE is:

\begin{align} \dot{\theta} = \frac{1}{\sqrt{\frac{\cos^{4}{\theta}}{\sin^{2}{\theta}} + \cos^{2}}{\theta}} \end{align}

After some simplification and consulting trigonometric identities, it becomes

\begin{align} \dot{\theta} = \pm\tan{\theta}, \end{align}

and since you start at $\theta(0) = -\pi/2$, your right-hand side is singular anyway, and your "solution" is dubious. (The reason for the ambiguity in sign via $\pm$ involves a trigonometric identity for sine in terms of cotangent, and the choice of sign depends on which quadrant $\theta$ is in.) The only reason you see a solution at all is due to quirks in floating-point arithmetic; when I evaluate your right-hand side at zero, I get 16331239353195370.0, probably a system- and hardware-dependent result.

For $\theta \in (-\pi/2, 0]$, you should have $\dot{\theta} = -\tan{\theta}$, and since $-\tan{\theta}$ is nonnegative on $(-\pi/2, 0]$ and $\tan{0} = 0$, any initial condition in $(-\pi/2, 0]$ should yield a monotonically increasing solution. Again, this property isn't guaranteed by the numerical methods you are using. Furthermore, if $\theta \in (0, \pi/2)$, your equation becomes $\dot{\theta} = \tan{\theta}$ to preserve positivity of the right-hand side (the original form of the equation ensures positivity), and $\tan{\pi/2}$ is of course, again, singular.

The main function of your tolerance here is to ensure that you get the right steady-state, and because floating-point arithmetic is mysterious, you somehow get what looks like a "correct" answer (for $\theta(0) = -\pi/2 + \varepsilon$ for some small $\varepsilon > 0$) despite the initial condition being at a singularity, and you see those deviations from the true steady state at $\theta = 0$ is because tolerances correspond to "roughly zero" in terms of estimating error in the computed solution. Due to quirks in floating-point arithmetic, your right-hand side won't necessarily evaluate to exactly zero, but because it will be sufficiently small relative to your tolerances, you should see no oscillations of amplitude larger than your tolerances (approximately, because this reasoning is heuristic, and not rigorous, but it should hold).

Aside from reasons of stability or necessary precision, I am skeptical of extremely small tolerances unless they are required due to physical or numerical reasons. Undue demands on accuracy will unnecessarily restrict time step size and increase simulation execution times without yielding any appreciable improvement in results. As a counterexample, combustion problems are notorious for requiring small tolerances for certain species to get approximately correct numerical results; if certain species are never present, certain reactions cannot occur, even if those species concentrations are calculated in quantities unlikely to be measured, so in that case, small absolute and relative tolerances are warranted. For your particular problem, I suspect that with sufficiently large tolerances, you will see wild oscillations because you will overshoot appreciably your steady state and reach the singularity at $\theta = \pi/2$, so the default tolerance is probably fine, but I see no good reason to decrease the tolerances here.

Answer 2

I suggest reducing the tolerance of odeint. Set the "rtol" or "atol" parameters to something smaller than their default which is 1.49012e-8. See the scipy documentation for odeint here.