diff --git a/src/hinit.jl b/src/hinit.jl index 0bf0ba7..78b7911 100644 --- a/src/hinit.jl +++ b/src/hinit.jl @@ -7,8 +7,6 @@ function euler_first_guess(solver::AbstractDPSolver{T}, hmax::T, posneg::T) wher (abs(f0i)/sk)^2, (abs(yi)/sk)^2 # dnf, dny end - println("euler_first_guess: dnf = $dnf, dny = $dny") - if (dnf <= 1.0e-10) || (dny <= 1.0e-10) h = 1.0e-6 else @@ -48,7 +46,6 @@ function hinit( the increment for explicit euler is small compared to the solution =# - println("hinit: y = $(solver.y), k1 = $(solver.k1), k2 = $(solver.k2), k3 = $(solver.k3)") h, dnf = euler_first_guess(solver, hmax, posneg) ###### Perform an explicit step diff --git a/test/errors.jl b/test/errors.jl index 250d8d2..b93e5ef 100644 --- a/test/errors.jl +++ b/test/errors.jl @@ -1,11 +1,13 @@ using Test using DormandPrince: DP5Solver, + DP8Solver, Options, LARGER_NMAX_NEEDED, STEP_SIZE_BECOMES_TOO_SMALL using DormandPrince. DP5: dopri5 +using DormandPrince. DP8: dop853 @testset "Larger nmax needed" begin solver = DP5Solver( @@ -19,6 +21,17 @@ using DormandPrince. DP5: dopri5 h, report = dopri5(solver, 2/0.0001, 0.1, 0.0) @test report.idid == LARGER_NMAX_NEEDED + solver = DP8Solver( + stiff_fcn, + 0.0, # start at 0.0 + [0.0001] # delta + ; maximum_allowed_steps=1 + ) + + + h, report = dop853(solver, 2/0.0001, 0.1, 0.0) + @test report.idid == LARGER_NMAX_NEEDED + end @testset "Step size becomes too small" begin @@ -32,4 +45,15 @@ end h, report = dopri5(solver, 2/0.0001, 0.1, 0.0) @test report.idid == STEP_SIZE_BECOMES_TOO_SMALL + + solver = DP8Solver( + stiff_fcn, + 0.0, + [0.0001] + ; uround=10000 + ) + + + h, report = dop853(solver, 2/0.0001, 0.1, 0.0) + @test report.idid == STEP_SIZE_BECOMES_TOO_SMALL end diff --git a/test/stiff.jl b/test/stiff.jl index bdfdf3f..b22c5b2 100644 --- a/test/stiff.jl +++ b/test/stiff.jl @@ -1,5 +1,5 @@ using Test -using DormandPrince: DP5Solver, integrate +using DormandPrince: DP5Solver, DP8Solver, integrate function stiff_fcn(x, y, f) f[1] = y[1]^2 - y[1]^3 @@ -7,12 +7,13 @@ end @testset "Stiff ODE" begin - solver = DP5Solver( - stiff_fcn, - 0.0, # start at 0.0 - [0.0001] # initial value of delta - ) - - integrate(solver, 2/0.0001) - + for SolverType in [DP5Solver, DP8Solver] + solver = SolverType( + stiff_fcn, + 0.0, + [0.0001] + ) + + integrate(solver, 2/0.0001) + end end