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?
