example for symplectic integration with ODE.verlet : harmonic oscillator

example_ho.jl

# symplectic integration
## harmonic oscillator
using Winston
using ODE

# rhs
function harmonic_oscillator(t, q)
	return -q
end

# fixed time-step integration
t, x, v = ODE.verlet_fixed(harmonic_oscillator, linspace(0.,2pi,11), 1., 0.)

tab = Table(2, 1)
# upper panel: total energy vs time
p1 = FramedPlot(ylabel="energy", xlabel="time")
tab[1,1] = p1
add(p1, Curve( t, (v.*v + x.*x)/2, color="blue"))
add(p1, Curve([0, 2pi], [(1.^2 + 0.^2)/2, (1.^2 + 0.^2)/2], color="black"))

# lower panel: trajectory in phase-space
p2 = FramedPlot(ylabel="velocity", xlabel="position")
tab[2,1] = p2
add(p2, Points(x, v, color="blue"))

# adaptive time-step integration
t, x, v = ODE.verlet_hh2(harmonic_oscillator, linspace(0.,2pi,11), 1., 0.)
add(p1, Points( t, (v.*v + x.*x)/2, color="red"))

add(p2, Points(x, v, color="red"))
alpha = linspace(0.,2pi,101)
xcirc = cos(alpha)
vcirc = sin(alpha)
add(p2, Curve(xcirc, vcirc, color="black"))

display(tab)