example for symplectic integration with ODE.verlet : pendulum

example_pendulum.jl

# symplectic integration
## pedulum
using Winston
using ODE

# rhs
function pendulum(t, q)
	return -sin(q)
end

# dim.less energy: epsilon = v^2/2 + 1 - cos(x)
# epsilon < 2: oscillations
# epsilon > 2: rotations
epsilon = 1.
v0 = 0.
x0 = acos(1. + v0^2/2 - epsilon)


t, x, v = ODE.verlet_fixed(pendulum, linspace(0.,4pi,11), x0, v0)

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/2 + 1. - cos(x), color="blue"))
add(p1, Curve([0, 2pi], [v0^2/2 + 1 - cos(x0), v0^2/2 + 1 - cos(x0)], color="black"))

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

# adaptive time-step integration
t, x, v = ODE.verlet_hh2(pendulum, [0.,4pi], x0,v0)
add(p1, Points( t, v.*v/2 + 1. - cos(x), color="red"))
add(p2, Points( mod(x+pi,2pi)-pi, v, color="red"))
alpha = linspace(0.,2pi,101)

display(tab)