Symplectic integration of harmonic oscillator 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?
 A: 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)

