Instantly share code, notes, and snippets.

# t-nissie/00ODE.en.md Last active Nov 13, 2016

Learn ordinary differential equation (ODE) solvers with a harmonic oscillator

# Learn ordinary differential equation (ODE) solvers with a harmonic oscillator

Solve some ordinary differential equations (ODE).

Winston is used for plotting, but `using Winston` takes some time. Github: https://github.com/nolta/Winston.jl Documents: https://winston.readthedocs.io/

## Comparison of some geometric integrators

Comparison of some geometric integrators, Euler, symplectic-Euler, velocity-Verlet, position-Verlet and leapfrog methods.

Package: GeometricIntegrators.jl

examples of usage:

Symplectic-Euler, velocity-Verlet, position-Verlet and leapfrog are so-called symplectic integrators. Test them with `h=0.4` and `n=10000`.

## Runge-Kutta method

Runge-Kutta: rk1.jl

### Energy is not conserved

Runge-Kutta method is not a symplectic integrator. Confirm that energy is not coserved with `h=0.4` and `n=10000`.

### Acceleration

Runge-Kutta method is good for simulations of acceleration: acc.jl

## Pretty print

``````\$ LANG=C a2ps --header='Printed by t-nissie' --media=a4 --prologue=color\
--portrait --columns=1 --margin=3 --borders=off -f 11.0 --pretty-print=julia\
-o RungeKutta.ps ~/.julia/v0.4/RungeKutta/src/RungeKutta.jl
\$ ~/scripts/PsDuplex RungeKutta.ps | lpr -l
``````
 #!/usr/bin/env julia # acceleration of an oscillator. # https://github.com/timothyrenner/RungeKutta.jl # Pkg.clone("git://github.com/timothyrenner/RungeKutta.jl.git") # Pkg.add("Winston") or Pkg.clone("git@github.com:nolta/Winston.jl.git") # Install these packages only for the first time. ## using RungeKutta m = 1.0 k0 = 1.0 k1 = 0.1 omega = sqrt(k0/m) f = [(t,x) -> x/m, (t,x) -> -x*(k0+k1*cos(2*omega*t))] x0 = [0.0, 2.0] t0 = 0.0 h = 0.1 n = 400 t,x = rk4f(f, x0, t0, h, n) energy = zeros(Float64, n+1) for i in 1:n+1 energy[i] = x[2,i]^2/2 + x[1,i]^2/2 end using Winston plot(t, x[1,:], t, x[2,:], t, energy[:]) xlim(-1,41) ylim(-4,8) savefig("acc.eps")
 module LeapFrog # Implements a leapfrog solver for systems of ODEs. # Computes the value of the next point in the leapfrog iteration. # # Args: # f: The array of functions for computing x_n+1. # x: The current value of the points in the space. # t: The current value of time. # h: The step size. # # Returns: The next point in the solution. function nextPoint{T<:AbstractFloat}( f::Array{Function,1}, x::Array{T,1}, t::AbstractFloat, h::AbstractFloat) q_half = x + h/2 .* map(f_each -> f_each(t, x) , f) p_next = x + h .* map(f_each -> f_each(t+h/2, q_half), f) q_next = q_half + h/2 .* map(f_each -> f_each(t+h, p_next), f) # t+h/2 ??? return t + h, vcat(q_next[1:length(x)÷2],p_next[length(x)÷2+1:length(x)]) end # Solves a system of ODEs using the leapfrog method. # # Args: # f: An array of functions defining the system of first-order ODEs such that # f[i](t, x_n) = x[i]_{n+1} # x0: The initial point. # t0: The initial time. # h: The step size. # n: The number of iterations. # # Throws: # ArgumentError: If the length of f is not equal to the length of x0. # ArgumentError: If n is not greater than zero. # # Returns: The vector of times, and the vector of points produced by the solver. function leapfrog{T<:AbstractFloat}(f::Array{Function,1}, x0::Array{T,1}, t0::AbstractFloat, h::AbstractFloat, n::Integer) #Validate that the lengths of f and x0 are the same. if length(f) != length(x0) throw(ArgumentError("There must be one function per element of x0.")); end #Validate that n is greater than zero. if n <= 0 throw(ArgumentError("Number of iterations must be positive.")); end #Validate that the step size is greater than zero. if h <= 0.0 throw(ArgumentError("Step size must be positive.")); end x = zeros(length(x0), n+1); t = zeros(1, n+1); #The first points are the initial conditions. x[:,1] = x0; t = t0; for ii in 2:(n+1) t[ii], x[:,ii] = nextPoint(f, x[:,ii-1], t[ii-1], h); end return t,x end export leapfrog end
 #!/usr/bin/env julia # Harmonic oscillator for a ordinary differential equation (ODE) solver ### include("LeapFrog.jl") using LeapFrog m = 1.0 k0 = 1.0 k1 = 0.1 omega = sqrt(k0/m) f = [(t,x) -> x/m, (t,x) -> -x*(k0+k1*cos(2*omega*t))] x0 = [0.0, 2.0] t0 = 0.0 h = 0.1 n = 400 t,x = leapfrog(f, x0, t0, h, n) energy = zeros(Float64, n+1) for i in 1:n+1 energy[i] = x[2,i]^2/2 + x[1,i]^2/2 end using Winston plot(t, x[1,:], t, x[2,:], t, energy[:]) xlim(-1,41) ylim(-4,8) savefig("lf3.eps")  