An implementation is presented and the advantages and disadvantages of this approach for dynamical systems is briefly discussed.

Available with Zymplectic v.0.4.1

The Hamiltonian of the spherical pendulum is presented most often in literature using Euler angles:
\(\displaystyle{\qquad H\left( {p,q} \right) = \frac{{p_2^2}}{{2m{r^2}\sin {{\left( {q_1} \right)}^2}}} + \frac{{p_1^2}}{{2m{r^2}}} + mgr\cos \left( {q_1} \right) }\)
In the following section, the Hamiltonian of the spherical pendulum will be derived using stereographic projection coordinates (SPC) yielding the Hamiltonian:
\(\displaystyle{\qquad H\left( {p,q} \right) = \frac{{p_1^2 + p_2^2}}{{2m{{\left( {\frac{{2r}}{{q_1^2 + q_2^2 + 1}}} \right)}^2}}} + gm\left( {r - \frac{{2r}}{{q_1^2 + q_2^2 + 1}}} \right) }\)
which is numerically far superior to the acclaimed Euler angle Hamiltonian as is discussed in the last section.#### Defining the spherical pendulum

Define first the variables
#### Deriving the Hamiltonian

The Cartesian coordinates (x,y,z) of the pendulum bob are retrieved directly from SPC, where the canonical coordinates q_{1} and q_{2} are defined as X and Z respectively.
\(\displaystyle{\qquad x = r\frac{{2{q_1}}}{{q_1^2 + q_2^2 + 1}} }\)
\(\displaystyle{\qquad y = r\frac{{ - 1 + q_1^2 + q_2^2}}{{1 + q_1^2 + q_2^2}} }\)
\(\displaystyle{\qquad z = r\frac{{2{q_2}}}{{q_1^2 + q_2^2 + 1}} }\)
where (q_{1},q_{2}) is the projected position of the pendulum bob from the north pole through the pendulum bob and onto the xz-plane.
Note that this projection is ill-defined at the north pole where both q_{1} and q_{2} approaches ∞ with sign that is also undefined.
Hence numerical simulation near the north pole will produce increased error from compactness, significant truncation errors from coordinate scaling and a discontinuity under certain conditions.
Noting that q_{1} and q_{2} are functions of time, the Cartesian velocities are obtained directly as the time derivatives.
The time derivatives can be simplified to:
\(\displaystyle{\qquad \dot x = {{\dot q}_1}a\left( {1 - \frac{{aq_1^2}}{r}} \right) - {{\dot q}_2}\frac{{{q_1}{q_2}{a^2}}}{r} }\)
\(\displaystyle{\qquad \dot y = \frac{{{a^2}}}{r}\left( {{q_1}{{\dot q}_1} + {q_2}{{\dot q}_2}} \right) }\)
\(\displaystyle{\qquad \dot z = {{\dot q}_2}a\left( {1 - \frac{{aq_2^2}}{r}} \right) - {{\dot q}_1}\frac{{{q_1}{q_2}{a^2}}}{r} }\)
where **a** is introduced for the sake of notational simplicity:
\(\displaystyle{\qquad a = \frac{{2r}}{{1 + q_1^2 + q_2^2}} }\)
the kinetic and potential energy are obtained:
\(\displaystyle{\qquad T = \frac{1}{2}m\left( {{{\dot x}^2} + {{\dot y}^2} + {{\dot z}^2}} \right) = \frac{1}{2}m\left( {\dot q_1^2 + \dot q_2^2} \right){a^2} }\)
\(\displaystyle{\qquad V = mgy = mg\left( {r - a} \right) }\)
from which the Lagrangian follows:
\(\displaystyle{\qquad L = T - V = \frac{1}{2}m\left( {\dot q_1^2 + \dot q_2^2} \right){a^2} - mg\left( {r - a} \right) }\)
yielding the canonical momenta:
\(\displaystyle{\qquad \frac{{\partial L}}{{{{\partial \dot q}_1}}} = p_1 = m{\dot q_1}{a^2} }\)
\(\displaystyle{\qquad \frac{{\partial L}}{{{{\partial \dot q}_2}}} = p_2 = m{\dot q_2}{a^2} }\)
The Hamiltonian can now by substitution be expressed explicitly in terms of the canonical coordinates p and q:
\(\displaystyle{\qquad H\left( {p,q} \right) = T + V = \frac{1}{{2m{a^2}}}\left( {p_1^2 + p_2^2} \right) + mg\left( {r - a} \right) }\)
Note that the Hamiltonian is nonseparable as it can not be separated into two functions T(p) and V(q).
It is also worth noting that the Hamiltonian diverges at the north pole where **a** approaches 0 as either q_{1} or q_{2} approaches ∞
The system is simulated using the gradient of the Hamiltonian:
\(\displaystyle{\qquad \frac{{\partial H}}{{\partial p_1}} = \frac{p_1}{m a^2} }\)
\(\displaystyle{\qquad \frac{{\partial H}}{{\partial p_2}} = \frac{p_2}{m a^2} }\)
\(\displaystyle{\qquad \frac{{\partial H}}{{\partial q_1}} = q_1 b }\)
\(\displaystyle{\qquad \frac{{\partial H}}{{\partial q_2}} = q_2 b }\)
where
\(\displaystyle{\qquad b = \frac{{p_1^2 + p_2^2}}{{mra}} + \frac{{mg{a^2}}}{r} }\)
#### Implementation

The following code simulates the spherical pendulum.
The pendulum is displayed in function DrawFunc by transforming the SPC to Cartesian coordinates.
The pendulum bob trajectory is likewise displayed by function IterFunc.
#### Discussion of stereographic projection for Hamiltonian systems

SPC provides a powerful way to derive and simulate Hamiltonian systems associated with spherical and rigid body motion.
The spherical pendulum is the simplest SPC system among the example codes provided with Zymplectic.
The major advantages and the few disadvantages of SPC can be summarized:
#### Remarks

It is remarkable that neither the Hamiltonian or Lagrangian of the spherical pendulum in SPC are to be found in literature, even when the resulting Hamiltonian is simpler, without critical singularity and is computationally more efficient for simulations.

- r as the constant length of the pendulum from the pendulum pivot to the bob
- m as the mass of the bob
- g is the gravitational constant

```
//Copy-paste into Zymplectic or load as file to run the simulation
const double m = 1.0; //mass of the pendulum bob
const double r = 0.5; //pendulum length
const double g = 1.0; //gravitational constant
double H(Arg) {
double a = 2*r/(q1*q1 + q2*q2 + 1);
return (p1*p1 + p2*p2)/(2*m*a*a) + g*m*(r - a);
}
void dH(Arg) {
double a = 2*r/(q1*q1 + q2*q2 + 1);
double b = (p1*p1 + p2*p2)/(m*r*a) + g*m*a*a/r;
dq1 = q1*b;
dq2 = q2*b;
dp1 = p1/(m*a*a);
dp2 = p2/(m*a*a);
}
void DrawFunc() {
double X = obj_q[0][0]; //X,Y is the projected position of the pendulum bob on the xz-plane
double Z = obj_q[0][1];
double a = 2*r/(X*X + Z*Z + 1);
double x = a*X; //x,y,z is the actual position of the pendulum bob
double y = 2*r - a;
double z = a*Z;
draw.rgb(0.1,0.8,0.1);
draw.sphere(x, y, z, 0.05);
draw.rod(x, y, z, 0, r, 0, 0.015);
draw.rgb(0.1,0.1,0.1); draw.rod(x, 0, z, x, -0.01, z, 0.05); //projection onto z-plane
draw.rgb(0.6,0.3,0.6); draw.rod(X, 0.01, Z, X, -0.01, Z, 0.05); //stereographic projection position
draw.rgb(1,0.1,0.1,0.5); draw.rod(X, 0, Z, 0, 2*r, 0, 0.01);
draw.rgb(0.6,0.3,0.6,0.5); draw.sphere(0,r,0,r);
}
void IterFunc() {
double X = obj_q[0][0];
double Z = obj_q[0][1];
double a = 2*r/(X*X + Z*Z + 1);
double x = a*X;
double y = 2*r - a;
double z = a*Z;
draw.track_add(x, y); //display pendulum bob trajectory
draw.trail(x, y, z);
}
double _q_init[] = {1, 0}; //initial position in the stereographic projection plane
double _p_init[] = {0, 0}; //initial momentum is set in main function
void main() {
double dqx = 0; //(dqx,dqy) is the velocity in the stereographic projection plane
double dqy = 1;
double a = 2*r/(_q_init[0]*_q_init[0] + _q_init[1]*_q_init[1] + 1);
_p_init[0] = dqx*m*a*a;
_p_init[1] = dqy*m*a*a;
Z.track = Z.trail = Z.graph3D = 1;
Z.OnDraw = DrawFunc; //draw the pendulum
Z.OnIter = IterFunc; //obtain bob trajectory
Z.H(H,dH,0);
}
```

- The spherical pendulum in Euler angles can be simulated as a separable system unlike the SPC Hamiltonian that is strictly nonseparable. This means that the SPC simulation requires roughly twice as many evaluations.
- SPC has only one singularity, making it more feasible for simulating spherical systems than Euler angles that generally yield two singularities
- The SPC Hamiltonian of the spherical pendulum is rational unlike the Hamiltonian with Euler angles. This enables exact analytical and rational simulations of the spherical pendulum in phase space
- SPC is considerably more efficient for simulating the spherical pendulum as computationally expensive trigonometric functions are avoided. In fact, Zymplectic simulates the SPC roughly three times faster than the Euler angle spherical pendulum.
- The SPC Hamiltonian has low compactness at the south pole. This means that simulations where the pendulum is generally located near the south pole or just the southern hemisphere, generally produce Hamiltonian errors that are easily an order of magnitude smaller than the same simulation in Euler angles.
- The derivation of the Hamiltonian is generally equally difficult as with other coordinates. While the derivation is simple for the spherical pendulum, the complexity increases greatly for systems with more complicated expressions for the canonical momenta, even when the resulting Hamiltonian may be greatly simplified.