From e4a7d2c4dd986ef1e725acdcde746a0fa45e59fb Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 22 Jun 2023 08:23:27 +0000 Subject: [PATCH 1/5] Added ```FineRKN4``` --- docs/src/dynamical/nystrom.md | 1 + src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 1 + src/algorithms.jl | 22 ++++ src/caches/rkn_caches.jl | 42 +++++++ src/perform_step/rkn_perform_step.jl | 114 +++++++++++++++++- src/tableaus/rkn_tableaus.jl | 101 ++++++++++++++++ .../partitioned_methods_tests.jl | 10 ++ 8 files changed, 290 insertions(+), 3 deletions(-) diff --git a/docs/src/dynamical/nystrom.md b/docs/src/dynamical/nystrom.md index dbedd6f5d4..ebea50d1e8 100644 --- a/docs/src/dynamical/nystrom.md +++ b/docs/src/dynamical/nystrom.md @@ -5,6 +5,7 @@ IRKN3 Nystrom4 Nystrom4VelocityIndependent IRKN4 +FineRKN4 FineRKN5 Nystrom5VelocityIndependent DPRKN6 diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index f64921afb7..5c0b57773b 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -437,7 +437,7 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, export SplitEuler -export Nystrom4, FineRKN5, Nystrom4VelocityIndependent, +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 0fec79643c..4c72d6ae85 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -502,6 +502,7 @@ alg_order(alg::SofSpa10) = 10 alg_order(alg::IRKN3) = 3 alg_order(alg::Nystrom4) = 4 +alg_order(alg::FineRKN4) = 4 alg_order(alg::FineRKN5) = 5 alg_order(alg::Nystrom4VelocityIndependent) = 4 alg_order(alg::IRKN4) = 4 diff --git a/src/algorithms.jl b/src/algorithms.jl index 40a1d6542c..0b588e8271 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5251,6 +5251,28 @@ E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equatio """ struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end +""" + FineRKN4() + +A 4th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. +In particular, this method allows the acceleration equation to depend on the velocity. + +## References +``` +@article{fine1987low, + title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, + author={Fine, Jerry Michael}, + journal={Computing}, + volume={38}, + number={4}, + pages={281--297}, + year={1987}, + publisher={Springer} +} +``` +""" +struct FineRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + """ FineRKN5() diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 074c7344c5..7ffa67bf23 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -36,6 +36,48 @@ end # alg_cache(alg::Nystrom4,u,rate_prototype,::Type{uEltypeNoUnits},::Type{uBottomEltypeNoUnits},::Type{tTypeNoUnits},uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false}) where {uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits} = Nystrom4ConstantCache(constvalue(uBottomEltypeNoUnits),constvalue(tTypeNoUnits)) +@cache struct FineRKN4Cache{uType, rateType, reducedRateType, uNoUnitsType, TabType} <: + OrdinaryDiffEqMutableCache + u::uType + uprev::uType + fsalfirst::rateType + k2::reducedRateType + k3::reducedRateType + k4::reducedRateType + k5::reducedRateType + k::rateType + utilde::uType + tmp::uType + atmp::uNoUnitsType + tab::TabType +end + +function alg_cache(alg::FineRKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + reduced_rate_prototype = rate_prototype.x[2] + tab = FineRKN4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) + k1 = zero(rate_prototype) + k2 = zero(reduced_rate_prototype) + k3 = zero(reduced_rate_prototype) + k4 = zero(reduced_rate_prototype) + k5 = zero(reduced_rate_prototype) + k = zero(rate_prototype) + utilde = zero(u) + atmp = similar(u, uEltypeNoUnits) + recursivefill!(atmp, false) + tmp = zero(u) + FineRKN4Cache(u, uprev, k1, k2, k3, k4, k5, k, utilde, tmp, atmp, tab) +end + +function alg_cache(alg::FineRKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + FineRKN4ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) +end + @cache struct FineRKN5Cache{uType, rateType, reducedRateType, uNoUnitsType, TabType} <: OrdinaryDiffEqMutableCache u::uType diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index 6995e655b9..4d93f3c9c6 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -4,7 +4,8 @@ ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ ## y'₁ = y'₀ + h∑bᵢk'ᵢ -const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN5ConstantCache, +const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN4ConstantCache, + FineRKN5ConstantCache, Nystrom4VelocityIndependentConstantCache, Nystrom5VelocityIndependentConstantCache, IRKN3ConstantCache, IRKN4ConstantCache, @@ -25,7 +26,7 @@ function initialize!(integrator, cache::NystromCCDefaultInitialization) integrator.fsalfirst = ArrayPartition((kdu, ku)) end -const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN5Cache, +const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN4Cache, FineRKN5Cache, Nystrom4VelocityIndependentCache, Nystrom5VelocityIndependentCache, IRKN3Cache, IRKN4Cache, @@ -122,6 +123,115 @@ end integrator.stats.nf2 += 1 end +@muladd function perform_step!(integrator, cache::FineRKN4ConstantCache, + repeat_step = false) + @unpack t, dt, f, p = integrator + duprev, uprev = integrator.uprev.x + @unpack c2, c3, c4, c5, a21, a31, a32, a41, a43, a51, + a52, a53, a54, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, btilde1, btilde3, btilde4, btilde5, bptilde1, + bptilde3, bptilde4, bptilde5 = cache + k1 = integrator.fsalfirst.x[1] + + ku = uprev + dt * (c2 * duprev + dt * (a21 * k1)) + kdu = duprev + dt * (abar21 * k1) + + k2 = f.f1(kdu, ku, p, t + dt * c2) + ku = uprev + dt * (c3 * duprev + dt * (a31 * k1 + a32 * k2)) + kdu = duprev + dt * (abar31 * k1 + abar32 * k2) + + k3 = f.f1(kdu, ku, p, t + dt * c3) + ku = uprev + dt * (c4 * duprev + dt * (a41 * k1 + a43 * k3)) # a42 = 0 + kdu = duprev + dt * (abar41 * k1 + abar42 * k2 + abar43 * k3) + + k4 = f.f1(kdu, ku, p, t + dt * c4) + ku = uprev + dt * (c5 * duprev + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4)) + kdu = duprev + dt * (abar51 * k1 + abar52 * k2 + abar53 * k3 + abar54 * k4) + + k5 = f.f1(kdu, ku, p, t + dt * c5) + + u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # b2 = 0 + du = duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5) # bbar2 = 0 + + integrator.u = ArrayPartition((du, u)) + integrator.fsallast = ArrayPartition((f.f1(du, u, p, t + dt), f.f2(du, u, p, t + dt))) + integrator.stats.nf += 5 + integrator.stats.nf2 += 1 + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast + + if integrator.opts.adaptive + dtsq = dt^2 + uhat = dtsq * (btilde1 * k1 + btilde3 * k3 + btilde4 * k4 + btilde5 * k5) # btilde2 = 0 + duhat = dt * (bptilde1 * k1 + bptilde3 * k3 + bptilde4 * k4 + bptilde5 * k5) # bptilde2 = 0 + utilde = ArrayPartition((duhat, uhat)) + atmp = calculate_residuals(utilde, integrator.uprev, integrator.u, + integrator.opts.abstol, integrator.opts.reltol, + integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(atmp, t) + end +end + +@muladd function perform_step!(integrator, cache::FineRKN4Cache, repeat_step = false) + @unpack t, dt, f, p = integrator + du, u = integrator.u.x + duprev, uprev = integrator.uprev.x + @unpack tmp, atmp, fsalfirst, k2, k3, k4, k5, k, utilde = cache + @unpack c2, c3, c4, c5, a21, a31, a32, a41, a43, a51, + a52, a53, a54, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, btilde1, btilde3, btilde4, btilde5, bptilde1, + bptilde3, bptilde4, bptilde5 = cache.tab + kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2] + uidx = eachindex(integrator.uprev.x[2]) + k1 = integrator.fsalfirst.x[1] + + @.. broadcast=false ku=uprev + dt * (c2 * duprev + dt * (a21 * k1)) + @.. broadcast=false kdu=duprev + dt * (abar21 * k1) + + f.f1(k2, kdu, ku, p, t + dt * c2) + @.. broadcast=false ku=uprev + dt * (c3 * duprev + dt * (a31 * k1 + a32 * k2)) + @.. broadcast=false kdu=duprev + dt * (abar31 * k1 + abar32 * k2) + + f.f1(k3, kdu, ku, p, t + dt * c3) + @.. broadcast=false ku=uprev + + dt * (c4 * duprev + dt * (a41 * k1 + a43 * k3)) # a42 = 0 + @.. broadcast=false kdu=duprev + dt * (abar41 * k1 + abar42 * k2 + abar43 * k3) + + f.f1(k4, kdu, ku, p, t + dt * c4) + @.. broadcast=false ku=uprev + + dt * + (c5 * duprev + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4)) + @.. broadcast=false kdu=duprev + + dt * (abar51 * k1 + abar52 * k2 + abar53 * k3 + abar54 * k4) + + f.f1(k5, kdu, ku, p, t + dt * c5) + @.. broadcast=false u=uprev + + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # b2 = 0 + @.. broadcast=false du=duprev + + dt * + (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5) # bbar2 = 0 + + f.f1(k.x[1], du, u, p, t + dt) + f.f2(k.x[2], du, u, p, t + dt) + integrator.stats.nf += 5 + integrator.stats.nf2 += 1 + if integrator.opts.adaptive + duhat, uhat = utilde.x + dtsq = dt^2 + @.. broadcast=false uhat=dtsq * + (btilde1 * k1 + btilde3 * k3 + btilde4 * k4 + + btilde5 * k5) # btilde2 = 0 + @.. broadcast=false duhat=dt * + (bptilde1 * k1 + bptilde3 * k3 + bptilde4 * k4 + + bptilde5 * k5) # bptilde2 = 0 + + calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u, + integrator.opts.abstol, integrator.opts.reltol, + integrator.opts.internalnorm, t) + integrator.EEst = integrator.opts.internalnorm(atmp, t) + end +end + @muladd function perform_step!(integrator, cache::FineRKN5ConstantCache, repeat_step = false) @unpack t, dt, f, p = integrator diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index aa2e7d7d44..d4e89d2a15 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -1,3 +1,104 @@ +struct FineRKN4ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache + c1::T2 + c2::T2 + c3::T2 + c4::T2 + c5::T2 + a21::T + a31::T + a32::T + a41::T + #a42::T + a43::T + a51::T + a52::T + a53::T + a54::T + abar21::T + abar31::T + abar32::T + abar41::T + abar42::T + abar43::T + abar51::T + abar52::T + abar53::T + abar54::T + b1::T + #b2::T + b3::T + b4::T + b5::T + bbar1::T + #bbar2::T + bbar3::T + bbar4::T + bbar5::T + btilde1::T + #btilde2::T + btilde3::T + btilde4::T + btilde5::T + bptilde1::T + #bptilde2::T + bptilde3::T + bptilde4::T + bptilde5::T +end + +function FineRKN4ConstantCache(T::Type, T2::Type) + c1 = convert(T2, 1 // 1) + c2 = convert(T2, 2 // 9) + c3 = convert(T2, 1 // 3) + c4 = convert(T2, 3 // 4) + c5 = convert(T2, 1 // 1) + a21 = convert(T, 2 // 81) + a31 = convert(T, 1 // 36) + a32 = convert(T, 1 // 36) + a41 = convert(T, 9 // 128) + #a42 = convert(T, 0 // 1) + a43 = convert(T, 27 // 128) + a51 = convert(T, 11 // 60) + a52 = convert(T, -3 // 20) + a53 = convert(T, 9 // 25) + a54 = convert(T, 8 // 75) + abar21 = convert(T, 2 // 9) + abar31 = convert(T, 1 // 12) + abar32 = convert(T, 1 // 4) + abar41 = convert(T, 69 // 128) + abar42 = convert(T, -243 // 128) + abar43 = convert(T, 135 // 64) + abar51 = convert(T, -17 // 12) + abar52 = convert(T, 27 // 4) + abar53 = convert(T, -27 // 5) + abar54 = convert(T, 16 // 15) + b1 = convert(T, 19 // 180) + #b2 = convert(T, 0 // 1) + b3 = convert(T, 63 // 200) + b4 = convert(T, 16 // 225) + b5 = convert(T, 1 // 120) + bbar1 = convert(T, 1 // 9) + #bbar2 = convert(T, 0 // 1) + bbar3 = convert(T, 9 // 20) + bbar4 = convert(T, 16 // 45) + bbar5 = convert(T, 1 // 12) + btilde1 = convert(T, 25 // 1116) + #btilde2 = convert(T, 0 // 1) + btilde3 = convert(T, -63 // 1240) + btilde4 = convert(T, 64 // 1395) + btilde5 = convert(T, -13 // 744) + bptilde1 = convert(T, 2 // 125) + #bptilde2 = convert(T, 0 // 1) + bptilde3 = convert(T, -27 // 625) + bptilde4 = convert(T, 32 // 625) + bptilde5 = convert(T, -3 // 125) + FineRKN4ConstantCache(c1, c2, c3, c4, c5, a21, a31, a32, a41, a43, a51, + a52, a53, a54, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, btilde1, + btilde3, btilde4, btilde5, bptilde1, + bptilde3, bptilde4, bptilde5) +end + struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache c1::T2 c2::T2 diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 3e1a88f9f1..11dd10566e 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -128,6 +128,9 @@ dts = 1.0 ./ 2.0 .^ (5:-1:0) sim = test_convergence(dts, prob, Nystrom5VelocityIndependent(), dense_errors = true) @test sim.𝒪est[:l2]≈5 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 +sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true) +@test sim.𝒪est[:l2]≈5 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true) @test sim.𝒪est[:l2]≈6 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -170,6 +173,8 @@ sim = test_convergence(dts, prob_big, ERKN7(), dense_errors = true) @test sim.𝒪est[:L2]≈4 rtol=1e-1 # Adaptive methods regression test +sol = solve(prob, FineRKN4()) +@test length(sol.u) < 16 sol = solve(prob, FineRKN5()) @test length(sol.u) < 14 sol = solve(prob, DPRKN4()) @@ -305,6 +310,9 @@ dts = 1.0 ./ 2.0 .^ (5:-1:0) sim = test_convergence(dts, prob, Nystrom5VelocityIndependent(), dense_errors = true) @test sim.𝒪est[:l2]≈5 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 +sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true) +@test sim.𝒪est[:l2]≈5 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true) @test sim.𝒪est[:l2]≈6 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -347,6 +355,8 @@ sim = test_convergence(dts, prob_big, ERKN7(), dense_errors = true) @test sim.𝒪est[:L2]≈4 rtol=1e-1 # Adaptive methods regression test +sol = solve(prob, FineRKN4()) +@test length(sol.u) < 16 sol = solve(prob, FineRKN5()) @test length(sol.u) < 14 sol = solve(prob, DPRKN4()) From 80c72a3acc79df63f98ff9a853c9666bb3bb704e Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 22 Jun 2023 09:31:05 +0000 Subject: [PATCH 2/5] Added convergence tests for the damped oscillator. --- .../partitioned_methods_tests.jl | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 11dd10566e..1a62aa5ff0 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -378,6 +378,75 @@ sol = solve(prob, ERKN5(), reltol = 1e-8) sol = solve(prob, ERKN7(), reltol = 1e-8) @test length(sol.u) < 38 +# Testing generalized Runge-Kutte-Nyström methods on velocity dependend ODEs with the damped oscillator +println("In Place") + +# Damped oscillator +prob = ODEProblem(DynamicalODEFunction{false}((du, u, p, t) -> -u - 0.5 * du, + (du, u, p, t) -> du, + analytic = (du0_u0, p, t) -> OrdinaryDiffEq.SciMLBase.ArrayPartition([ + exp(-t / 4) / 15 * (15 * du0_u0[1] * cos(sqrt(15) * t / 4) - + sqrt(15) * (du0_u0[1] + 4 * du0_u0[2]) * sin(sqrt(15) * t / 4)), + ], # du + [ + exp(-t / 4) / 15 * (15 * du0_u0[2] * cos(sqrt(15) * t / 4) + + sqrt(15) * (4 * du0_u0[1] + du0_u0[2]) * sin(sqrt(15) * t / 4)), + ])), + OrdinaryDiffEq.SciMLBase.ArrayPartition([0.0], [1.0]), # du0, u0 + (0.0, 10.0), # tspan + DiffEqBase.NullParameters(), # p + SecondOrderODEProblem{false}()) +dts = 1.0 ./ 2.0 .^ (5:-1:0) +sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true); +@test sim.𝒪est[:l2]≈4 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 +sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true); +@test sim.𝒪est[:l2]≈4 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 +sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true); +@test sim.𝒪est[:l2]≈5 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 + +# Adaptive methods regression test +sol = solve(prob, FineRKN4()) +@test length(sol.u) < 28 +sol = solve(prob, FineRKN5()) +@test length(sol.u) < 20 + +println("Out of Place") +# Damped oscillator +prob = ODEProblem(DynamicalODEFunction{true}((d_du, du, u, p, t) -> @.(d_du=-u - 0.5 * du), + (d_u, du, u, p, t) -> d_u .= du, + analytic = (du0_u0, p, t) -> OrdinaryDiffEq.SciMLBase.ArrayPartition([ + exp(-t / 4) / 15 * (15 * du0_u0[1] * cos(sqrt(15) * t / 4) - + sqrt(15) * (du0_u0[1] + 4 * du0_u0[2]) * sin(sqrt(15) * t / 4)), + ], # du + [ + exp(-t / 4) / 15 * (15 * du0_u0[2] * cos(sqrt(15) * t / 4) + + sqrt(15) * (4 * du0_u0[1] + du0_u0[2]) * sin(sqrt(15) * t / 4)), + ])), + OrdinaryDiffEq.SciMLBase.ArrayPartition([0.0], [1.0]), # du0, u0 + (0.0, 10.0), # tspan + DiffEqBase.NullParameters(), # p + SecondOrderODEProblem{false}()) + +dts = 1.0 ./ 2.0 .^ (5:-1:0) +sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true); +@test sim.𝒪est[:l2]≈4 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 +sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true); +@test sim.𝒪est[:l2]≈4 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 +sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true); +@test sim.𝒪est[:l2]≈5 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 + +# Adaptive methods regression test +sol = solve(prob, FineRKN4()) +@test length(sol.u) < 28 +sol = solve(prob, FineRKN5()) +@test length(sol.u) < 20 + # Compare in-place and out-of-place versions function damped_oscillator(du, u, p, t) return -u - 0.5 * du @@ -409,6 +478,26 @@ end @test abs(sol_i.destats.nf - 4 * sol_i.destats.naccept) < 4 end + @testset "FineRKN4" begin + alg = FineRKN4() + dt = 0.5 + # fixed time step + sol_i = solve(ode_i, alg, adaptive = false, dt = dt); + sol_o = solve(ode_o, alg, adaptive = false, dt = dt); + @test sol_i.t ≈ sol_o.t + @test sol_i.u ≈ sol_o.u + @test sol_i.destats.nf == sol_o.destats.nf + @test sol_i.destats.nf2 == sol_o.destats.nf2 + @test sol_i.destats.naccept == sol_o.destats.naccept + @test 19 <= sol_i.destats.naccept <= 21 + @test abs(sol_i.destats.nf - 5 * sol_i.destats.naccept) < 4 + # adaptive time step + sol_i = solve(ode_i, alg) + sol_o = solve(ode_o, alg) + @test sol_i.t ≈ sol_o.t + @test sol_i.u ≈ sol_o.u + end + @testset "FineRKN5" begin alg = FineRKN5() dt = 0.5 From 045654bb28c2c2e9b74e265bc9a1f796c1af206a Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 22 Jun 2023 10:06:58 +0000 Subject: [PATCH 3/5] typo --- test/algconvergence/partitioned_methods_tests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 1a62aa5ff0..097f167b74 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -396,6 +396,7 @@ prob = ODEProblem(DynamicalODEFunction{false}((du, u, p, t) -> -u - 0.5 * du, (0.0, 10.0), # tspan DiffEqBase.NullParameters(), # p SecondOrderODEProblem{false}()) + dts = 1.0 ./ 2.0 .^ (5:-1:0) sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true); @test sim.𝒪est[:l2]≈4 rtol=1e-1 From 059682d56a2d225302b750086550d245976974e8 Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 22 Jun 2023 10:14:22 +0000 Subject: [PATCH 4/5] typo --- .../partitioned_methods_tests.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 097f167b74..eee3ba10db 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -396,15 +396,15 @@ prob = ODEProblem(DynamicalODEFunction{false}((du, u, p, t) -> -u - 0.5 * du, (0.0, 10.0), # tspan DiffEqBase.NullParameters(), # p SecondOrderODEProblem{false}()) - + dts = 1.0 ./ 2.0 .^ (5:-1:0) -sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true); +sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 -sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true); +sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 -sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true); +sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true) @test sim.𝒪est[:l2]≈5 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -432,13 +432,13 @@ prob = ODEProblem(DynamicalODEFunction{true}((d_du, du, u, p, t) -> @.(d_du=-u - SecondOrderODEProblem{false}()) dts = 1.0 ./ 2.0 .^ (5:-1:0) -sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true); +sim = test_convergence(dts, prob, Nystrom4(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 -sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true); +sim = test_convergence(dts, prob, FineRKN4(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 -sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true); +sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true) @test sim.𝒪est[:l2]≈5 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -483,8 +483,8 @@ end alg = FineRKN4() dt = 0.5 # fixed time step - sol_i = solve(ode_i, alg, adaptive = false, dt = dt); - sol_o = solve(ode_o, alg, adaptive = false, dt = dt); + sol_i = solve(ode_i, alg, adaptive = false, dt = dt) + sol_o = solve(ode_o, alg, adaptive = false, dt = dt) @test sol_i.t ≈ sol_o.t @test sol_i.u ≈ sol_o.u @test sol_i.destats.nf == sol_o.destats.nf From 5cefbb8d5f3286d76b3999d1de7961641506d2e5 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Thu, 22 Jun 2023 17:10:31 +0200 Subject: [PATCH 5/5] format --- test/algconvergence/partitioned_methods_tests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index eee3ba10db..09444262bf 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -498,7 +498,7 @@ end @test sol_i.t ≈ sol_o.t @test sol_i.u ≈ sol_o.u end - + @testset "FineRKN5" begin alg = FineRKN5() dt = 0.5