From ad9851207c862086df43e10e21d7ac5854211263 Mon Sep 17 00:00:00 2001 From: acroy Date: Mon, 27 Jan 2014 11:23:08 +0100 Subject: [PATCH 1/4] Added symplectic solver based on velocity verlet. --- src/ODE.jl | 4 +- src/symplectic.jl | 119 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 122 insertions(+), 1 deletion(-) create mode 100644 src/symplectic.jl diff --git a/src/ODE.jl b/src/ODE.jl index c2243831f..f71c2021c 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -5,7 +5,7 @@ module ODE using Polynomial ## minimal function export list -export ode23, ode4, ode45, ode4s, ode4ms +export ode23, ode4, ode45, ode4s, ode4ms, verlet ## complete function export list #export ode23, ode4, @@ -13,6 +13,8 @@ export ode23, ode4, ode45, ode4s, ode4ms # oderosenbrock, ode4s, ode4s_kr, ode4s_s, # ode4ms, ode_ms +include("symplectic.jl") + #ODE23 Solve non-stiff differential equations. # diff --git a/src/symplectic.jl b/src/symplectic.jl new file mode 100644 index 000000000..9daab84d7 --- /dev/null +++ b/src/symplectic.jl @@ -0,0 +1,119 @@ +# symplectic integration methods + +## velocity verlet with fixed step-size +## integrates the system +## dx/dt = v(t) +## dv/dt = a(t, x) = F(t, x)/m +## with initial conditions x_0 and v_0 using the velocity-verlet method. +## +## The function takes as a first argument a function, a(t,x), +## of time and position (the acceleration). The second argument +## specifies a list of times for which a solution is requested. The last +## two arguments are the initial conditions x_0 and v_0. +## +## The output, (t,x,v), consists of the solutions x(t), v(t) at times t=tspan. +## +function verlet_fixed(a, tspan, x_0, v_0) + + if length(x_0) != length(v_0) + error("Intial data x_0 and v_0 must have equal length.") + end + + x = Array(eltype(x_0), length(tspan), length(x_0)) + v = Array(eltype(v_0), length(tspan), length(v_0)) + + x[1,:] = x_0' + v[1,:] = v_0' + + atmp = a(tspan[1], x_0).' + + for i = 1:(length(tspan)-1) + h = tspan[i+1] - tspan[i] + + v[i+1,:] = v[i,:] + atmp*h/2 + x[i+1,:] = x[i,:] + v[i+1,:]*h + atmp = a(tspan[i]+h, x[i+1,:].').' + v[i+1,:] = v[i+1,:] + atmp*h/2 + end + + return tspan, x, v +end + +## velocity verlet with adaptive step-size +## integrates the system +## dx/dt = v(t) +## dv/dt = a(t, x) = F(t, x)/m +## with initial conditions x_0 and v_0 using the velocity-verlet method. +## The step-size is adapted accoording to an error estimate based on +## comparing the results of two half steps and one full step. +## (for some background see: http://www.theorphys.science.ru.nl/people/fasolino/sub_java/pendula/computational-en.shtml) +## +## The function takes as a first argument a function, a(t,x), +## of time and position (the acceleration). The second argument +## specifies a list of times for which a solution is requested. The last +## two arguments are the initial conditions x_0 and v_0. +## +## The output, (t,x,v), consists of the solutions x(t), v(t) at the +## intermediate integration times t including tspan[1] and tspan[end]. +## +function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm) + if length(x_0) != length(v_0) + error("Intial data x_0 and v_0 must have equal length.") + end + + tc = tspan[1] + tf = tspan[end] + h = tf - tc + v = v_0 + x = x_0 + + tout = tc + xout = x_0.' + vout = v_0.' + + while tc < tf + vo = v + xo = x + a_tmp = a(tc, xo) + + # try a full step + v = vo + a_tmp*h/2 + x1 = xo + v*h + v1 = v + a(tc+h, x)*h/2 + + # try two half steps + h = h/2 + + v = vo + a_tmp*h/2 + x = xo + v*h + a_tmp = a(tc+h, x) + v = v + a_tmp*h/2 + + v = v + a_tmp*h/2 + x = x + v*h + v = v + a(tc+2*h, x)*h/2 + + # error estimate + err = (8/7)*norm(x1 - x) + hmax = 2*h*abs(atol/err)^(1/4) + + if hmax < h + # repeat step with smaller step size + h = 0.9*2*hmax + v = vo + x = xo + else + # accept last step + tc = tc + 2*h + h = 2*h + + tout = [tout; tc] + xout = [xout; x.'] + vout = [vout; v.'] + end + end + return tout, xout, vout +end + +# use adaptive method by default +verlet(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm)= verlet_hh2(a, tspan, x_0, v_0; atol=atol, norm=norm) From 53a7016ad1ebabef7740d2ba790034f80d3054cf Mon Sep 17 00:00:00 2001 From: acroy Date: Mon, 27 Jan 2014 21:36:17 +0100 Subject: [PATCH 2/4] Fixed indentation & typos. Added test. --- src/symplectic.jl | 150 +++++++++++++++++++++++----------------------- test/tests.jl | 5 ++ 2 files changed, 80 insertions(+), 75 deletions(-) diff --git a/src/symplectic.jl b/src/symplectic.jl index 9daab84d7..f03386012 100644 --- a/src/symplectic.jl +++ b/src/symplectic.jl @@ -14,29 +14,29 @@ ## The output, (t,x,v), consists of the solutions x(t), v(t) at times t=tspan. ## function verlet_fixed(a, tspan, x_0, v_0) - - if length(x_0) != length(v_0) - error("Intial data x_0 and v_0 must have equal length.") - end - - x = Array(eltype(x_0), length(tspan), length(x_0)) - v = Array(eltype(v_0), length(tspan), length(v_0)) - - x[1,:] = x_0' - v[1,:] = v_0' - - atmp = a(tspan[1], x_0).' - + + if length(x_0) != length(v_0) + error("Initial data x_0 and v_0 must have equal length.") + end + + x = Array(eltype(x_0), length(tspan), length(x_0)) + v = Array(eltype(v_0), length(tspan), length(v_0)) + + x[1,:] = x_0' + v[1,:] = v_0' + + atmp = a(tspan[1], x_0).' + for i = 1:(length(tspan)-1) - h = tspan[i+1] - tspan[i] - - v[i+1,:] = v[i,:] + atmp*h/2 - x[i+1,:] = x[i,:] + v[i+1,:]*h - atmp = a(tspan[i]+h, x[i+1,:].').' - v[i+1,:] = v[i+1,:] + atmp*h/2 - end - - return tspan, x, v + h = tspan[i+1] - tspan[i] + + v[i+1,:] = v[i,:] + atmp*h/2 + x[i+1,:] = x[i,:] + v[i+1,:]*h + atmp = a(tspan[i]+h, x[i+1,:].').' + v[i+1,:] = v[i+1,:] + atmp*h/2 + end + + return tspan, x, v end ## velocity verlet with adaptive step-size @@ -57,62 +57,62 @@ end ## intermediate integration times t including tspan[1] and tspan[end]. ## function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm) - if length(x_0) != length(v_0) - error("Intial data x_0 and v_0 must have equal length.") - end + if length(x_0) != length(v_0) + error("Initial data x_0 and v_0 must have equal length.") + end - tc = tspan[1] - tf = tspan[end] - h = tf - tc - v = v_0 - x = x_0 - - tout = tc - xout = x_0.' - vout = v_0.' + tc = tspan[1] + tf = tspan[end] + h = tf - tc + v = v_0 + x = x_0 + + tout = tc + xout = x_0.' + vout = v_0.' - while tc < tf - vo = v - xo = x - a_tmp = a(tc, xo) - - # try a full step - v = vo + a_tmp*h/2 - x1 = xo + v*h - v1 = v + a(tc+h, x)*h/2 - - # try two half steps - h = h/2 - - v = vo + a_tmp*h/2 - x = xo + v*h - a_tmp = a(tc+h, x) - v = v + a_tmp*h/2 + while tc < tf + vo = v + xo = x + a_tmp = a(tc, xo) + + # try a full step + v = vo + a_tmp*h/2 + x1 = xo + v*h + v1 = v + a(tc+h, x)*h/2 + + # try two half steps + h = h/2 + + v = vo + a_tmp*h/2 + x = xo + v*h + a_tmp = a(tc+h, x) + v = v + a_tmp*h/2 - v = v + a_tmp*h/2 - x = x + v*h - v = v + a(tc+2*h, x)*h/2 - - # error estimate - err = (8/7)*norm(x1 - x) - hmax = 2*h*abs(atol/err)^(1/4) - - if hmax < h - # repeat step with smaller step size - h = 0.9*2*hmax - v = vo - x = xo - else - # accept last step - tc = tc + 2*h - h = 2*h - - tout = [tout; tc] - xout = [xout; x.'] - vout = [vout; v.'] - end - end - return tout, xout, vout + v = v + a_tmp*h/2 + x = x + v*h + v = v + a(tc+2*h, x)*h/2 + + # error estimate + err = (8/7)*norm(x1 - x) + hmax = 2*h*abs(atol/err)^(1/4) + + if hmax < h + # repeat step with smaller step size + h = 0.9*2*hmax + v = vo + x = xo + else + # accept last step + tc = tc + 2*h + h = 2*h + + tout = [tout; tc] + xout = [xout; x.'] + vout = [vout; v.'] + end + end + return tout, xout, vout end # use adaptive method by default diff --git a/test/tests.jl b/test/tests.jl index 87df8384e..7f730e2ef 100644 --- a/test/tests.jl +++ b/test/tests.jl @@ -43,4 +43,9 @@ for solver in solvers @test maximum(abs(y-[cos(t)-2*sin(t) 2*cos(t)+sin(t)])) < tol end +println("using ODE.verlet") +t,x,v=ODE.verlet((t,x)-> -x, [0,2pi], 2., 1.) +@test maximum(abs([v - cos(t)+2*sin(t), x - 2*cos(t)-sin(t)])) < tol + + println("All looks OK") From 5ad0b0cd666165430c9d2db7fa1000adfc97aa7d Mon Sep 17 00:00:00 2001 From: acroy Date: Fri, 14 Feb 2014 14:08:22 +0100 Subject: [PATCH 3/4] Keyword points=:all|:specified added to adaptive solver. --- src/symplectic.jl | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/symplectic.jl b/src/symplectic.jl index f03386012..aaff63cd1 100644 --- a/src/symplectic.jl +++ b/src/symplectic.jl @@ -53,10 +53,14 @@ end ## specifies a list of times for which a solution is requested. The last ## two arguments are the initial conditions x_0 and v_0. ## +## Optional keywords are atol (error tolerance), norm (function to calculate the error) +## and points=:all|:specified (flag to indicate type of output). +## ## The output, (t,x,v), consists of the solutions x(t), v(t) at the -## intermediate integration times t including tspan[1] and tspan[end]. +## intermediate integration times t including tspan[1] and tspan[end] for points==:all +## or at tspan[1] and tspan[end] for points==:specified ## -function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm) +function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Symbol=:all) if length(x_0) != length(v_0) error("Initial data x_0 and v_0 must have equal length.") end @@ -107,13 +111,21 @@ function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm) tc = tc + 2*h h = 2*h - tout = [tout; tc] - xout = [xout; x.'] - vout = [vout; v.'] + if points == :all + tout = [tout; tc] + xout = [xout; x.'] + vout = [vout; v.'] + end end end + + if points == :specified + tout = [tout; tc] + xout = [xout; x.'] + vout = [vout; v.'] + end return tout, xout, vout end # use adaptive method by default -verlet(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm)= verlet_hh2(a, tspan, x_0, v_0; atol=atol, norm=norm) +verlet(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Symbol=:all)= verlet_hh2(a, tspan, x_0, v_0; atol=atol, norm=norm, points=points) From 8525b8b28c2d3abf504c59576a97ba69c1af3f04 Mon Sep 17 00:00:00 2001 From: acroy Date: Sun, 16 Feb 2014 21:14:29 +0100 Subject: [PATCH 4/4] Support of user specified intermediate points. --- src/symplectic.jl | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/symplectic.jl b/src/symplectic.jl index aaff63cd1..4f806966d 100644 --- a/src/symplectic.jl +++ b/src/symplectic.jl @@ -50,23 +50,22 @@ end ## ## The function takes as a first argument a function, a(t,x), ## of time and position (the acceleration). The second argument -## specifies a list of times for which a solution is requested. The last +## specifies an intial time t_0 and a final time tf for which a solution is requested. The last ## two arguments are the initial conditions x_0 and v_0. ## ## Optional keywords are atol (error tolerance), norm (function to calculate the error) ## and points=:all|:specified (flag to indicate type of output). ## ## The output, (t,x,v), consists of the solutions x(t), v(t) at the -## intermediate integration times t including tspan[1] and tspan[end] for points==:all -## or at tspan[1] and tspan[end] for points==:specified +## intermediate integration times t including ti and tf for points==:all +## or at ti and tf only for points==:specified ## -function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Symbol=:all) +function verlet_hh2(a, t_0, tf, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Symbol=:all) if length(x_0) != length(v_0) error("Initial data x_0 and v_0 must have equal length.") end - tc = tspan[1] - tf = tspan[end] + tc = t_0 h = tf - tc v = v_0 x = x_0 @@ -128,4 +127,25 @@ function verlet_hh2(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Sym end # use adaptive method by default -verlet(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Symbol=:all)= verlet_hh2(a, tspan, x_0, v_0; atol=atol, norm=norm, points=points) +function verlet(a, tspan, x_0, v_0; atol = 1e-5, norm=Base.norm, points::Symbol=:all) + if length(tspan) < 2 + error("Argument tspan must have at least two elements.") + end + + tout, xout, vout = verlet_hh2(a, tspan[1], tspan[2], x_0, v_0, atol=atol, norm=norm, points=points) + for i = 2:(length(tspan)-1) + t, x, v = verlet_hh2(a, tout[end], tspan[i+1], xout[end,:].', vout[end,:].', atol=atol, norm=norm, points=points) + + if points == :all + tout = [tout; t[2:end]] + xout = [xout; x[2:end,:]] + vout = [vout; v[2:end,:]] + elseif points == :specified + tout = [tout; t[end]] + xout = [xout; x[end,:]] + vout = [vout; v[end,:]] + end + end + + return tout, xout, vout +end