From af2c6a6628b5da57d5db3dee404d57c038589a9f Mon Sep 17 00:00:00 2001 From: acroy Date: Sun, 2 Mar 2014 20:04:27 +0100 Subject: [PATCH 1/7] Generalize oderkf to support orders != 5. Add ode23_bs based on oderkf. --- src/ODE.jl | 31 ++++++++++++++++++++++++++----- 1 file changed, 26 insertions(+), 5 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 116a0fcbe..b9c73a86c 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -183,9 +183,9 @@ end # ode23 # modified: 17 January 2001 function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) # see p.91 in the Ascher & Petzold reference for more infomation. - pow = 1/5 + pow = 1/p # use the higher order to estimate the next step size - c = sum(a, 2) + c = sum(a, 2) # consistency condition # Initialization t = tspan[1] @@ -255,13 +255,28 @@ function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) return tout, xout end +# Bogacki–Shampine coefficients +const bs_coefficients = (3, + [ 0 0 0 0 + 1/2 0 0 0 + 0 3/4 0 0 + 2/9 1/3 4/9 0], + # 2nd order b-coefficients + [7/24 1/4 1/3 1/8], + # 3rd order b-coefficients + [2/9 1/3 4/9 0], + ) +ode23_bs(F, tspan, x0) = oderkf(F, tspan, x0, bs_coefficients...) + + # Both the Dormand-Prince and Fehlberg 4(5) coefficients are from a tableau in # U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations # and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics # (SIAM), Philadelphia, 1998 # # Dormand-Prince coefficients -const dp_coefficients = ([ 0 0 0 0 0 0 +const dp_coefficients = (5, + [ 0 0 0 0 0 0 1/5 0 0 0 0 0 3/40 9/40 0 0 0 0 44/45 -56/15 32/9 0 0 0 @@ -276,7 +291,8 @@ const dp_coefficients = ([ 0 0 0 0 0 ode45_dp(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, dp_coefficients...; kwargs...) # Fehlberg coefficients -const fb_coefficients = ([ 0 0 0 0 0 +const fb_coefficients = (5, + [ 0 0 0 0 0 1/4 0 0 0 0 3/32 9/32 0 0 0 1932/2197 -7200/2197 7296/2197 0 0 @@ -291,7 +307,8 @@ ode45_fb(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, fb_coefficients...; kwa # Cash-Karp coefficients # Numerical Recipes in Fortran 77 -const ck_coefficients = ([ 0 0 0 0 0 +const ck_coefficients = (5, + [ 0 0 0 0 0 1/5 0 0 0 0 3/40 9/40 0 0 0 3/10 -9/10 6/5 0 0 @@ -307,6 +324,10 @@ ode45_ck(F, x0, tspan; kwargs...) = oderkf(F, x0, tspan, ck_coefficients...; kwa # Use Dormand Prince version of ode45 by default const ode45 = ode45_dp +# some higher-order embedded methods can be found in: +# P.J. Prince and J.R.Dormand: High order embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 7(1), 1981. + + #ODE4 Solve non-stiff differential equations, fourth order # fixed-step Runge-Kutta method. # From be5ca6e960b05b7ffc5cc61a5367ad7d217766d3 Mon Sep 17 00:00:00 2001 From: acroy Date: Sun, 2 Mar 2014 20:53:46 +0100 Subject: [PATCH 2/7] Add performance comparison for adaptive solvers. --- test/perf_test.jl | 29 +++++++++++++++++++++++++++++ test/runtests.jl | 4 +++- 2 files changed, 32 insertions(+), 1 deletion(-) create mode 100644 test/perf_test.jl diff --git a/test/perf_test.jl b/test/perf_test.jl new file mode 100644 index 000000000..cf66ecf90 --- /dev/null +++ b/test/perf_test.jl @@ -0,0 +1,29 @@ +# comparison of adaptive methods + +# adaptive solvers +solvers = [ + ODE.ode23, + ODE.ode23_bs, + + ODE.ode45_dp, + ODE.ode45_fb, + ODE.ode45_ck + ] + +for solver in solvers + println("testing method $(string(solver))") + + # problem A3 in DETEST + t,y = solver((t,y)->y*cos(t), [0.,20.], [1.]); + @time begin + for i=1:10 + t,y = solver((t,y)->y*cos(t), [0.,20.], [1.]); + end + end + + println("* number of steps: $(length(t))") + println("* minimal step: $(minimum(diff(t)))") + println("* maximal step: $(maximum(diff(t)))") + + println("* maximal deviation from known solution: $(maximum(abs(y[:]-exp(sin(t)))))\n") +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 6ad84f4a4..e23461fb9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,11 +5,13 @@ tol = 1e-2 solvers = [ ODE.ode23, + ODE.ode23_bs, + ODE.ode4, ODE.ode45_dp, ODE.ode45_fb, ODE.ode45_ck, - + ODE.ode4ms, ODE.ode4s_s, ODE.ode4s_kr] From dcc26615a713ad07157ecf156897a41584abdcca Mon Sep 17 00:00:00 2001 From: acroy Date: Tue, 4 Mar 2014 20:55:03 +0100 Subject: [PATCH 3/7] Some cosmetics & new test added (B5). --- src/ODE.jl | 6 +++--- test/perf_test.jl | 39 ++++++++++++++++++++++++++++++++++++--- 2 files changed, 39 insertions(+), 6 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index b9c73a86c..aa1f61a50 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -181,7 +181,7 @@ end # ode23 # CompereM@asme.org # created : 06 October 1999 # modified: 17 January 2001 -function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) +function oderkf(F, x0, tspan, p, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) # see p.91 in the Ascher & Petzold reference for more infomation. pow = 1/p # use the higher order to estimate the next step size @@ -217,7 +217,7 @@ function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) x5 = x + h.*(b5*k)[1] # estimate the local truncation error - gamma1 = x5 - x4 + gamma1 = x4 - x5 # Estimate the error and the acceptable error delta = norm(gamma1, Inf) # actual error @@ -233,7 +233,7 @@ function oderkf(F, x0, tspan, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) # Compute the slopes by computing the k[:,j+1]'th column based on the previous k[:,1:j] columns # notes: k needs to end up as an Nxs, a is 7x6, which is s by (s-1), # s is the number of intermediate RK stages on [t (t+h)] (Dormand-Prince has s=7 stages) - if c[end] == 1 && t != tspan[1] + if c[end] == 1 # Assign the last stage for x(k) as the first stage for computing x[k+1]. # This is part of the Dormand-Prince pair caveat. # k[:,7] has already been computed, so use it instead of recomputing it diff --git a/test/perf_test.jl b/test/perf_test.jl index cf66ecf90..42b97a2e4 100644 --- a/test/perf_test.jl +++ b/test/perf_test.jl @@ -1,5 +1,6 @@ # comparison of adaptive methods + # adaptive solvers solvers = [ ODE.ode23, @@ -10,14 +11,15 @@ solvers = [ ODE.ode45_ck ] +println("\n== problem A3 in DETEST ==\n") for solver in solvers println("testing method $(string(solver))") # problem A3 in DETEST - t,y = solver((t,y)->y*cos(t), [0.,20.], [1.]); + t,y = solver((t,y)->y*cos(t), 1., [0.,20.]); @time begin for i=1:10 - t,y = solver((t,y)->y*cos(t), [0.,20.], [1.]); + t,y = solver((t,y)->y*cos(t), 1., [0.,20.]); end end @@ -26,4 +28,35 @@ for solver in solvers println("* maximal step: $(maximum(diff(t)))") println("* maximal deviation from known solution: $(maximum(abs(y[:]-exp(sin(t)))))\n") -end \ No newline at end of file +end + +# problem B5 in DETEST +# motion of a rigid body without external forces +# (see also Example 1 from http://www.mathworks.se/help/matlab/ref/ode45.html) +function rigid(t,y) + dy = copy(y) + dy[1] = y[2] * y[3] + dy[2] = -y[1] * y[3] + dy[3] = -0.51 * y[1] * y[2] + + return dy +end + +println("\n== problem B5 in DETEST ==\n") +for solver in solvers + println("testing method $(string(solver))") + + # problem B5 in DETEST + t,y = solver(rigid, [0., 1., 1.], [0.,12.]); + @time begin + for i=1:10 + t,y = solver(rigid, [0., 1., 1.], [0.,12.]); + end + end + + println("* number of steps: $(length(t))") + println("* minimal step: $(minimum(diff(t)))") + println("* maximal step: $(maximum(diff(t)))\n") + + # println("* maximal deviation from known solution: $(maximum(abs(y[:]-exp(sin(t)))))\n") +end From 4e1c417a15da79b6b660ff65502873e6c82f533d Mon Sep 17 00:00:00 2001 From: acroy Date: Wed, 12 Mar 2014 22:18:39 +0100 Subject: [PATCH 4/7] Devectorization in oderkf, oderosenbrock and ode_ms (fixes problems related to JuliaLang/julia#6118). --- src/ODE.jl | 48 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index aa1f61a50..356768654 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -181,7 +181,7 @@ end # ode23 # CompereM@asme.org # created : 06 October 1999 # modified: 17 January 2001 -function oderkf(F, x0, tspan, p, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) +function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8) # see p.91 in the Ascher & Petzold reference for more infomation. pow = 1/p # use the higher order to estimate the next step size @@ -206,18 +206,24 @@ function oderkf(F, x0, tspan, p, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) h = tfinal - t end + #(p-1)th and pth order estimates + xs = x + h*bs[1]*k[1] + xp = x + h*bp[1]*k[1] for j = 2:length(c) - k[j] = F(t + h.*c[j], x + h.*(a[j,1:j-1]*k[1:j-1])[1]) - end - - # compute the 4th order estimate - x4 = x + h.*(b4*k)[1] + dx = a[j,1]*k[1] + for i = 2:j-1 + dx += a[j,i]*k[i] + end + k[j] = F(t + h*c[j], x + h*dx) - # compute the 5th order estimate - x5 = x + h.*(b5*k)[1] + # compute the (p-1)th order estimate + xs = xs + h*bs[j]*k[j] + # compute the pth order estimate + xp = xp + h*bp[j]*k[j] + end # estimate the local truncation error - gamma1 = x4 - x5 + gamma1 = xs - xp # Estimate the error and the acceptable error delta = norm(gamma1, Inf) # actual error @@ -226,7 +232,7 @@ function oderkf(F, x0, tspan, p, a, b4, b5; reltol = 1.0e-5, abstol = 1.0e-8) # Update the solution only if the error is acceptable if delta <= tau t = t + h - x = x5 # <-- using the higher order estimate is called 'local extrapolation' + x = xp # <-- using the higher order estimate is called 'local extrapolation' tout = [tout; t] push!(xout, x) @@ -404,15 +410,23 @@ function oderosenbrock(F, x0, tspan, gamma, a, b, c; jacobian=nothing) if size(dFdx,1) == 1 jac = 1/gamma/hs - dFdx[1] else - jac = eye(dFdx)./gamma./hs - dFdx + jac = eye(dFdx)/gamma/hs - dFdx end g = Array(typeof(x0), size(a,1)) - g[1] = (jac \ F(ts + b[1].*hs, xs)) + g[1] = (jac \ F(ts + b[1]*hs, xs)) + x[solstep+1] = x[solstep] + b[1]*g[1] + for i = 2:size(a,1) - g[i] = (jac \ (F(ts + b[i].*hs, xs + (a[i,1:i-1]*g[1:i-1])[1]) + (c[i,1:i-1]*g[1:i-1])[1]./hs)) + dx = zero(x0) + dF = zero(x0/hs) + for j = 1:i-1 + dx += a[i,j]*g[j] + dF += c[i,j]*g[j] + end + g[i] = (jac \ (F(ts + b[i]*hs, xs + dx) + dF/hs)) + x[solstep+1] += b[i]*g[i] end - x[solstep+1] = x[solstep] + (b*g)[1] solstep += 1 end return [tspan], x @@ -478,7 +492,11 @@ function ode_ms(F, x0, tspan, order::Integer) # Need to run the first several steps at reduced order steporder = min(i, order) xdot[i] = F(tspan[i], x[i]) - x[i+1] = x[i] + (b[steporder,1:steporder]*xdot[i-(steporder-1):i])[1].*h[i] + + x[i+1] = x[i] + for j = 1:steporder + x[i+1] += h[i]*b[steporder, j]*xdot[i-(steporder-1) + (j-1)] + end end return [tspan], x end From f1f00410986206dd699ff7ca664f449bef4356dc Mon Sep 17 00:00:00 2001 From: acroy Date: Thu, 13 Mar 2014 13:13:48 +0100 Subject: [PATCH 5/7] Cleaned up performance test. --- test/perf_test.jl | 79 +++++++++++++++++++++-------------------------- 1 file changed, 36 insertions(+), 43 deletions(-) diff --git a/test/perf_test.jl b/test/perf_test.jl index 42b97a2e4..7db8a1d53 100644 --- a/test/perf_test.jl +++ b/test/perf_test.jl @@ -1,62 +1,55 @@ # comparison of adaptive methods +# problem B5 in DETEST +# motion of a rigid body without external forces +# (see also Example 1 from http://www.mathworks.se/help/matlab/ref/ode45.html) +function rigid(t,y) + dy = similar(y) + dy[1] = y[2] * y[3] + dy[2] = -y[1] * y[3] + dy[3] = -0.51 * y[1] * y[2] + + return dy +end + +# selected problems from DETEST +problems = (["A3", true, (t,y)->y*cos(t), 1., 0., 20., (t)->exp(sin(t))], + {"B5", false, rigid, [0., 1., 1.], 0., 12., nothing},) # adaptive solvers solvers = [ ODE.ode23, ODE.ode23_bs, - + ODE.ode45_dp, ODE.ode45_fb, ODE.ode45_ck ] -println("\n== problem A3 in DETEST ==\n") -for solver in solvers - println("testing method $(string(solver))") - - # problem A3 in DETEST - t,y = solver((t,y)->y*cos(t), 1., [0.,20.]); - @time begin - for i=1:10 - t,y = solver((t,y)->y*cos(t), 1., [0.,20.]); - end - end +for problem in problems - println("* number of steps: $(length(t))") - println("* minimal step: $(minimum(diff(t)))") - println("* maximal step: $(maximum(diff(t)))") - - println("* maximal deviation from known solution: $(maximum(abs(y[:]-exp(sin(t)))))\n") -end + (pname, hassol, F, y0, t0, tf, soly) = problem + println("\n== problem $(pname) in DETEST ==\n") -# problem B5 in DETEST -# motion of a rigid body without external forces -# (see also Example 1 from http://www.mathworks.se/help/matlab/ref/ode45.html) -function rigid(t,y) - dy = copy(y) - dy[1] = y[2] * y[3] - dy[2] = -y[1] * y[3] - dy[3] = -0.51 * y[1] * y[2] + for solver in solvers + println("testing method $(string(solver))") - return dy -end + t,y = solver(F, y0, [t0, tf]); + tau = @elapsed begin + for i=1:10 + t,y = solver(F, y0, [t0, tf]); + end + end -println("\n== problem B5 in DETEST ==\n") -for solver in solvers - println("testing method $(string(solver))") - - # problem B5 in DETEST - t,y = solver(rigid, [0., 1., 1.], [0.,12.]); - @time begin - for i=1:10 - t,y = solver(rigid, [0., 1., 1.], [0.,12.]); + println("* elapsed time: $(tau)") + println("* number of steps: $(length(t))") + println("* minimal step: $(minimum(diff(t)))") + println("* maximal step: $(maximum(diff(t)))") + + if hassol + println("* maximal deviation from known solution: $(maximum(abs(y[:]-soly(t))))\n") + else + println("\n") end end - - println("* number of steps: $(length(t))") - println("* minimal step: $(minimum(diff(t)))") - println("* maximal step: $(maximum(diff(t)))\n") - - # println("* maximal deviation from known solution: $(maximum(abs(y[:]-exp(sin(t)))))\n") end From c7dc86f8d73e5eff9d6035751da3cc6b34af852d Mon Sep 17 00:00:00 2001 From: acroy Date: Fri, 14 Mar 2014 14:00:07 +0100 Subject: [PATCH 6/7] Added backward propagation in oderkf (part of #2). --- src/ODE.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 356768654..b02362e91 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -190,9 +190,10 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8) # Initialization t = tspan[1] tfinal = tspan[end] - hmax = (tfinal - t)/2.5 - hmin = (tfinal - t)/1e9 - h = (tfinal - t)/100 # initial guess at a step size + tdir = sign(tfinal - t) + hmax = abs(tfinal - t)/2.5 + hmin = abs(tfinal - t)/1e9 + h = tdir*abs(tfinal - t)/100 # initial guess at a step size x = x0 tout = t # first output time xout = Array(typeof(x0), 1) @@ -201,8 +202,8 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8) k = Array(typeof(x0), length(c)) k[1] = F(t,x) # first stage - while t < tfinal && h >= hmin - if t + h > tfinal + while abs(t) != abs(tfinal) && abs(h) >= hmin + if abs(h) > abs(tfinal-t) h = tfinal - t end @@ -254,8 +255,8 @@ function oderkf(F, x0, tspan, p, a, bs, bp; reltol = 1.0e-5, abstol = 1.0e-8) h = min(hmax, 0.8*h*(tau/delta)^pow) end # while (t < tfinal) & (h >= hmin) - if t < tfinal - println("Step size grew too small. t=", t, ", h=", h, ", x=", x) + if abs(t) < abs(tfinal) + println("Step size grew too small. t=", t, ", h=", abs(h), ", x=", x) end return tout, xout From 17a11e2878cf098b9dc73321bda85ff3ab3b6ed7 Mon Sep 17 00:00:00 2001 From: acroy Date: Mon, 5 May 2014 21:25:24 +0200 Subject: [PATCH 7/7] Remove performance test for now. --- test/perf_test.jl | 55 ----------------------------------------------- 1 file changed, 55 deletions(-) delete mode 100644 test/perf_test.jl diff --git a/test/perf_test.jl b/test/perf_test.jl deleted file mode 100644 index 7db8a1d53..000000000 --- a/test/perf_test.jl +++ /dev/null @@ -1,55 +0,0 @@ -# comparison of adaptive methods - -# problem B5 in DETEST -# motion of a rigid body without external forces -# (see also Example 1 from http://www.mathworks.se/help/matlab/ref/ode45.html) -function rigid(t,y) - dy = similar(y) - dy[1] = y[2] * y[3] - dy[2] = -y[1] * y[3] - dy[3] = -0.51 * y[1] * y[2] - - return dy -end - -# selected problems from DETEST -problems = (["A3", true, (t,y)->y*cos(t), 1., 0., 20., (t)->exp(sin(t))], - {"B5", false, rigid, [0., 1., 1.], 0., 12., nothing},) - -# adaptive solvers -solvers = [ - ODE.ode23, - ODE.ode23_bs, - - ODE.ode45_dp, - ODE.ode45_fb, - ODE.ode45_ck - ] - -for problem in problems - - (pname, hassol, F, y0, t0, tf, soly) = problem - println("\n== problem $(pname) in DETEST ==\n") - - for solver in solvers - println("testing method $(string(solver))") - - t,y = solver(F, y0, [t0, tf]); - tau = @elapsed begin - for i=1:10 - t,y = solver(F, y0, [t0, tf]); - end - end - - println("* elapsed time: $(tau)") - println("* number of steps: $(length(t))") - println("* minimal step: $(minimum(diff(t)))") - println("* maximal step: $(maximum(diff(t)))") - - if hassol - println("* maximal deviation from known solution: $(maximum(abs(y[:]-soly(t))))\n") - else - println("\n") - end - end -end