We consider here the simplest case of explicit symplectic integrators where the Hamiltonian is separable and autonomous:
$$H(\mathbf{p},\mathbf{q}) = T(\mathbf{p}) + V(\mathbf{q})$$
The implementations described in this section can be applied to a wide range of physical systems, but are here shown specifically for the Henon-Heiles system, which is expressed by:
$$T(\mathbf{p}) = \frac{1}{2}(p_1^2 + p_2^2)$$
$$V(\mathbf{q}) = \frac{1}{2}(q_1^2 + q_2^2) + \lambda(q_1^2 q_2 -\frac{q_2^3}{3})$$
where λ is here set to 1. The potential is visualized and first part of the simulation is shown below:
#### Symplectic euler (method 0)

As this method is not particularly accurate or useful, we list this method only as background of the methods to follow.
The symplectic time evolution of the coordinates
$$(q_{old},p_{old}) \rightarrow (q_{new},p_{new})$$
over the time step τ can be obtained by the symplectic Euler method:
$$q_{new}=q_{old} + \tau\frac{\partial T}{\partial p}(p_{old})$$
$$p_{new}=p_{old} - \tau\frac{\partial V}{\partial q}(q_{new})$$
where q and p are interpreted as vectors in the multidimensional case. This procedure can be repeated any number of iterations without generating any dissipative error of the Hamiltonian even if the symplectic Euler method is only of first order. It should be noted that the time evolution is only symplectic if τ is reasonably small and the second derivatives of the Hamiltonian are smooth in the relevant area of integration.
#### Vanilla explicit symplectic integrator (method 1)

Higher order explicit symplectic integrators are obtained by splitting the time evolution by multiple intermediate steps:
$$q_{1}=q_{old} + \tau c_1\frac{\partial T}{\partial p}(p_{old}) $$
$$p_{1}=p_{old} - \tau d_1\frac{\partial V}{\partial q}(q_{1})$$
$$q_{2}=q_{1} + \tau c_2\frac{\partial T}{\partial p}(p_{1}) $$
$$p_{2}=p_{1} - \tau d_2\frac{\partial V}{\partial q}(q_{2})$$
$$q_{3}=q_{2} + \tau c_3\frac{\partial T}{\partial p}(p_{2}) $$
$$p_{3}=p_{2} - \tau d_3\frac{\partial V}{\partial q}(q_{3})$$
$$...$$
$$q_{new}=q_{s-1} + \tau c_{s}\frac{\partial T}{\partial p}(p_{s-1}) $$
$$p_{new}=p_{s-1} - \tau d_{s}\frac{\partial V}{\partial q}(q_{s})$$
where c and d are arrays containing the splitting coefficients which each sum up to 1, and s is the number of "stages" associated with the length of these arrays. The symplectic Euler method is considered a 1-stage integrator where c_{1}=d_{1}=1. It is important to note that only q_{new}, p_{new}, q_{old}, p_{old} are meaningful and that the intermediate steps provide no useful information about the time evolution. The algorithm is implemented in the following MATLAB code, which simulates the Henon-Heiles system and displays the Hamiltonian error as a function of time.

```
% We consider the two-dimensional Hamiltonian system "Henon-Heiles"
% The Hamiltonian is separable and autonomous
lambda = 1.0;
Hamilton = @(p,q) (0.5*(p(1)^2+p(2)^2)+0.5*(q(1)^2+q(2)^2) + lambda*(q(1)^2*q(2)-q(2)^3/3));
dV = @(q) ([q(1)+2*lambda*q(1)*q(2), q(2)+lambda*(q(1)^2-q(2)^2)]); %Gradient of V
dT = @(p) ([p(1), p(2)]); %Gradient of P
% Integrator coefficients can be obtained from the text file "integrators.zym"
% The following coefficients "Ruth 1989" provide a 4-stage 4th order integrator
c(1) = 0.675603595979828817;
c(2) =-0.175603595979828817;
c(3) =-0.175603595979828817;
c(4) = 0.675603595979828817;
d(1) = 1.351207191959657634;
d(2) =-1.702414383919315268;
d(3) = 1.351207191959657634;
d(4) = 0;
dt = 0.01; %time step size also denoted tau or h in other literature
IterCount = 5000; %number of iterations
H = zeros(IterCount,1);
T = dt*(0:(IterCount-1));
q = [0.3 0.0];
p = [0.0 0.4];
H(1) = Hamilton(p(1,:),q(1,:)); %Initial Hamiltonian
for n = 1:IterCount
for i = 1:length(c)
q = q + c(i)*dt*dT(p);
p = p - d(i)*dt*dV(q);
end
H(n) = Hamilton(p,q);
end
plot(T,H - H(1)) %display graphically the error of the Hamiltonian
```

```
lambda = 1.0;
Hamilton = @(p,q) (0.5*(p(1)^2+p(2)^2)+0.5*(q(1)^2+q(2)^2) + lambda*(q(1)^2*q(2)-q(2)^3/3));
dV = @(q) ([q(1)+2*lambda*q(1)*q(2), q(2)+lambda*(q(1)^2-q(2)^2)]);
dT = @(p) ([p(1), p(2)]);
c(1) = 0.675603595979828817; %"Ruth 1989"
c(2) =-0.175603595979828817;
c(3) =-0.175603595979828817;
c(4) = 0.675603595979828817;
d(1) = 1.351207191959657634;
d(2) =-1.702414383919315268;
d(3) = 1.351207191959657634; %c(4) is now undeclared and unused
dt = 0.01;
IterCount = 5000;
H = zeros(IterCount,1);
T = dt*(0:(IterCount-1));
q = [0.3 0.0];
p = [0.0 0.4];
H(1) = Hamilton(p(1,:),q(1,:));
FSAL = dT(p); %First same as last
for n = 1:IterCount
q = q + c(1)*dt*FSAL; %Calculated in the end of each step
p = p - d(1)*dt*dV(q);
for i = 2:(length(c)-1)
q = q + c(i)*dt*dT(p);
p = p - d(i)*dt*dV(q);
end
FSAL = dT(p); %Used at the beginning of next step
q = q + c(length(c))*dt*FSAL;
H(n) = Hamilton(p,q);
end
plot(T,H - H(1))
```

```
lambda = 1.0;
Hamilton = @(p,q) (0.5*(p(1)^2+p(2)^2)+0.5*(q(1)^2+q(2)^2) + lambda*(q(1)^2*q(2)-q(2)^3/3));
dV = @(q) ([q(1)+2*lambda*q(1)*q(2), q(2)+lambda*(q(1)^2-q(2)^2)]);
dT = @(p) ([p(1), p(2)]);
% The following coefficients provide a 3-stage 3rd order integrator
c(1) = 1;
c(2) =-2/3;
c(3) = 2/3;
d(1) =-1/24;
d(2) = 3/4;
d(3) = 7/24;
dt = 0.01;
IterCount = 5000;
H = zeros(IterCount,1);
T = dt*(0:(IterCount-1));
q = [0.3 0.0];
p = [0.0 0.4];
H(1) = Hamilton(p(1,:),q(1,:)); %Initial Hamiltonian
for n = 1:IterCount
if (mod(n,2) == 0) %Alternating between forward and backward implementation
for i = 1:length(d)
q = q + c(i)*dt*dT(p);
p = p - d(i)*dt*dV(q);
end
else
for i = length(d):-1:1
p = p - d(i)*dt*dV(q);
q = q + c(i)*dt*dT(p);
end
end
H(n) = Hamilton(p,q);
end
plot(T,H - H(1))
```

```
lambda = 1.0;
Hamilton = @(p,q) (0.5*(p(1)^2+p(2)^2)+0.5*(q(1)^2+q(2)^2) + lambda*(q(1)^2*q(2)-q(2)^3/3));
dV = @(q) ([q(1)+2*lambda*q(1)*q(2), q(2)+lambda*(q(1)^2-q(2)^2)]);
dT = @(p) ([p(1), p(2)]);
% The following coefficients provide a 3-stage 3rd order integrator
c(1) = 1;
c(2) =-2/3;
c(3) = 2/3;
d(1) =-1/24;
d(2) = 3/4;
d(3) = 7/24;
dt = 0.01;
IterCount = 5000;
H = zeros(IterCount,1);
T = dt*(0:(IterCount-1));
q = [0.3 0.0];
p = [0.0 0.4];
H(1) = Hamilton(p(1,:),q(1,:)); %Initial Hamiltonian
FSAL = dT(p);
for n = 1:IterCount
if (mod(n,2) == 0) %Alternating between forward and backward implementation
q = q + c(1)*dt*FSAL;
p = p - d(1)*dt*dV(q);
for i = 2:(length(d)-1)
q = q + c(i)*dt*dT(p);
p = p - d(i)*dt*dV(q);
end
q = q + c(length(d))*dt*dT(p);
FSAL = dV(q);
p = p - d(length(d))*dt*FSAL;
else
p = p - d(length(d))*dt*FSAL;
q = q + c(length(d))*dt*dT(p);
for i = (length(d)-1):-1:2
p = p - d(i)*dt*dV(q);
q = q + c(i)*dt*dT(p);
end
p = p - d(1)*dt*dV(q);
FSAL = dT(p);
q = q + c(1)*dt*FSAL;
end
H(n) = Hamilton(p,q);
end
plot(T,H - H(1))
```

When comparing the computational cost of method 2 and 4, it is found that the number of gradient evaluations is proportional to the total number of coefficients subtracted by 2. It is left to the reader to obtain an optimal method for integrators where d(last)=0, although it should be noted that symmetric coefficients implemented with method 2 and 4 provide exactly the same performance and at the same cost.
It is also worth noting that the optimized variants discussed here may not perform notably faster due to increased overhead time. It is left to the reader to implement the methods with emphasis on performance over readability.