Skip to content

Instantly share code, notes, and snippets.

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

Embed
What would you like to do?
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[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")
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
#!/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")
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
#!/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")
/* -*-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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.