Add more tests

Author: tansongchen

bench/adaptive_order.jl

@@ -18,6 +18,18 @@ end
 prob = ODEProblem{true, SciMLBase.FullSpecialize}(f, [100.0, 100.0], (0.0, 10.0), nothing)
 sol = solve(prob, Vern7(), abstol = 1 / 10^14, reltol = 1 / 10^14);
 
+function f2(du, u, p, t)
+    t0 = 5.0
+    δ = 1
+    κ = 10.0
+    y = u[1]
+    ϕ1 = κ / δ / cosh((t - t0) / δ)^2
+    du[1] = ϕ1 * (y - y^2)
+end
+
+prob2 = ODEProblem{true, SciMLBase.FullSpecialize}(f2, [0.5], (0.0, 10.0), nothing)
+sol2 = solve(prob2, Vern7(), abstol = 1 / 10^14, reltol = 1 / 10^14);
+
 # Profile
 function profile1()
     cache1 = init(prob, ExplicitTaylor(order = Val(8)))
@@ -35,6 +47,6 @@ for order in 6:2:12
     push!(names, "Taylor $(order)")
     push!(setups, Dict(:alg => ExplicitTaylor(order = Val(order + 1))))
 end
-wp = WorkPrecisionSet([prob], abstols, reltols, setups; names = names, appxsol = [sol],
+wp = WorkPrecisionSet([prob2], abstols, reltols, setups; names = names, appxsol = [sol2],
     save_everystep = false, numruns = 100, maxiters = 10000)
 plot(wp)

bench/non_stiff.jl

@@ -1,20 +1,71 @@
 using OrdinaryDiffEq, ParameterizedFunctions, DiffEqDevTools, StaticArrays,
       OrdinaryDiffEqTaylorSeries
-using Plots, BenchmarkTools
+using Plots, BenchmarkTools, Random, LinearAlgebra
 gr();
 
-abstols = 1.0 ./ 10.0 .^ (6:15)
-reltols = 1.0 ./ 10.0 .^ (3:12);
-
 function lotka(du, u, p, t)
     x = u[1]
     y = u[2]
     du[1] = p[1] * x - p[2] * x * y
     du[2] = -p[3] * y + p[4] * x * y
 end
-prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+prob_lotka = ODEProblem{true, SciMLBase.FullSpecialize}(
     lotka, [1.0, 1.0], (0.0, 10.0), [1.5, 1.0, 3.0, 1.0])
-sol = solve(prob, Vern7(), abstol = 1 / 10^14, reltol = 1 / 10^14)
+sol_lotka = solve(prob_lotka, Vern7(), abstol = 1e-14, reltol = 1e-14);
+
+function fitzhugh(du, u, p, t)
+    v = u[1]
+    w = u[2]
+    a = p[1]
+    b = p[2]
+    τinv = p[3]
+    l = p[4]
+    du[1] = v - v^3 / 3 - w + l
+    du[2] = τinv * (v + a - b * w)
+end
+prob_fitzhugh = ODEProblem{true, SciMLBase.FullSpecialize}(
+    fitzhugh, [1.0; 1.0], (0.0, 10.0), [0.7, 0.8, 1 / 12.5, 0.5])
+sol_fitzhugh = solve(prob_fitzhugh, Vern7(), abstol = 1e-14, reltol = 1e-14);
+
+function rigid(du, u, p, t)
+    I₁ = p[1]
+    I₂ = p[2]
+    I₃ = p[3]
+    du[1] = I₁ * u[2] * u[3]
+    du[2] = I₂ * u[1] * u[3]
+    du[3] = I₃ * u[1] * u[2] + 0.25 * sin(t)^2
+end
+prob_rigid = ODEProblem{true, SciMLBase.FullSpecialize}(
+    rigid, [1.0; 0.0; 0.9], (0.0, 10.0), [-2.0, 1.25, -0.5])
+sol_rigid = solve(prob_rigid, Vern7(), abstol = 1e-14, reltol = 1e-14);
+
+function make_problem_random_dense(N; seed = 1)
+    Random.seed!(seed)
+    λ = -rand(N) .- 1   # eigenvalues in [-2,-1]
+    D = Diagonal(λ)
+    Q, _ = qr(randn(N, N))            # random orthogonal matrix
+    A = Q * D * Q'                    # similar transformation
+    f(du, u, p, t) = du .= A * u
+    u0 = rand(N)
+    ODEProblem{true, SciMLBase.FullSpecialize}(f, u0, (0.0, 10.0))
+end
+
+function make_problem_random_sparse(N; seed = 1, density = 0.1)
+    Random.seed!(seed)
+    λ = -rand(N) .- 1   # eigenvalues in [-2,-1]
+    D = Diagonal(λ)
+    Q, _ = qr(randn(N, N))            # random orthogonal matrix
+    A_dense = Q * D * Q'               # similar transformation
+    A_sparse = sparse(A_dense)
+    A_sparse .= sprand(N, N, density) .* A_sparse  # make sparse but keep structure
+    f(du, u, p, t) = du .= A_sparse * u
+    u0 = rand(N)
+    ODEProblem{true, SciMLBase.FullSpecialize}(f, u0, (0.0, 10.0))
+end
+
+# prob_dense = make_problem_random_dense(64)
+prob_dense = make_problem_poisson(16)
+sol_dense = solve(prob_dense, Vern7(), abstol = 1e-14, reltol = 1e-14);
 
 # Make work-precision plot
 setups = [Dict(:alg => DP5())
@@ -26,6 +77,17 @@ for order in 6:2:12
     push!(names, "Taylor $(order)")
     push!(setups, Dict(:alg => ExplicitTaylor(order = Val(order + 1))))
 end
-wp = WorkPrecisionSet([prob], abstols, reltols, setups; names = names, appxsol = [sol],
-    save_everystep = false, numruns = 100, maxiters = 10000)
-plot(wp)
+
+function make_plot(
+        prob, sol; numruns = 100, maxiters = 1000, absrange = (-13:-6), relrange = (-10:-3))
+    abstols = 10.0 .^ absrange
+    reltols = 10.0 .^ relrange
+    wp = WorkPrecisionSet([prob], abstols, reltols, setups; names = names, appxsol = [sol],
+        save_everystep = false, numruns = numruns, maxiters = maxiters)
+    p = plot(wp)
+    return p
+end
+p_lotka = make_plot(prob_lotka, sol_lotka)
+p_fitzhugh = make_plot(prob_fitzhugh, sol_fitzhugh)
+p_rigid = make_plot(prob_rigid, sol_rigid)
+p_dense = make_plot(prob_dense, sol_dense)

bench/threebody.jl

@@ -31,12 +31,30 @@ function py!(q0, J0)
 end
 q0 = [-0.8, 0.0, 0.0, 0.0]
 py!(q0, J0)
-tspan = (0.0, 2000.0)
+tspan = (0.0, 100.0)
 p = SA[μ]
 prob = ODEProblem{true, SciMLBase.FullSpecialize}(pcr3bp!, q0, tspan, p)
 
-ref = solve(prob, Vern9(), abstol = 1e-10, reltol = 1e-10)
+ref = solve(prob, TaylorMethod(32), abstol = 1e-20)
 
 @benchmark solve($prob, $(TaylorMethod(25)), abstol = 1e-15)
-@benchmark solve($prob, $(ExplicitTaylor(order = Val(25))), abstol = 1e-15, reltol = 1e-15)
-@benchmark solve($prob, $(Vern9()), abstol = 1e-15, reltol = 1e-15)
+@benchmark solve($prob, $(ExplicitTaylor(order = Val(11))), abstol = 1e-10, reltol = 1e-10)
+@benchmark solve($prob, $(Vern9()), abstol = 1e-10, reltol = 1e-10)
+
+setups = [Dict(:alg => Vern7())
+          Dict(:alg => Vern9())
+          Dict(:alg => ExplicitTaylor(order = Val(13)))
+          Dict(:alg => ExplicitTaylor(order = Val(17)))
+          Dict(:alg => ExplicitTaylor(order = Val(21)))
+          Dict(:alg => ExplicitTaylor(order = Val(25)))
+          Dict(:alg => TaylorMethod(12))
+          Dict(:alg => TaylorMethod(16))
+          Dict(:alg => TaylorMethod(20))
+          Dict(:alg => TaylorMethod(24))]
+abstols = 10.0 .^ (-19:-15)
+reltols = 10.0 .^ (-15:-11)
+names = ["Vern7", "Vern9", "Taylor12", "Taylor16", "Taylor20", "Taylor24", "TI12", "TI16",
+    "TI20", "TI24"]
+wp = WorkPrecisionSet([prob], abstols, reltols, setups; names = names, appxsol = [ref],
+    save_everystep = false, numruns = 100)
+plot(wp)