Numerical integration of ODEs: Why does higher accuracy and precision not lead to convergence?

Question

I'm trying to integrate numerically a non-linear ODE that, on paper, is simple (just the damped, forced pendulum), but that, as well known, its dynamics displays a chaotic behavior and is very sensitive to initial conditions. On longer time spans, it becomes difficult to obtain a convergent answer. This is expected, due to the nature of the dynamical system. But my problem is that increasing the numerical accuracy and precision, even by several orders of magnitude, leads only to a slightly better convergence. And I'm wondering why?

Some details follow.

This is the second order ODE of the damped and forced pendulum: \begin{equation} \ddot{\theta} + 2\beta \dot{\theta} + \omega_{0}^2 \sin \theta = \gamma \omega{_0}^2 \cos{\omega t}, \end{equation} with $\theta$ the angular excursion, the term on the rhs is the harmonic forcing, and where $\omega=2\pi$, $\omega_0=\frac{3}{2}\,\omega=3\pi$, $\beta=\frac{\omega_0}{4}=\frac{3\pi}{4}$, $\gamma=1.16$.

This can be reduced to a first-order ODE by the substitution ($y_{1}=\theta $): $$y'_{1}=y_2.$$

Then it becomes: \begin{equation} y'_{2} + 2\beta y_{2} + \omega{_0}^2 \sin y_{1} = \gamma \omega{_0}^2 \cos{\omega t}, \end{equation}

with initial conditions $y_{1}(t_0)=\theta(t_0)$ and $y_{2}(t_0)=y'_{1}(t_0)=\dot{\theta}(t_0)$.

I use a Bulirsch-Stoer integrator in Matlab which is supposedly the most adequate for high precision numerical integrations (except for functions with singularities, which isn't the case here.) Indeed, it performs better than a RK, but not that much as expected. The picture shows the numerical solutions for different tolerances (the accuracy of the solution at the grid-points.) Until about $27$ seconds (dynamical time of the pendulum), all solutions are the same, but then begin to diverge for the different tolerances. Of course, this is expected. For different accuracies, the solutions must turn out differently beyond a certain time interval. What was not expected that no matter what tolerance is applied, they all begin to diverge at about the same time. I expected that applying lower tolerances (e.g., here four tolerances, from $10^-20$ to $10^-23$) one could obtain accurate solutions going beyond the $27$ second limit. Not so. One factor is that double precision isn't sufficient. So I implemented a multiprecision tool that allows to calculate with arbitrary numbers of digits. I went so far in using 300 digits!! This allows to go further by about ten seconds, but then that's it. No matter how much accuracy and precision one applies, convergence doesn't work beyond a $40$ seconds barrier, either. I wondered whether the integrator works and does its job. Solving ODEs with known analytic solutions, I could check that it does indeed do the right thing. So, I think that there is something fundamentally flawed in my procedure (or my thinking.) My question is: Why does one not observe converge the solutions from a certain level of accuracy and precision onwards? What could be the underlying cause for such a behavior, or what is my mistake?

Answer 1

This is the nature of chaos. the number of digits of precision you need to obtain convergence to time t is often exponential in t.

That said, using better integrators, it is possible to do better.

using OrdinaryDiffEq, Plots
function pend(dy, y, p, t)
    (;gamma, w, w0, B) = p
    dy[1] = y[2]
    dy[2] = gamma*w0^2*cospi(w*t) - w0^2*sin(y[1]) - 2*B*y[2]
end
# guessing at initial conditions because you didn't give them in the post
prob = ODEProblem(pend, BigFloat[0, 0], (big"0.0",50), (gamma=1.16, B=3*big(pi)/4, w0=3*big(pi), w=2))
sol2 = solve(prob, Vern7(), abstol=1e-20, reltol=1e-20);
sol2 = solve(prob, Vern7(), abstol=1e-22, reltol=1e-22);
plot(sol1, idxs=2) # blue
plot!(sol2, idxs=2) # orange

Answer 2

As Oscar states, one of the definitions of chaos is that errors in a trajectory are magnified exponentially in time. However, note that although you lose accuracy on the trajectory, the strange attractor itself is stable so you will still simulate a trajectory on the attractor (or morely likely, some approximation of the attractor). This can be used to compute properties of the attractor such as time-averages, delay embeddings, Lyapunov exponents, etc. This is done even on mission critical simulations such as turbulent fluid simulations for aeronautics.

Answer 3

At some point, assuming your y scale is in radians, two differing solutions get arbitrarily close to an unstable inverted configuration. The solutions diverge noticeably when one solution tips back to rotate before going over center, while the other one makes it over. But even before that, the two solutions are slowly diverging in a less obvious way.

The time horizon for solution divergence in a chaotic system is on the order of

$\frac{1}{\lambda} ln\frac{a}{|\delta_0|}$ where

$a$ is a measure of error tolerance, $\lambda$ is a Lyapunov exponent that governs the growth of the error, and $\delta_0$ is a measure of the amplitude of the intial error.

So, as others have noted, it's an unfortunate mathematical fact that vastly improving integration accuracy does not result in a very large gain in "accurate" solution time. Because its all under the natural log.

Source for horizon equation and further discussion: "Nonlinear Dynamics and Chaos" by Steve Strogatz.