Let's start out using notation similar to the example you linked to:
$$
\ddot{y}+b\dot{y}+\sin(y)=a\cos(ct)+d
$$
As in the example, we'll write this in autonomous form:
$$
\mathbf{x}=\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}=\begin{pmatrix}y\\\dot{y}\\ct\end{pmatrix}
$$
the differential equation can now be written:
$$
\mathbf{\dot{x}}=\begin{pmatrix}\dot{x_1}\\\dot{x_2}\\\dot{x_3}\end{pmatrix}=\mathbf{f(x)}=\begin{pmatrix}x_2\\-bx_2-\sin{x_1}+a\cos{x_3}+d\\c\end{pmatrix}
$$
Now, suppose we study the evolution of two different starting points in this system that are arbitraily close. Let's denote the trajectory of one of them with $\mathbf{x}(t)$ and the other one with $\mathbf{x}(t)+\mathbf{z}(t)$, where $\mathbf{z}(t)$ is the vector difference between the two points. The time derivative of $\mathbf{z}(t)$ can be written:
$$
\mathbf{\dot{z}}=\mathbf{f(x+z)}-\mathbf{f(x)}
$$
If $\mathbf{z}$ is sufficiently small this can be written:
$$
\mathbf{\dot{z}}=\mathbf{Jz}
$$
where $\mathbf{J(x)}$ is the Jacobian of $\mathbf{f(x)}$:
$$
\mathbf{J}=\begin{pmatrix} 0 & 1 & 0\\ -\cos{x_1} & -b & -a\sin{x_3}\\ 0 & 0 & 0 \end{pmatrix}
$$
Over a short period $\delta t$, the change in $\mathbf{z}$ is given by:
$$
\mathbf{z}(t+\delta t)=\mathbf{z}(t)+\delta t\mathbf{Jz}(t)=(\mathbf{I} +\delta t\mathbf{J})\mathbf{z}(t)
$$
Assuming $\mathbf{z}$ and $\delta t$ are small enough for the linearizations to hold we can decribe the change in $\mathbf{z}$ from time $t_1$ to $t_2$ with:
$$
\mathbf{z}(t_2)=\left[\prod_i(\mathbf{I} +\delta t_i\mathbf{J}(\mathbf{x}_i))\right]\mathbf{z}(t_1)
$$
where $\delta t_i$ are sufficiently small time intervals from $t_1$ to $t_2$ and $\mathbf{x}_i$ is $\mathbf{x}$ at those times.
What you need to do to calculate the Lyaponov exponents of the this system is to numerically solve the trajectory $\mathbf{x}(t)$ for an arbitrary starting point $\mathbf{x}_0$ over a long time interval $T$. Then you need to calculate the matrix $\mathbf{A}$:
$$
\mathbf{A}=\prod_i(\mathbf{I} +\delta t_i\mathbf{J}(\mathbf{x}_i))
$$
using the values for $\delta t_i$ and $\mathbf{x}_i$ from the solved trajectory. You then need to calculate the eigenvalues of $\mathbf{A}$. Let's call them $\alpha_1$, $\alpha_2$ and $\alpha_3$. The lyaponov exponents are given by:
$$
\lambda_i=\frac{1}{T}\log{\alpha_i}
$$
assuming $T$ is long enough. You need to experiment with different $\mathbf{x}_0$ and $T$. Hopefully, the variations in the Lyaponov exponents you calculate for different $\mathbf{x}_0$ will decrease with larger $T$ and converge for sufficiently large $T$.
One problem you will almost certainly run into is that the values of $\mathbf{A}$ will become too big to handle numerically. The way to get around this is to keep track of the magnitude of $\mathbf{A}$ as it is iteratively calculated and normalize it whenever its highest value reaches some threshold. The normalizing factors needs to be saved each time this is done. Assuming this has been done you will in the end have a matrix $\mathbf{B}$ and a series of normalizing factors $N_j$ (in my notation $N_j$ are the values you have divided by when normalizing) relating to $\mathbf{A}$ as:
$$
\mathbf{A}=\left(\prod_j N_j\right)\mathbf{B}
$$
If $\beta_i$ are the eigenvalues of $\mathbf{B}$, the Lyaponov exponents can be calculated as:
$$
\lambda_i=\frac{1}{T}\left(\log{\beta_i}+\sum_j\log{N_j}\right)
$$
It is usually the largest of the Lyaponov exponents that is of interest as it will dominate the exponential growth in $\mathbf{z}$. At least one of the Lyaponov exponents must be positive for the system to be chaotic.
It is important to understand that $\mathbf{z}$ should be thought of as a vector that remains very small throughout the whole period $T$. This might be a bit counter-intuitive since it grows exponentially. However, for any finite $T$, the starting value of $\mathbf{z}$ can always be chosen small enough to make sure it reamins sufficiently small thoughout the time period.
First. Yes, it is possible for the curves in Poincare section to cross. (I am assuming you mean generally). Remember, that Poincare section is a 2D projection of a 3D section of a 4D phase space. Regular, non-chaotic dynamics correspond to the winding of a $T_2$ torus embedded in this 4D space. Sections of such torus could be the couple of closed curves in the 3D hyperplane. Projecting this hyperplane on the 2D space could produce intersections. Additionally, Poincare section may include only points with a given directions of a hypersurface crossing, then the number of resulting curves would be reduced. Probably even then there still is a possibility for the intersecting curves in the Poincare section, although those would require highly nonlinear regimes (possibly after at least one bifurcation).
Second. Why would you expect the symmetry with respect to sign of the momentum? The Hamiltonian for the double pendulum contains the term proportional to the product of both momenta $l_\alpha l_\beta$. Therefore Hamiltonian and the torus in the 4D phase space are invariant only with respect to change of signs of both momenta,
but not with respect to a change of sign of a single momentum. After making cross-section and projecting it onto 2D plane, no mirror symmetry w.r.t. momentum sign is to be expected in the resulting Poincare section.
Addition Upon discussion in the comments, I understand that the one unresolved issue is the set of initial data which generates the figure from the book.
In the book the figure is accompanied by the comment:
To produce Figure 11.4 and all subsequent Poincare sections,
the initial conditions are varied. The total energy has the value $E = H(\alpha(0), 0, -\alpha(0), 0)$.
So, all those curves correspond to different sections arising from varied ICs with the fixed energy $E_0$. Here is my attempt at reproducing the figure:
Various colors correspond to an individual Poincare section defined by its own initial point. The point are chosen so that:
$$
\alpha=0, \quad \beta=0, \quad \dot{\alpha}>0, \quad H(0,l_\alpha,0,l_\beta)=E_0.
$$
The value of $l_\beta$ is chosen so that curves are evenly spaced between two limiting points corresponding to normal modes of oscillation. The last two conditions and known $l_\beta$ allow us to unambiguously calculate $l_\alpha$.
Finally, the limiting curve is defined by the equations:
$$
\alpha=0, \quad \dot{\alpha}=0, \quad H(0,l_\alpha,0,l_\beta)=E_0.
$$
The code for the picture is based on the linked question from Mathematica SE and can be seen here.
Best Answer
You are making to the following mistakes:
Your initial displacement $d_0$ is 5 % of the variation of the respective variable, while it should be orders of magnitude lower. The Lyapunov exponent is defined for the limit of infinitesimal displacements, i.e., $d_0→0$.
Your observation time is far too short. You are looking at four oscillations of your dynamics, when you could have thousands. The Lyapunov exponent is defined as the average over the whole attractor/trajectory.
You do not consider multiple initial displacements or rescale the displacement vector. If you do not do this, the displacement will eventually grow to the size of the attractor and become meaningless.
In your specific example, the following happens:
You seem to be using a pendulum without small-angle approximation and thus the frequency of your pendulum depends on the intial condition.
Thus, your displaced system (red) has a slightly lower frequency than the original one and thus both systems dephase over time (at least that’s all they do in your short observation time). You see this as an increase of the distance which leads to your positive Lyapunov exponent.
If you had used an infinitesimal excitation, the different frequencies would had little impact and you would have to be very precise to see them. For methods calculating the Lyapunov exponent in tangent space (e.g., Benettin et al.’s), this would have been no problem at all.
If you had averaged over multiple initial displacements and, you would also have had displacements in which the phases would first align rather than diverge and I would intuitively say that on average you would have obtained a zero Lyapunov exponent (I cannot prove this though and it’s a rather meaningless excercise given the other mistakes).
The single pendulum should have a zero (largest) Lyapunov exponent like all limit cycle dynamics. You can understand this as follows: A displacement along the trajectory will neither shrink nor grow in time (on average). The two time series will be exactly the same except for an offset corresponding to the displacement.
The absolute value of the Lyapunov exponent is rather meaningless, as it depends on how you scale time. Without knowing the detailed parameters of your double pendulum, it’s impossible to say what you should expect. The only thing that is meaningful is the sign of the largest Lyapunov exponent.