Merge pull request #66 from SciML/myb/rootedtrees

Add 4-th order Ralston and RootedTrees.jl integration

Project.toml

@@ -1,7 +1,7 @@
 name = "DiffEqDevTools"
 uuid = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "2.22.0"
+version = "2.23.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -12,6 +12,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+RootedTrees = "47965b36-3f3e-11e9-0dcf-4570dfd42a8c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
@@ -22,6 +23,7 @@ DiffEqProblemLibrary = "4.5"
 NLsolve = "4.2"
 RecipesBase = "0.7, 0.8, 1.0"
 RecursiveArrayTools = "2"
+RootedTrees = "1"
 julia = "1"
 
 [extras]

src/DiffEqDevTools.jl

@@ -40,7 +40,7 @@ export test_convergence, analyticless_test_convergence, appxtrue!, appxtrue
 export get_sample_errors
 
 #Plot Functions
-export stability_region
+export stability_region, residual_order_condition
 
 #Tableaus
 export deduce_Butcher_tableau
@@ -51,7 +51,7 @@ export constructEuler, constructKutta3, constructRK4, constructRK438Rule,
        constructLobattoIIICStar4, constructLobattoIIID2, constructLobattoIIID4,
        constructRadauIA3, constructRadauIA5,
        constructRadauIIA3, constructRadauIIA5,
-       constructRalston, constructHeun, constructRKF5, constructBogakiShampine3,
+       constructRalston, constructRalston4, constructHeun, constructRKF5, constructBogakiShampine3,
        constructCashKarp, constructRKF8, constructDormandPrince8,
        constructMSRI1,constructFeagin10, constructFeagin12, constructFeagin14,
        constructDormandPrince8_64bit, constructRKF5, constructRungeFirst5,

src/ode_tableaus.jl

@@ -231,6 +231,35 @@ function constructRK438Rule(T::Type = Float64)
   return(ExplicitRKTableau(A,c,α,4))
 end
 
+"""
+Ralston's Order 4 method with minimum truncation error.
+"""
+function constructRalston4(T::Type = Float64)
+  sqrt5 = sqrt(convert(T, 5))
+  a21 = 4//10
+  a31 = (-2889 + 1428 * sqrt5) / 1024
+  a32 = (3785 - 1620 * sqrt5) / 1024
+  a41 = (-3365 + 2094 * sqrt5) / 6040
+  a42 = (-975 - 3046 * sqrt5) / 2552
+  a43 = (467040 + 203968 * sqrt5) / 240845
+  A = [0     0     0      0
+       a21   0     0      0
+       a31   a32   0      0
+       a41   a42   a43    0]
+  b2 = 4//10
+  b3 = (14 - 3 * sqrt5) / 16
+  c = [0, b2, b3, 1]
+  b1 = (263 + 24 * sqrt5) / 1812
+  b2 = (125 - 1000 * sqrt5) / 3828
+  b3 = 1024 * (3346 + 1623 * sqrt5) / 5924787
+  b4 = (30 - 4 * sqrt5) / 123
+  α = [b1, b2, b3, b4]
+  A = map(T,A)
+  α = map(T,α)
+  c = map(T,c)
+  return ExplicitRKTableau(A,c,α,4)
+end
+
 """
 Explicit SSP method of order 2 using 2 stages.
 """

src/tableau_info.jl

@@ -1,4 +1,4 @@
-using NLsolve, LinearAlgebra
+using NLsolve, LinearAlgebra, RootedTrees
 """
 `Base.length(tab::ODERKTableau)`
 
@@ -29,3 +29,16 @@ function stability_region(tab::ODERKTableau; initial_guess=-3.0)
   sol = nlsolve(residual!, [initial_guess])
   sol.zero[1]
 end
+
+function RootedTrees.residual_order_condition(tab::ODERKTableau, order::Int, reducer=nothing, mapper=x->x^2)
+  if reducer === nothing
+    resid = map(RootedTreeIterator(order)) do t
+      residual_order_condition(t, tab.A, tab.α, tab.c)
+    end
+  else
+    resid = mapreduce(reducer, RootedTreeIterator(order)) do t
+      mapper(residual_order_condition(t, tab.A, tab.α, tab.c))
+    end
+  end
+  return resid
+end

test/analyticless_convergence_tests.jl

@@ -1,5 +1,5 @@
 using OrdinaryDiffEq, ParameterizedFunctions, Test, Random
-
+using ParameterizedFunctions.ModelingToolkit # macro hygiene
 f = @ode_def LotkaVolterra begin
   dx = 1.5x - x*y
   dy = -3y + x*y

test/ode_tableau_convergence_tests.jl

@@ -18,6 +18,8 @@ dts = 1 .//2 .^(8:-1:4)
 testTol = 0.3
 superduperbool = Vector{Bool}(undef, 2)
 
+@test all(i -> residual_order_condition(constructRalston4(), i, +, abs) < 10eps(1.0), 1:4)
+
 for i = 1:2 # 1 = num, 2 = ExplicitRK
   global dts
   if i>1
@@ -63,6 +65,10 @@ for i = 1:2 # 1 = num, 2 = ExplicitRK
   sim = test_convergence(dts,prob,tabalg)
   @test abs(sim.𝒪est[:l∞]-4) < testTol
 
+  tabalg = ExplicitRK(tableau=constructRalston4())
+  sim = test_convergence(dts,prob,tabalg)
+  @test abs(sim.𝒪est[:l∞]-4) < testTol
+
   tabalg = ExplicitRK(tableau=constructSSPRK104())
   sim = test_convergence(dts,prob,tabalg)
   @test abs(sim.𝒪est[:l∞]-4) < testTol