In [1]:
using Plots

Symplectic Integrator¶

Standard methods of solving ODE, e.g. explicit Euler's method or RK4, suffer from the problem of "energy drift, i.e. energy is not conserved due to finite step size, not good for large t evolution.

Symplectic integrator adjusts the Hamiltonian in such a way that energy is conserved.

We study a simple harmonic oscillator as illustration:

$$ H(p, q) = \frac{1}{2} p^2 + \frac{1}{2} q^2 $$

References¶

Boudreau & Swanson pp. 383

Stickler & Schachinger pp. 59

H. Yoshida: Symplectic Integrators for Hamiltonian Systems: Basic Theory

H. Yoshida: Construction of higher order symplectic integrators

In [2]:
N = 50
dt = 1.8 * 2*pi / N

Ham = (q, p) -> 0.5 * (q^2 + p^2)

qv = [1.0]
pv = [0.0]
Ev = [Ham(qv[1], pv[1])]

for i1 in 1:N

    q1 = qv[end]
    p1 = pv[end]

    # explicit Euler
    q2 = q1 + p1 * dt
    p2 = p1 - q1 * dt

    push!(qv, q2)
    push!(pv, p2)
    push!(Ev, Ham(q2, p2))
    
end

p1 = scatter(qv, pv, label="Explicit Euler", xlabel="q", ylabel="p")
dqv = qv[2:end] .- qv[1:end-1]
dpv = pv[2:end] .- pv[1:end-1]
quiver!(qv[2:end], pv[2:end], quiver=(dqv, dpv))
p2 = plot(1:N+1, Ev, label="energy", xlabel="n")

plot(p1, p2)
Out[2]:
In [3]:
function RK4_solver(ff, tv, yinit)
    # return the full solution yv(tv, [y, y'])

    Nt = length(tv)
    yv = zeros(Nt, length(yinit))

    yv[1, :] = yinit

    for i1 = 1:Nt-1

        tnow = tv[i1]
        dt = tv[i1+1]-tv[i1]

        ynow = yv[i1, :]

        F1 = dt * ff(tnow, ynow)
        F2 = dt * ff(tnow + dt / 2.0, ynow + F1 / 2.0)
        F3 = dt * ff(tnow + dt / 2.0, ynow + F2 / 2.0)
        F4 = dt * ff(tnow + dt, ynow + F3)

        # RK 4th order
        yv[i1+1, :] = yv[i1, :] + (F1 + 2 * F2 + 2 * F3 + F4) / 6

        # 3/8 rules
        # yv[i1+1, :] = ynow + (F1 + 3*F2 + 3*F3 + F4)/8
    end

    return yv
end

# simple 1D problem
# H = 1/2 * q^2 + 1/2 * p^2
function ff(t, y)
    ret = zeros(2)
    ret[1] = y[2]
    ret[2] = -y[1]    
    return ret
end

N = 25
dt = 2.8 * 2*pi / N
tv = [(i1-1) * dt for i1 in 1:N]
Ham = (q, p) -> 0.5 * (q^2 + p^2)

_qv = RK4_solver(ff, tv, [1.0, 0.0])

qv = _qv[:, 1]
pv = _qv[:, 2]
Ev =[Ham(qv[i1], pv[i1]) for i1 in 1:N] 

p1 = scatter(qv, pv, label="Explicit Euler", xlabel="q", ylabel="p")
dqv = qv[2:end] .- qv[1:end-1]
dpv = pv[2:end] .- pv[1:end-1]
quiver!(qv[2:end], pv[2:end], quiver=(dqv, dpv))
p2 = plot(1:N, Ev, label="energy", xlabel="n")

plot(p1, p2)
Out[3]:

Implementation Notes¶

qdot = dH/dp pdot = -dH/dq

assume separable H(p,q) = T(p) + V(q);

then:

qdot(p)

pdot(q)

order 1: Euler-Cromer¶

first, update q1

q2 = q1 + qdot(p1) * dt

then update p1 with updated q2

p2 = p1 + pdot(q2) * dt

order 2: Verlet¶

q2 = q1 + qdot(p1) * 0.5 * dt

p2 = p1 + pdot(q2) * dt

then update q2 again

q2 = q2 + qdot(p2) * 0.5 * dt

In [4]:
N = 40
dt = 1.6 * 2*pi/N

# for this simple case
qdot = (q, p) -> p
pdot = (q, p) -> -1 * q

Ham = (q, p) -> 0.5 * (q^2 + p^2)

qv = [1.0]
pv = [0.0]
Ev = [Ham(qv[1], pv[1])]

for i1 in 1:N

    q1 = qv[end]
    p1 = pv[end]

    q2 = q1 + p1*dt
    p2 = p1 - dt*q2

    push!(qv, q2)
    push!(pv, p2)
    push!(Ev, Ham(q2, p2))
    
end

p1 = plot(qv, pv, label="symplectic")


dqv = qv[2:end] .- qv[1:end-1]
dpv = pv[2:end] .- pv[1:end-1]

quiver!(
    qv[2:end], pv[2:end],
    quiver=(dqv,dpv)
)

p2 = plot(1:N+1, Ev, label="energy", xlabel="n")

plot(p1, p2)
Out[4]:
In [5]:
N = 40
dt = 1.6 * 2*pi/N

# for this simple case
qdot = (q, p) -> p
pdot = (q, p) -> -1.0 * q

Ham = (q, p) -> 0.5 * (q^2 + p^2)

qv = [1.0]
pv = [0.0]
Ev = [Ham(qv[1], pv[1])]

for i1 in 1:N

    q1 = qv[end]
    p1 = pv[end]

    # Verlet method
    q2 = q1 + 0.5*dt * p1
    p2 = p1 - dt*q2
    q2 = q2 + 0.5*dt * p2

    push!(qv, q2)
    push!(pv, p2)
    push!(Ev, Ham(q2, p2))
    
end

p1 = plot(qv, pv, label="symplectic")


dqv = qv[2:end] .- qv[1:end-1]
dpv = pv[2:end] .- pv[1:end-1]

quiver!(
    qv[2:end], pv[2:end],
    quiver=(dqv,dpv)
)

p2 = plot(1:N+1, Ev, label="energy", xlabel="n")

plot(p1, p2)
Out[5]: