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

numerical methodsordinary 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).

Best Answer

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

plot of the orbits

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

enter image description here

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
Related Question