Symplectic integration of harmonic oscillator

Question

I try to get numerical solution of ordinary harmonic oscillator with symplectic integrator. The problem is that what I obtain doesn't conserve energy (but symplectic integration should do).

I consider 2 simplectic integrators - symplectic Euler method and symplectic 4th order method.

The equations of motion are:

$$ \dot{p} = -q \\ \dot{q} = p $$

So for Euler method we have formulas:

$$q_{n+1} = q_n + \Delta t \ p_n \\ p_{n+1} = p_n - \Delta t \ q_{n+1} $$

For 4th order symplectic method one has sequence of intermediate points: $$ p_{n+1/4} = p_n - \frac{2}{2 (2-2^{1/3})} \Delta t \ q_n \\ q_{n+1/4} = q_n + \frac{1}{2 (2-2^{1/3})} \Delta t \ p_{n+1/4} \\ \\ p_{n+2/4} = p_{n+1/4} + \frac{2^{1/3}}{2 - 2^{1/3}} \Delta t \ q_{n+1/4} \\ q_{n+2/4} = q_{n+1/4} + \frac{1-2^{1/3}}{2 (2-2^{1/3})} \Delta t \ p_{n+2/4} \\ \\ p_{n+3/4} = p_{n+2/4} - \frac{2}{2 (2-2^{1/3})} \Delta t \ q_{n+2/4} \\ q_{n+3/4} = q_{n+2/4} + \frac{1-2^{1/3}}{2 (2-2^{1/3})} \Delta t \ p_{n+3/4} \\ \\ p_{n+1} = p_{n+3/4} \\ q_{n+1} = q_{n+3/4} + \frac{1}{2 (2-2^{1/3})} \Delta t \ p_{n+1} $$

But both methods (Euler and 4th order) does not conserve energy. I choose initial conditions as $$ q_0 = 1.2 \ , \ \ p_0 = 0 $$

But what I have with both methods is not oscillations with conserved amplitude - it's oscillations with growing and growing (sufficiently fast) amplitude.

So, what might be a problem?

Answer 1

The energy should oscillate, but not grow. There is a -- for Euler O(h) -- modified energy function that should be constant. Indeed, $$ H(p,q)=p^2+q^2+Δt·pq $$ is constant for the symplectic Euler integration, test with

q0 = 1.2
p0 = 0
E0 = p0**2+q0**2

T = 6.2
N = 20
dt = T/N

q = q0
p = p0
t=0

for k in range(N+1):
    q = q + p*dt
    p = p - q*dt
    t = t + dt
    E = p**2+q**2

    print "%6.3f : p=%12.8f q=%12.8f : E=%12.8f, E-E0+dt*p*q=%12.6e\n" % (t,q,p,E,E-E0+dt*p*q)

And for composite Verlet4, the following test confirms that the energy is constant up to an oscillating term of order $O(Δt^4)$, indeed the modification $$ E(q,p)=p^2+(1 + \tfrac{1}{24}·b_0^2·Δt^4)·q^2 $$ ($b_0=$ b0 as below) is constant up to errors $O(Δt^6)$.

def Verlet(q,p,dt):
    p = p - 0.5*dt*q
    q = q + dt*p
    p = p - 0.5*dt*q
    return q,p

def Verlet4(q,p,dt):
    b0 = b0=1/(2-2**(1./3))
    b1=2*b0-1
    q,p=Verlet(q,p, b0*dt)
    q,p=Verlet(q,p,-b1*dt)
    q,p=Verlet(q,p, b0*dt)
    return q,p

q0 = 1.2
p0 = 0
E0 = p0**2+q0**2

T = 20.0
N = 120
dt = T/N

q = q0
p = p0
t=0

for k in range(N+1):
    q,p = Verlet4(q,p,dt)
    t = t + dt
    E = p**2+q**2

    print "%6.3f : p=%12.8f q=%12.8f : E=%12.8f, (E-E0)/dt^4=%12.6e p*q=%12.6e" % (t,q,p,E,(E-E0)/dt**4,p*q)