From 6624457898d6ce6df09a58d12f8a520509c5843e Mon Sep 17 00:00:00 2001 From: Phillip Weinberg Date: Fri, 22 Dec 2023 09:11:46 -0500 Subject: [PATCH] finish impl of dp8 + DRY --- src/dp5/helpers.jl | 34 +-------------------- src/dp5/mod.jl | 2 +- src/dp5/solver.jl | 40 ------------------------- src/dp8/helpers.jl | 44 ++++----------------------- src/dp8/mod.jl | 2 +- src/hinit.jl | 74 ++++++++++++++++++++++++++++++++++++++++++++++ src/types.jl | 1 + 7 files changed, 83 insertions(+), 114 deletions(-) create mode 100644 src/hinit.jl diff --git a/src/dp5/helpers.jl b/src/dp5/helpers.jl index f0b02a1..403e1a6 100644 --- a/src/dp5/helpers.jl +++ b/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 \ No newline at end of file diff --git a/src/dp5/mod.jl b/src/dp5/mod.jl index 502799e..d1983e1 100644 --- a/src/dp5/mod.jl +++ b/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") diff --git a/src/dp5/solver.jl b/src/dp5/solver.jl index b70fead..b558375 100644 --- a/src/dp5/solver.jl +++ b/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 \ No newline at end of file diff --git a/src/dp8/helpers.jl b/src/dp8/helpers.jl index 6f7ec24..f8fde6e 100644 --- a/src/dp8/helpers.jl +++ b/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 \ No newline at end of file diff --git a/src/dp8/mod.jl b/src/dp8/mod.jl index 71ee7cf..0143398 100644 --- a/src/dp8/mod.jl +++ b/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") diff --git a/src/hinit.jl b/src/hinit.jl new file mode 100644 index 0000000..638eb32 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index 14fc3d4..dbb5d05 100644 --- a/src/types.jl +++ b/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