finish impl of dp8 + DRY

src/dp5/helpers.jl

@@ -37,19 +37,6 @@ function error_estimation(solver)
     return err 
 end
 
-function estimate_second_derivative(solver, h)
-
-    der2 = mapreduce(+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k2, solver.k1, solver.y) do atoli, rtoli, f1i, f0i, yi
-        sk = atoli + rtoli*abs(yi)
-        ((f1i-f0i)/sk)^2
-    end
-
-    der2 = sqrt(der2)/h
-
-    return der2
-
-end
-
 function stiffness_detection!(solver, naccpt, h)
     if (mod(naccpt, solver.options.stiffness_test_activation_step) == 0) || (solver.vars.iasti > 0)
         #stnum = 0.0
@@ -63,7 +50,7 @@ function stiffness_detection!(solver, naccpt, h)
         end
 
         if stden > 0.0
-            solver.vars.hlamb = h*sqrt(stnum/stden)
+            solver.vars.hlamb = abs(h)*sqrt(stnum/stden)
         else
             solver.vars.hlamb = Inf
         end
@@ -82,22 +69,3 @@ function stiffness_detection!(solver, naccpt, h)
         end
     end
 end
-
-function euler_first_guess(solver, hmax, posneg)
-
-    dnf, dny = mapreduce(.+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k1, solver.y) do atoli, rtoli, f0i, yi
-        sk = atoli + rtoli*abs(yi)
-        abs(f0i/sk)^2, abs(yi/sk)^2 # dnf, dny
-    end
-
-
-    if (dnf <= 1.0e-10) || (dny <= 1.0e-10)
-        h = 1.0e-6
-    else
-        h = 0.01*sqrt(dny/dnf)
-    end
-    h = min(h, hmax)
-    h = h * Base.sign(posneg)
-
-    return h, dnf
-end

src/dp5/mod.jl

@@ -1,6 +1,6 @@
 module  DP5
 
-using ..DormandPrince: DormandPrince, DP5Solver, Vars, Consts, Options, Report
+using ..DormandPrince: DormandPrince, DP5Solver, Vars, Consts, Options, Report, hinit
 
 include("params.jl")
 include("helpers.jl")

src/dp5/solver.jl

@@ -194,43 +194,3 @@ function dopcor(
     )
 
 end
-
-function hinit(
-    solver, 
-    posneg, 
-    iord,
-    hmax
-    # f0 arg is k1 from dopcor
-    # f1 arg is k2 from dopcor
-    # y1 arg is k3 from dopcor
-)
-    #=
-    Compute a first guess for explicit euler as
-        h = 0.01 * norm (y0) / norm (f0)
-    the increment for explicit euler is small
-    compared to the solution
-    =#
-    h, dnf = euler_first_guess(solver, hmax, posneg)
-
-    ###### Perform an explicit step
-    #y1 = y + h*f0
-    #fcn(n, x+h, y1, f1)
-    # copyto!(solver.y1, solver.y + h*solver.k1)
-    solver.y1 .= solver.y .+ h .*solver.k1
-    solver.f(solver.vars.x + h, solver.k3, solver.k2)
-
-    ###### Estimate the second derivative of the solution
-    der2 = estimate_second_derivative(solver, h)
-
-    ##### Step size is computed such that
-    ##### H**IORD * MAX ( NORM(F0), NORM(F1), DER2 ) = 0.01
-    der12 = max(abs(der2), sqrt(dnf))
-    if der12 <= 1e-15
-        h1 = max(1.0e-6, abs(h)*1.0e-3)
-    else
-        h1 = (0.01/der12)^(1.0/iord)
-    end
-    h = min(100*abs(h), h1, hmax)
-    return h * Base.sign(posneg)
-
-end

src/dp8/helpers.jl

@@ -145,39 +145,24 @@ function error_estimation(solver)
     return err 
 end
 
-function estimate_second_derivative(solver, h)
-
-    der2 = mapreduce(+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k2, solver.k1, solver.y) do atoli, rtoli, f1i, f0i, yi
-        sk = atoli + rtoli*abs(yi)
-        ((f1i-f0i)/sk)^2
-    end
-
-    der2 = sqrt(der2)/h
-
-    return der2
-
-end
-
 function stiffness_detection!(solver, naccpt, h)
+
     if (mod(naccpt, solver.options.stiffness_test_activation_step) == 0) || (solver.vars.iasti > 0)
         #stnum = 0.0
         #stden = 0.0
 
-        stnum, stden = mapreduce(.+, solver.k2, solver.k6, solver.y1, solver.ysti) do k2i, k6i, y1i, ystii
-            #stnum = abs(k2i-k6i)^2 
-            #stden = abs(y1i-ystii)^2
-            abs(k2i-k6i)^2, abs(y1i-ystii)^2
-            # stnum, stden
+        stnum, stden = mapreduce(.+, solver.k3, solver.k4, solver.k5, solver.y1) do k3i, k4i, k5i, y1i
+            abs(k4i-k3i)^2, abs(k5i-y1i)^2
         end
 
         if stden > 0.0
-            solver.vars.hlamb = h*sqrt(stnum/stden)
+            solver.vars.hlamb = abs(h)*sqrt(stnum/stden)
         else
             solver.vars.hlamb = Inf
         end
 
 
-        if solver.vars.hlamb > 3.25
+        if solver.vars.hlamb > 6.1
             solver.vars.iasti += 1
             if solver.vars.iasti == 15
                 @debug "The problem seems to become stiff at $x" 
@@ -190,22 +175,3 @@ function stiffness_detection!(solver, naccpt, h)
         end
     end
 end
-
-function euler_first_guess(solver, hmax, posneg)
-
-    dnf, dny = mapreduce(.+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k1, solver.y) do atoli, rtoli, f0i, yi
-        sk = atoli + rtoli*abs(yi)
-        abs(f0i/sk)^2, abs(yi/sk)^2 # dnf, dny
-    end
-
-
-    if (dnf <= 1.0e-10) || (dny <= 1.0e-10)
-        h = 1.0e-6
-    else
-        h = 0.01*sqrt(dny/dnf)
-    end
-    h = min(h, hmax)
-    h = h * Base.sign(posneg)
-
-    return h, dnf
-end

src/dp8/mod.jl

@@ -1,6 +1,6 @@
 module  DP8
 
-using ..DormandPrince: DormandPrince, DP8Solver, Vars, Consts, Options, Report
+using ..DormandPrince: DormandPrince, DP8Solver, Vars, Consts, Options, Report, hinit
 
 include("params.jl")
 include("helpers.jl")

src/hinit.jl

@@ -0,0 +1,74 @@
+
+
+function euler_first_guess(solver::AbstractDPSolver{T}, hmax::T, posneg::T) where T
+
+    dnf, dny = mapreduce(.+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k1, solver.y) do atoli, rtoli, f0i, yi
+        sk = atoli + rtoli*abs(yi)
+        abs(f0i/sk)^2, abs(yi/sk)^2 # dnf, dny
+    end
+
+
+    if (dnf <= 1.0e-10) || (dny <= 1.0e-10)
+        h = 1.0e-6
+    else
+        h = 0.01*sqrt(dny/dnf)
+    end
+    h = min(h, hmax)
+    h = h * Base.sign(posneg)
+
+    return h, dnf
+end
+
+
+function estimate_second_derivative(solver::AbstractDPSolver{T}, h::T) where {T}
+
+    der2 = mapreduce(+, solver.consts.atol_iter, solver.consts.rtol_iter, solver.k2, solver.k1, solver.y) do atoli, rtoli, f1i, f0i, yi
+        sk = atoli + rtoli*abs(yi)
+        ((f1i-f0i)/sk)^2
+    end
+
+    der2 = sqrt(der2)/h
+
+    return der2
+
+end
+
+function hinit(
+    solver::AbstractDPSolver{T}, 
+    posneg::T, 
+    iord::Int,
+    hmax::T
+    # f0 arg is k1 from dopcor
+    # f1 arg is k2 from dopcor
+    # y1 arg is k3 from dopcor
+) where T
+    #=
+    Compute a first guess for explicit euler as
+        h = 0.01 * norm (y0) / norm (f0)
+    the increment for explicit euler is small
+    compared to the solution
+    =#
+    h, dnf = euler_first_guess(solver, hmax, posneg)
+
+    ###### Perform an explicit step
+    #y1 = y + h*f0
+    #fcn(n, x+h, y1, f1)
+    # copyto!(solver.y1, solver.y + h*solver.k1)
+    solver.y1 .= solver.y .+ h .*solver.k1
+    solver.f(solver.vars.x + h, solver.k3, solver.k2)
+
+    ###### Estimate the second derivative of the solution
+    der2 = estimate_second_derivative(solver, h)
+
+    ##### Step size is computed such that
+    ##### H**IORD * MAX ( NORM(F0), NORM(F1), DER2 ) = 0.01
+    der12 = max(abs(der2), sqrt(dnf))
+    if der12 <= 1e-15
+        h1 = max(1.0e-6, abs(h)*1.0e-3)
+    else
+        h1 = (0.01/der12)^(1.0/iord)
+    end
+    h = min(100*abs(h), h1, hmax)
+    return h * Base.sign(posneg)
+
+end

src/types.jl

@@ -68,6 +68,7 @@ end
 
 @kwdef mutable struct Vars{T <: Real}
     x::T = zero(T)
+    xph::T = zero(T)
     h::T = zero(T)
     facold::T = 1e-4
     iasti::Int = 0