# Energy-preserving numerical method for system of coupled 2nd order differential equations

I'm a physicist who, due to Covid-19, switched to programming simulations. I would like to know if there is a numerical method that preserves the energy of the system for a second order coupled differential equation of the form:

\begin{aligned} \ddot x &= a x + b \dot y\\ \ddot y &= a y - b \dot x \end{aligned}

where $$a$$ and $$b$$ are real constants.

I'm searching for something equivalent to the leapfrog integration, which seems to work only for non-coupled systems where $$\ddot x$$ does not depend on $$\dot x$$. The goal is to preserve the Hamiltonian of the system over time (i.e., constrain its error).

• Are $a,b$ constants ? If yes, why don't you use an analytical solution ?
– user65203
Jun 25, 2020 at 15:40
• That can indeed be solved explicitly, you just need to diagonalize $\begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ a & 0 & 0 & b \\ 0 & a & -b & 0 \end{bmatrix}$, which looks like it can be done without any horrible quartic formula.
– Ian
Jun 25, 2020 at 15:41
• Yes, $a$ and $b$ are real constants. I don't use the analytical solution because this is only the "base" system where additional forces will be added (which partially have random time behavior, but are not a function of $x$, $y$ or the derivatives). Jun 25, 2020 at 15:43
• Why not using a Runge-Kutta integrator ?
– user65203
Jun 25, 2020 at 15:47
• If you want that, you should look into symplectic integration, for example the Verlet method...but when you start adding non-conservative inputs, energy conservation will not perform as well as it does when the underlying system conserves energy.
– Ian
Jun 25, 2020 at 15:51

The given system, as a special case of the Lorentz system of a charged particle in an electromagnetic field, has in fact a perfectly serviceable Hamiltonian $$H=\frac12[(p+\tilde by)^2+(q-\tilde bx)^2]-\frac a2(x^2+y^2).$$ where $$(p,q)$$ is the impulse vector to the position vector $$(x,y)$$. Formulating any symplectic method for the resulting Hamiltonian system \begin{align} \dot x=H_p&=p+\tilde by, \\ \dot y=H_q&=q-\tilde bx, \\ \dot p=-H_x&=-\tilde b(\tilde bx-q)+ax, \\ \dot q = -H_y &= -\tilde b(\tilde by+p)+ay, \end{align} results in the conservation of a modified Hamiltonian as energy functional. In other words, while the error in the energy still has the order of the method, it is completely dependent on the state. So if the trajectory is bounded in position and impulse, the energy error will be likewise bounded. If the trajectory is periodic or quasi-periodic, the same applies to the error.

Unlike the usual mechanical systems with a separable Hamiltonian, the step equations will be all implicit here, that means that one has to solve some linear systems.

Remark: To test that the Hamiltonian is correct, compute the second derivative and eliminate the impulses \begin{align} \ddot x &= \dot p + \tilde b\dot y = ax+2\tilde b\dot y \\ \ddot y &= \dot y- \tilde b\dot x = ay-2\tilde b\dot x \end{align} so that $$\tilde b=\frac b2$$.

For instance, the Verlet method step of step size (for convenience) $$h=2\Delta t$$ consists of two opposite symplectic Euler steps of half the step size $$Δt$$, usually first a space-explicit, impulse-implicit step and then a space-implicit, impulse-explicit step \begin{align} \vec x_{n+1/2}&=\vec x_n + H_{\vec p}(\vec x_n,\vec p_{n+1/2})Δt \\ \vec p_{n+1/2}&=\vec p_n - H_{\vec x}(\vec x_n,\vec p_{n+1/2})Δt \\ \hline \vec x_{n+1}&=\vec x_{n+1/2} + H_{\vec p}(\vec x_{n+1},\vec p_{n+1/2})Δt \\ \vec p_{n+1}&=\vec p_{n+1/2} - H_{\vec x}(\vec x_{n+1},\vec p_{n+1/2})Δt \end{align}

def semi_Euler_A(u,dt):
x0,y0,p0,q0 = u
# p1 = p0 + dt*( b*(q1-b*x0)+a*x0)
# q1 = q0 + dt*(-b*(p1+b*y0)+a*y0)
p1,q1 = solve([[1,-dt*b],[dt*b,1]], [p0+dt*(a-b*b)*x0, q0+dt*(a-b*b)*y0])
x1 = x0 + dt*(p1+b*y0)
y1 = y0 + dt*(q1-b*x0)
return x1,y1,p1,q1

def semi_Euler_B(u,dt):
x1,y1,p1,q1 = u
# x2 = x1 + dt*(p1+b*y2)
# y2 = y1 + dt*(q1-b*x2)
x2,y2 = solve([[1,-dt*b],[dt*b,1]],[x1+dt*p1, y1+dt*q1])
p2 = p1 + dt*( b*(q1-b*x2)+a*x2)
q2 = q1 + dt*(-b*(p1+b*y2)+a*y2)
return x2,y2,p2,q2

def Verlet_step(u0,h):
u1 = semi_Euler_A(u0,h/2)
u2 = semi_Euler_B(u1,h/2)
return u2


To get bounded energy levels and thus bounded trajectories select $$a$$ negative, for the following plots I chose $$a=-4$$, $$\tilde b=1$$, $$(x_0,y_0)=(0,2)$$, $$(\dot x_0,\dot y_0)=(1,0)$$. The step size is $$h=2\Delta t=0.2$$ (in the above sense, $$h$$ for the RK4 step size). Then the numerical trajectories in the $$(x,y)$$ plane for $$t\in[0,10]$$ are

One can see that the 4th order RK4 method is closer to the much more precise odeint solution than the 2nd order Verlet method. For the energy extend the interval to $$t\in[0,30]$$ to get the differences to the initial energy as

As one can see, the initial error of the Verlet method is larger. But as the Verlet error terms in the energy expansion depend on the state, not on the time (except for the drift due to accumulation of floating point noise), the Verlet error is periodic and thus stays bounded as the solution stays bounded. In contrast the RK4 energy error shows a significant slope and increases steadily. Of course, for smaller step sizes the cross-over point will be later, as the oscillation in the Verlet error is $$O(h^2)$$ while the RK4 error is $$O(th^4)$$ (for medium large $$t$$, for huge $$t$$ it becomes exponential in $$t$$). Empirically, the RK4 error is even about $$70th^5$$, the Verlet error $$9.5h^2$$, so that the cross-over point is at about $$t=0.14/h^3$$.

One can raise the order of the Verlet energy error (and of the whole method) to 4 by Ruth-Forest-Yoshida, which is just an extrapolation eliminating the quadratic error terms, that is, the cubic local error terms, in a sequence of Verlet steps of length $$b_0h,b_1h,b_0h$$. This preserves the time symmetry. The parameters have to satisfy $$2b_0+b_1=1$$ and $$2b_0^3+b_1^3=0$$ leading to the modified main loop

    b0 = 1/(2-2**(1/3)) # =  1.35120719196
b1 = 1-2*b0         # = -1.70241438392
for i in range(len(t)-1):
u[i+1] = Verlet_step(u[i],b0*h)
u[i+1] = Verlet_step(u[i+1],b1*h)
u[i+1] = Verlet_step(u[i+1],b0*h)
return t,u

• If I recall correctly, symplectic methods preserve the symplectic form, not the Hamiltonian. Jun 27, 2020 at 10:29
• If you want to get into details, symplectic methods preserve a symplectic form that is related to the form of the continuous problem. Also, one can find a modified Hamiltonian for fixed step size that is preserved exactly. For periodic solutions the deviation to the exact Hamiltonian is also periodic. Jun 27, 2020 at 10:38
• Correct, but from a physical.perspective (without differential geometry, so to speak) it may be misleading to present symplectic methods as "energy preserving". Jun 27, 2020 at 10:45
• @Miguel : Now the formulation should be better, supported by the example. Jun 29, 2020 at 6:30