Created
January 27, 2014 11:23
-
-
Save acroy/8646927 to your computer and use it in GitHub Desktop.
example for symplectic integration with ODE.verlet : harmonic oscillator
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment