Some crude benchmarking for ODE.jl

perf_result.txt

== problem A3 in DETEST ==

testing method ode23
elapsed time: 0.069644091 seconds (12007840 bytes allocated)
*  number of steps: 307
*  minimal step: 0.017235477520255074
*  maximal step: 0.13251841618852822
*  maximal deviation from known solution: 0.00024150339669093412

testing method ode23_bs
elapsed time: 0.049506498 seconds (9628320 bytes allocated)
*  number of steps: 309
*  minimal step: 0.030318943004498777
*  maximal step: 0.13334865624485204
*  maximal deviation from known solution: 0.0002443491465893288

testing method ode45_dp
elapsed time: 0.010797216 seconds (1163120 bytes allocated)
*  number of steps: 52
*  minimal step: 0.14907708023597976
*  maximal step: 0.5277942154717685
*  maximal deviation from known solution: 1.4358013608717357e-5

testing method ode45_fb
elapsed time: 0.009037712 seconds (1255360 bytes allocated)
*  number of steps: 57
*  minimal step: 0.2
*  maximal step: 0.6019190646768298
*  maximal deviation from known solution: 0.00029615990161868666

testing method ode45_ck
elapsed time: 0.007782929 seconds (1071120 bytes allocated)
*  number of steps: 50
*  minimal step: 0.2
*  maximal step: 0.6412708946258334
*  maximal deviation from known solution: 0.00015482199537864005

testing method ode78_fb
elapsed time: 0.004410885 seconds (642960 bytes allocated)
*  number of steps: 26
*  minimal step: 0.06587922675161195
*  maximal step: 1.0422081956865448
*  maximal deviation from known solution: 1.1710040716206294e-5


== problem B5 in DETEST ==

testing method ode23
elapsed time: 0.073563663 seconds (15847520 bytes allocated)
*  number of steps: 238
*  minimal step: 1.7235477520255075e-5
*  maximal step: 0.13058852624073225


testing method ode23_bs
elapsed time: 0.043170629 seconds (10955760 bytes allocated)
*  number of steps: 172
*  minimal step: 0.0007839590369780325
*  maximal step: 0.18195874964159842


testing method ode45_dp
elapsed time: 0.011750347 seconds (3670320 bytes allocated)
*  number of steps: 34
*  minimal step: 0.12
*  maximal step: 0.5079142619969694


testing method ode45_fb
elapsed time: 0.011377746 seconds (3475760 bytes allocated)
*  number of steps: 36
*  minimal step: 0.12
*  maximal step: 0.5526493932284717


testing method ode45_ck
elapsed time: 0.010097592 seconds (2588080 bytes allocated)
*  number of steps: 28
*  minimal step: 0.12
*  maximal step: 0.548455619189153


testing method ode78_fb
elapsed time: 0.014098016 seconds (5003920 bytes allocated)
*  number of steps: 16
*  minimal step: 0.12
*  maximal step: 1.1468526118328999



== problem Gravity in DETEST ==

testing method ode23
elapsed time: 1.476745924 seconds (699598960 bytes allocated)
*  number of steps: 3458
*  minimal step: 1.7969148347567307
*  maximal step: 39.06078198352043


testing method ode23_bs
elapsed time: 0.386518477 seconds (202950000 bytes allocated)
*  number of steps: 1609
*  minimal step: 24.995537553550093
*  maximal step: 58.107348672667285


testing method ode45_dp
elapsed time: 0.095624438 seconds (41533040 bytes allocated)
*  number of steps: 304
*  minimal step: 232.63351500181307
*  maximal step: 296.21785556010127


testing method ode45_fb
elapsed time: 0.090328791 seconds (38245520 bytes allocated)
*  number of steps: 334
*  minimal step: 199.54950236153672
*  maximal step: 269.97606116933457


testing method ode45_ck
elapsed time: 0.067431822 seconds (26644400 bytes allocated)
*  number of steps: 240
*  minimal step: 283.06228400213877
*  maximal step: 378.17764195700147


testing method ode78_fb
elapsed time: 0.107533446 seconds (45735280 bytes allocated)
*  number of steps: 141
*  minimal step: 196.05248931393726
*  maximal step: 635.4255672814579

perf_test.jl

# comparison of adaptive methods

function runproblem(problem, solver; nruns::Int=10)

  (pname, hassol, F, y0, t0, tf, soly) = problem
  # println("\n== problem $(pname) in DETEST ==\n")

  println("testing method $(string(solver))")

  t,y = solver(F, y0, [t0, tf]);
  gc_disable()
  tau = @time begin
      for i=1:nruns
          t,y = solver(F, y0, [t0, tf]);
      end
  end
  gc_enable()

  # println("*  elapsed time ($(nruns) runs): $(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


# 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

function gravity(t::Float64, y::Vector{Float64})
    const mu = 398600.4415 # [km^3/s^2] Earth's gravity
    f = similar(y)
    r = sqrt(y[1]*y[1]+y[2]*y[2]+y[3]*y[3])
    r3 = r*r*r
    f[1] = y[4]
    f[2] = y[5]
    f[3] = y[6]
    f[4] = -mu*y[1]/r3
    f[5] = -mu*y[2]/r3
    f[6] = -mu*y[3]/r3

    return f
end

const s0 = Float64[-1814.0, -3708.0, 5153.0, 6.512, -4.229, -0.744]
const tp = 5402.582703094263*16

# selected problems from DETEST
# problem consists of {name, has_solution flag, rhs, x0, t0, tend, solution}
problems = ({"A3", true, (t,y)->y*cos(t), 1., 0., 20., (t)->exp(sin(t))},
            {"B5", false, rigid, [0., 1., 1.], 0., 12., nothing},
            {"Gravity", false, gravity, s0, 0., tp, nothing},
            )

# adaptive solvers
solvers = [
    ODE.ode23,
    ODE.ode23_bs,

    ODE.ode45_dp,
    ODE.ode45_fb,
    ODE.ode45_ck,
    ODE.ode78_fb
    ]

for problem in problems
  (pname, hassol, F, y0, t0, tf, soly) = problem
  println("\n== problem $(pname) in DETEST ==\n")

  for solver in solvers
    runproblem(problem, solver)
  end
end