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

00ODE.en.md

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:

• harmonic oscillator
• two body problem

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

acc.jl

#!/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[2]/m,
     (t,x) -> -x[1]*(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")

LeapFrog.jl

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[1] = 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

lf3.jl

#!/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[2]/m,
     (t,x) -> -x[1]*(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")

pv-0.4-10000.pdf

rk-0.4-10000.pdf

rk1.jl

#!/usr/bin/env julia
# rk1.jl
# 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

# H = p^2/(2m) + k*q^2/2
m = 1.0
k = 1.0
f = [(t,x) ->  x[2]/m,
     (t,x) -> -x[1]*k]   # array of anonymous functions
x0 = [0.0, 2.0]
t0 = 0.0
h = 0.01
n = 1000

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[:])
savefig("rk1.eps")

style.css

/* -*-CSS-*-
 * style.css for README.html of feram
 * Time-stamp: <2013-11-30 12:19:00 takeshi>
 * Author: Takeshi NISHIMATSU
 */
body {
  color:       black;
  font-family: verdana, arial, helvetica, sans-serif;
}

h1, h2, h3, h4, h6 {
  font-family: verdana, arial, helvetica, sans-serif;
}

h1 {
  color:            #dd0000;
  background-color: #fff0f0;
  font-size:        240%;
}

h2 {
  border-top:       red 5px solid;
  border-bottom:    red 1px solid;
  padding-left:     8px;
  background-color: #fff0f0;
}

h3 {
  border-top:    red 2px solid;
  border-bottom: red 1px solid;
  padding-left:  4px;
}

h4 {
  border-top:    red 1px solid;
  padding-left:  4px;
  background-color: #fff0f0;
}

h5 {
  font-size:     larger;
  font-family:   courier, verdana, arial, helvetica, sans-serif;
  padding-top:   10px;
  color:         darkred;
}

pre {
  font-family:      monospace, courier, verdana, arial, helvetica, sans-serif;
  padding-right:    0.5em;
  padding-left:     0.5em;
  padding-top:      0.1ex;
  padding-bottom:   0.1ex;
  margin-left:      0.5em;
  margin-right:     1.0em;
  white-space:      pre;
  color:            darkred;
  background-color: #f3f3f3;
}

div.figure img {
    width:50%;
    margin: auto;
    display: block;
}

div.figure div.figcaption {
    width: 60%;
    margin: auto;
    display: block;
}

div.navi {
 text-align: right;
 margin-right: 1.0em;
}

div.contents {
    margin-left: 10%;
}

figure img{
    width: 50%;
    margin: auto;
    margin-top:     3.0em;
    display: block;
}

figure figcaption{
    width: 60%;
    margin: auto;
    margin-bottom:     3.0em;
    display: block;
}

table {
    border: blue 2px solid;
    text-align: center;
    margin: auto;
}