I am dealing with a Double Pendulum with two masses of size $1$ linked with inelastic, weightless rods. There is assumed to be no air resistance too.
The Hamiltonian $H$ satisfies the following system:
$$ \dot{\theta}_1 = \frac{\partial H}{\partial p_1} \ \ \ \ \ \ \dot{\theta}_2 = \frac{\partial H}{\partial p_2} \ \ \ \ \ \ \dot{p}_1 = -\frac{\partial H}{\partial \theta_1} \ \ \ \ \ \ \dot{p}_2 = -\frac{\partial H}{\partial \theta_2} $$
where $p_1$ and $p_2$ are the momenta of the two masses. All of the above are functions that take in inputs $\theta_1, \theta_2, p_1$ and $p_2$. These functions can be written as one vector function i.e. $ \ \ \ \textbf{F}(t,\textbf{y}) = \left[\dot{\theta}_1 \ \ \dot{\theta}_2 \ \ \dot{p}_1 \ \ \dot{p}_2\right]^\top$ that takes inputs time $t$ and the parameters $\theta_1, \theta_2, p_1$ and $p_2$ in a form of a vector, $\textbf{y}$. Take
$$ \textbf{y}_n = \begin{bmatrix} \theta^n_1 \\ \theta^n_2 \\ p^n_1 \\ p^n_2 \end{bmatrix} $$
where $\theta^n_1$ would be the approximation of $\theta_1$ at time $t_n$ and the same principle applies to the other variables. Now the Implicit Midpoint Rule says:
$$ \textbf{y}_{n+1} = \textbf{y}_{n} + \Delta t \textbf{F}\left(t_n + \frac{\Delta t}{2}, \frac{\textbf{y}_{n+1}+\textbf{y}_{n}}{2} \right) $$
As this is an Implicit Method, I have tried using Fixed Point Iteration to eventually get values $\theta^{n+1}_1, \theta^{n+1}_2, p^{n+1}_1$ and $p^{n+1}_2$ at the same time, that is
$$ \begin{bmatrix} \theta^{n+1}_1 \\ \theta^{n+1}_2 \\ p^{n+1}_1 \\ p^{n+1}_2 \\ \end{bmatrix} = \begin{bmatrix} \theta^{n}_1 \\ \theta^{n}_2 \\ p^{n}_1 \\ p^{n}_2 \\ \end{bmatrix} + \textbf{F}\left( t_n + \frac{\Delta t}{2}, \ \ \frac{1}{2} \begin{bmatrix} \theta^{n+1}_1 + \theta^{n}_1 \\ \theta^{n+1}_2 + \theta^{n}_2\\ p^{n+1}_1 + p^{n}_1 \\ p^{n+1}_2 + p^{n}_2 \\ \end{bmatrix} \right) $$
i.e. looping through 4 systems and inputting the newly acquired approximations (to each system) each time until a tolerance is reached.
But this has proven unsuccessful as some parameters shoot up to infinity.
Ultimately, I am investigating the effect of the initial value of $\theta_2$ to the time it takes for $\theta_2\geq \pi$ or $\theta_2\leq -\pi$ (If anyone can provide me with what I should be expecting, that would be great!), so I tried applying a Fixed Point Iteration for only $\theta_2$ when applying the Implicit Midpoint Rule, but this too has proven unsuccessful and also seems like an erroneous approach.
What would be the correct approach to this method in order to produce correct results?
NOTE: I have purposefully omitted the actual functions of the partial derivatives and $H$ itself as I think they are irrelevant to my question and might make this too long. I have tested the functions with my (C++) code and they work fine. I will not include my code, however, will be happy to provide it (I have made it readable).
