From de350a3c0c948cf2cdd3daec1369085b80feabde Mon Sep 17 00:00:00 2001 From: Henry Date: Thu, 18 May 2023 11:39:27 +0200 Subject: [PATCH 01/15] Caches, Documentation and all relevant information for a new Nystrom5 Algorithm added. The perform_step! function is not yet written. MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Nystrom5 algorithm is based on "Fine, Jerry Michael. "Low order practical Runge-Kutta-Nyström methods." Computing 38.4 (1987): 281-297". --- src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 1 + src/algorithms.jl | 22 ++++++++++++++++++++ src/caches/rkn_caches.jl | 45 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 69 insertions(+), 1 deletion(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 0c89f1cfe6..76dc985fda 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -423,7 +423,7 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, export SplitEuler -export Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, +export Nystrom4, Nystrom5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7 export ROCK2, ROCK4, RKC, IRKC, ESERK4, ESERK5, SERK2 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 92cb2c49c3..4514b2acb3 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -487,6 +487,7 @@ alg_order(alg::SofSpa10) = 10 alg_order(alg::IRKN3) = 3 alg_order(alg::Nystrom4) = 4 +alg_order(alg::Nystrom4) = 5 alg_order(alg::Nystrom4VelocityIndependent) = 4 alg_order(alg::IRKN4) = 4 alg_order(alg::Nystrom5VelocityIndependent) = 5 diff --git a/src/algorithms.jl b/src/algorithms.jl index 2725aeb3cf..2457568229 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5246,6 +5246,28 @@ E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equatio """ struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end +""" + Nystrom5 + +A 5th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. + +In case the ODE Problem is not dependent on the first derivative consider using +[`Nystrom5VelocityIndependent`](@ref) to increase performance. + +## References +@article{fine1987low, + title={Low order practical Runge-Kutta-Nystr{\"o}m methods}, + author={Fine, Jerry Michael}, + journal={Computing}, + volume={38}, + number={4}, + pages={281--297}, + year={1987}, + publisher={Springer} +} +""" +struct Nystrom5 <: OrdinaryDiffEqPartitionedAlgorithm end + """ Nystrom4VelocityIdependent diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index dd39b3321f..d04ab79829 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -34,8 +34,53 @@ function alg_cache(alg::Nystrom4, u, rate_prototype, ::Type{uEltypeNoUnits}, Nystrom4ConstantCache() 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 Nystrom5Cache{uType, rateType, reducedRateType} <: OrdinaryDiffEqMutableCache + u::uType + uprev::uType + fsalfirst::rateType + k2::reducedRateType + k3::reducedRateType + k4::reducedRateType + k5::reducedRateType + k6::reducedRateType + k7::reducedRateType + k::rateType + tmp::uType +end + +# struct Nystrom5ConstantCache <: OrdinaryDiffEqConstantCache end + +function alg_cache(alg::Nystrom5, 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] + k1 = zero(rate_prototype) + k2 = zero(reduced_rate_prototype) + k3 = zero(reduced_rate_prototype) + k4 = zero(reduced_rate_prototype) + k5 = zero(reduced_rate_prototype) + k6 = zero(reduced_rate_prototype) + k7 = zero(reduced_rate_prototype) + k = zero(rate_prototype) + tmp = zero(u) + Nystrom5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp) +end + +struct Nystrom5ConstantCache <: OrdinaryDiffEqConstantCache end + +function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, + ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, + dt, reltol, p, calck, + ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} + Nystrom5ConstantCache() +end + + + @cache struct Nystrom4VelocityIndependentCache{uType, rateType, reducedRateType} <: OrdinaryDiffEqMutableCache u::uType From fa44553160c104801315e1d6869fa0a351830f41 Mon Sep 17 00:00:00 2001 From: Henry Date: Tue, 23 May 2023 10:45:06 +0200 Subject: [PATCH 02/15] Added tableau structure for Nystrom5, values not yet added --- src/caches/rkn_caches.jl | 4 - src/tableaus/rkn_tableaus.jl | 202 +++++++++++++++++++++++++++++++++++ 2 files changed, 202 insertions(+), 4 deletions(-) diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index d04ab79829..dd9a0d3748 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -51,8 +51,6 @@ end tmp::uType end -# struct Nystrom5ConstantCache <: OrdinaryDiffEqConstantCache end - function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, @@ -70,8 +68,6 @@ function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, Nystrom5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp) end -struct Nystrom5ConstantCache <: OrdinaryDiffEqConstantCache end - function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 725ab9b256..9962875788 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -1,3 +1,205 @@ +struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache + c1::T2 + c2::T2 + c3::T2 + c4::T2 + c5::T2 + c6::T2 + c7::T2 + a21::T + a31::T + a32::T + a41::T + a42::T + a43::T + a51::T + a52::T + a53::T + a54::T + a61::T + a62::T + a63::T + a64::T + a65::T + a71::T + a72::T + a73::T + a74::T + a75::T + a76::T + abar21::T + abar31::T + abar32::T + abar41::T + abar42::T + abar43::T + abar51::T + abar52::T + abar53::T + abar54::T + abar61::T + abar62::T + abar63::T + abar64::T + abar65::T + abar71::T + abar72::T + abar73::T + abar74::T + abar75::T + abar76::T + b1::T + b2::T + b3::T + b4::T + b5::T + b6::T + b7::T + bbar1::T + bbar2::T + bbar3::T + bbar4::T + bbar5::T + bbar6::T + bbar7::T +end + +function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats}) + c1 = convert(T2, \\ ) + c2 = convert(T2, \\ ) + c3 = convert(T2, \\ ) + c4 = convert(T2, \\ ) + c5 = convert(T2, \\ ) + c6 = convert(T2, \\ ) + c7 = convert(T2, \\ ) + a21 = convert(T, \\ ) + a31 = convert(T, \\ ) + a32 = convert(T, \\ ) + a41 = convert(T, \\ ) + a42 = convert(T, \\ ) + a43 = convert(T, \\ ) + a51 = convert(T, \\ ) + a52 = convert(T, \\ ) + a53 = convert(T, \\ ) + a54 = convert(T, \\ ) + a61 = convert(T, \\ ) + a62 = convert(T, \\ ) + a63 = convert(T, \\ ) + a64 = convert(T, \\ ) + a65 = convert(T, \\ ) + a71 = convert(T, \\ ) + a72 = convert(T, \\ ) + a73 = convert(T, \\ ) + a74 = convert(T, \\ ) + a75 = convert(T, \\ ) + a76 = convert(T, \\ ) + abar21 = convert(T, \\ ) + abar31 = convert(T, \\ ) + abar32 = convert(T, \\ ) + abar41 = convert(T, \\ ) + abar42 = convert(T, \\ ) + abar43 = convert(T, \\ ) + abar51 = convert(T, \\ ) + abar52 = convert(T, \\ ) + abar53 = convert(T, \\ ) + abar54 = convert(T, \\ ) + abar61 = convert(T, \\ ) + abar62 = convert(T, \\ ) + abar63 = convert(T, \\ ) + abar64 = convert(T, \\ ) + abar65 = convert(T, \\ ) + abar71 = convert(T, \\ ) + abar72 = convert(T, \\ ) + abar73 = convert(T, \\ ) + abar74 = convert(T, \\ ) + abar75 = convert(T, \\ ) + abar76 = convert(T, \\ ) + b1 = convert(T, \\ ) + b2 = convert(T, \\ ) + b3 = convert(T, \\ ) + b4 = convert(T, \\ ) + b5 = convert(T, \\ ) + b6 = convert(T, \\ ) + b7 = convert(T, \\ ) + bbar1 = convert(T, \\ ) + bbar2 = convert(T, \\ ) + bbar3 = convert(T, \\ ) + bbar4 = convert(T, \\ ) + bbar5 = convert(T, \\ ) + bbar6 = convert(T, \\ ) + bbar7 = convert(T, \\ ) + Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) +end + +function Nystrom5ConstantCache(T::Type, T2::Type) + c1 = convert(T2, \\ ) + c2 = convert(T2, \\ ) + c3 = convert(T2, \\ ) + c4 = convert(T2, \\ ) + c5 = convert(T2, \\ ) + c6 = convert(T2, \\ ) + c7 = convert(T2, \\ ) + a21 = convert(T, \\ ) + a31 = convert(T, \\ ) + a32 = convert(T, \\ ) + a41 = convert(T, \\ ) + a42 = convert(T, \\ ) + a43 = convert(T, \\ ) + a51 = convert(T, \\ ) + a52 = convert(T, \\ ) + a53 = convert(T, \\ ) + a54 = convert(T, \\ ) + a61 = convert(T, \\ ) + a62 = convert(T, \\ ) + a63 = convert(T, \\ ) + a64 = convert(T, \\ ) + a65 = convert(T, \\ ) + a71 = convert(T, \\ ) + a72 = convert(T, \\ ) + a73 = convert(T, \\ ) + a74 = convert(T, \\ ) + a75 = convert(T, \\ ) + a76 = convert(T, \\ ) + abar21 = convert(T, \\ ) + abar31 = convert(T, \\ ) + abar32 = convert(T, \\ ) + abar41 = convert(T, \\ ) + abar42 = convert(T, \\ ) + abar43 = convert(T, \\ ) + abar51 = convert(T, \\ ) + abar52 = convert(T, \\ ) + abar53 = convert(T, \\ ) + abar54 = convert(T, \\ ) + abar61 = convert(T, \\ ) + abar62 = convert(T, \\ ) + abar63 = convert(T, \\ ) + abar64 = convert(T, \\ ) + abar65 = convert(T, \\ ) + abar71 = convert(T, \\ ) + abar72 = convert(T, \\ ) + abar73 = convert(T, \\ ) + abar74 = convert(T, \\ ) + abar75 = convert(T, \\ ) + abar76 = convert(T, \\ ) + b1 = convert(T, \\ ) + b2 = convert(T, \\ ) + b3 = convert(T, \\ ) + b4 = convert(T, \\ ) + b5 = convert(T, \\ ) + b6 = convert(T, \\ ) + b7 = convert(T, \\ ) + bbar1 = convert(T, \\ ) + bbar2 = convert(T, \\ ) + bbar3 = convert(T, \\ ) + bbar4 = convert(T, \\ ) + bbar5 = convert(T, \\ ) + bbar6 = convert(T, \\ ) + bbar7 = convert(T, \\ ) + Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) +end + +################################################################################# + struct IRKN3ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache bconst1::T bconst2::T From 5d259360b59582f06fe536a895be2d2e4578c5e6 Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 26 May 2023 09:29:51 +0200 Subject: [PATCH 03/15] Perform_set and tableau for Nystrom5 added. Nystrom5 is working on in place and out of place ODEs however the convergence test is not passt for order >1. --- src/caches/rkn_caches.jl | 10 +- src/perform_step/rkn_perform_step.jl | 93 ++++++- src/tableaus/rkn_tableaus.jl | 252 +++++++++--------- .../partitioned_methods_tests.jl | 6 + 4 files changed, 227 insertions(+), 134 deletions(-) diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index dd9a0d3748..39c4b8d630 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -34,13 +34,13 @@ function alg_cache(alg::Nystrom4, u, rate_prototype, ::Type{uEltypeNoUnits}, Nystrom4ConstantCache() 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 Nystrom5Cache{uType, rateType, reducedRateType} <: OrdinaryDiffEqMutableCache +@cache struct Nystrom5Cache{uType, rateType, reducedRateType, TabType} <: OrdinaryDiffEqMutableCache u::uType uprev::uType fsalfirst::rateType + #k1::rateType k2::reducedRateType k3::reducedRateType k4::reducedRateType @@ -49,6 +49,7 @@ end k7::reducedRateType k::rateType tmp::uType + tab::TabType end function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, @@ -56,6 +57,7 @@ function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} reduced_rate_prototype = rate_prototype.x[2] + tab = Nystrom5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) k1 = zero(rate_prototype) k2 = zero(reduced_rate_prototype) k3 = zero(reduced_rate_prototype) @@ -65,14 +67,14 @@ function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, k7 = zero(reduced_rate_prototype) k = zero(rate_prototype) tmp = zero(u) - Nystrom5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp) + Nystrom5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp, tab) end function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - Nystrom5ConstantCache() + Nystrom5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index 75ea187d10..19f2cd53f6 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -4,7 +4,7 @@ ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ ## y'₁ = y'₀ + h∑bᵢk'ᵢ -const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, +const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, Nystrom5ConstantCache, Nystrom4VelocityIndependentConstantCache, Nystrom5VelocityIndependentConstantCache, IRKN3ConstantCache, IRKN4ConstantCache, @@ -25,7 +25,7 @@ function initialize!(integrator, cache::NystromCCDefaultInitialization) integrator.fsalfirst = ArrayPartition((kdu, ku)) end -const NystromDefaultInitialization = Union{Nystrom4Cache, +const NystromDefaultInitialization = Union{Nystrom4Cache, Nystrom5Cache, Nystrom4VelocityIndependentCache, Nystrom5VelocityIndependentCache, IRKN3Cache, IRKN4Cache, @@ -67,7 +67,7 @@ end kdu = duprev + halfdt * k₁ k₂ = f.f1(kdu, ku, p, ttmp) - ku = uprev + halfdt * duprev + eighth_dtsq * k₁ + ku = uprev + halfdt * duprev + eighth_dtsq * k₁ kdu = duprev + halfdt * k₂ k₃ = f.f1(kdu, ku, p, ttmp) @@ -105,7 +105,7 @@ end @.. broadcast=false kdu=duprev + halfdt * k₁ f.f1(k₂, kdu, ku, p, ttmp) - @.. broadcast=false ku=uprev + halfdt * duprev + eighth_dtsq * k₁ + @.. broadcast=false ku=uprev + halfdt * duprev + eighth_dtsq * k₁ @.. broadcast=false kdu=duprev + halfdt * k₂ f.f1(k₃, kdu, ku, p, ttmp) @@ -122,6 +122,91 @@ end integrator.stats.nf2 += 1 end +@muladd function perform_step!(integrator, cache::Nystrom5ConstantCache, + repeat_step = false) + @unpack t, dt, f, p = integrator + duprev, uprev = integrator.uprev.x + @unpack c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b3, b4, b5, b6, b7, bbar1, bbar3, bbar4, bbar5, bbar6, bbar7 = cache + k1 = integrator.fsalfirst.x[1] + + ku = uprev + (c2 * dt * duprev) + (dt^2 * 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 + a42 * k2 + a43 * k3)) + 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(duprev, ku, p, t + dt * c5) + ku = uprev + dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) + kdu = duprev + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) + + k6 = f.f1(duprev, ku, p, t + dt * c6) + ku = uprev + dt * (c7 * duprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)) + kdu = duprev + dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + abar75 * k5 + abar76 * k6) + + k7 = f.f1(duprev, ku, p, t + dt * c7) + u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # no b6, b7 + du = duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) # no b2, b7 + + 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 += 4 + integrator.stats.nf2 += 1 + integrator.k[1] = integrator.fsalfirst + integrator.k[2] = integrator.fsallast +end + +@muladd function perform_step!(integrator, cache::Nystrom5Cache, repeat_step = false) + @unpack t, dt, f, p = integrator + du, u = integrator.u.x + duprev, uprev = integrator.uprev.x + @unpack tmp, fsalfirst, k2, k3, k4, k5, k6, k7, k = cache + @unpack c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b3, b4, b5, b6, b7, bbar1, bbar3, bbar4, bbar5, bbar6, bbar7 = cache.tab #cache.tab? + + kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.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 + a42 * k2 + a43 * k3)) + @.. 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 ku=uprev + dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) + @.. broadcast=false kdu=duprev + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) + + f.f1(k6, kdu, ku, p, t + dt * c6) + @.. broadcast=false ku=uprev + dt * (c7 * duprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)) + @.. broadcast=false kdu=duprev + dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + abar75 * k5 + abar76 * k6) + + f.f1(k7, kdu, ku, p, t + dt * c7) + @.. broadcast=false u=uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) + @.. broadcast=false du=duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) + + f.f1(k.x[1], du, u, p, t + dt) + f.f2(k.x[2], du, u, p, t + dt) + integrator.stats.nf += 1 + integrator.stats.nf2 += 1 +end + @muladd function perform_step!(integrator, cache::Nystrom4VelocityIndependentConstantCache, repeat_step = false) @unpack t, dt, f, p = integrator diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 9962875788..2c2dd2a7c1 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -65,136 +65,136 @@ struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache end function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats}) - c1 = convert(T2, \\ ) - c2 = convert(T2, \\ ) - c3 = convert(T2, \\ ) - c4 = convert(T2, \\ ) - c5 = convert(T2, \\ ) - c6 = convert(T2, \\ ) - c7 = convert(T2, \\ ) - a21 = convert(T, \\ ) - a31 = convert(T, \\ ) - a32 = convert(T, \\ ) - a41 = convert(T, \\ ) - a42 = convert(T, \\ ) - a43 = convert(T, \\ ) - a51 = convert(T, \\ ) - a52 = convert(T, \\ ) - a53 = convert(T, \\ ) - a54 = convert(T, \\ ) - a61 = convert(T, \\ ) - a62 = convert(T, \\ ) - a63 = convert(T, \\ ) - a64 = convert(T, \\ ) - a65 = convert(T, \\ ) - a71 = convert(T, \\ ) - a72 = convert(T, \\ ) - a73 = convert(T, \\ ) - a74 = convert(T, \\ ) - a75 = convert(T, \\ ) - a76 = convert(T, \\ ) - abar21 = convert(T, \\ ) - abar31 = convert(T, \\ ) - abar32 = convert(T, \\ ) - abar41 = convert(T, \\ ) - abar42 = convert(T, \\ ) - abar43 = convert(T, \\ ) - abar51 = convert(T, \\ ) - abar52 = convert(T, \\ ) - abar53 = convert(T, \\ ) - abar54 = convert(T, \\ ) - abar61 = convert(T, \\ ) - abar62 = convert(T, \\ ) - abar63 = convert(T, \\ ) - abar64 = convert(T, \\ ) - abar65 = convert(T, \\ ) - abar71 = convert(T, \\ ) - abar72 = convert(T, \\ ) - abar73 = convert(T, \\ ) - abar74 = convert(T, \\ ) - abar75 = convert(T, \\ ) - abar76 = convert(T, \\ ) - b1 = convert(T, \\ ) - b2 = convert(T, \\ ) - b3 = convert(T, \\ ) - b4 = convert(T, \\ ) - b5 = convert(T, \\ ) - b6 = convert(T, \\ ) - b7 = convert(T, \\ ) - bbar1 = convert(T, \\ ) - bbar2 = convert(T, \\ ) - bbar3 = convert(T, \\ ) - bbar4 = convert(T, \\ ) - bbar5 = convert(T, \\ ) - bbar6 = convert(T, \\ ) - bbar7 = convert(T, \\ ) + c1 = convert(T2, 1.0) + c2 = convert(T2, 0.2051282051282051) + c3 = convert(T2, 0.3076923076923077) + c4 = convert(T2, 0.8333333333333334) + c5 = convert(T2, 0.9148936170212766) + c6 = convert(T2, 0.9999725756910926) + c7 = convert(T2, 1.0) + a21 = convert(T, 0.0210387902695595) + a31 = convert(T, 0.0236686390532544) + a32 = convert(T, 0.0236686390532544) + a41 = convert(T, 0.0337577160493827) + a42 = convert(T, 0.0) + a43 = convert(T, 0.31346450617284) + a51 = convert(T, -0.060954498687391) + a52 = convert(T, 0.145802190676171) + a53 = convert(T, 0.34307079575472) + a54 = convert(T, -0.0094033225103623) + a61 = convert(T, -0.135737843227489) + a62 = convert(T, 0.352968246663599) + a63 = convert(T, 0.268963798128547) + a64 = convert(T, 0.0138057984353428) + a65 = convert(T, 0.0) + a71 = convert(T, 0.0933527131782946) + a72 = convert(T, 0.0) + a73 = convert(T, 0.319562323318651) + a74 = convert(T, 0.138960763520679) + a75 = convert(T, -0.0518758000176242) + a76 = convert(T, 0.0) + abar21 = convert(T, 0.205128205128205) + abar31 = convert(T, 0.0769230769230769) + abar32 = convert(T, 0.230769230769231) + abar41 = convert(T, 1.06843171296296) + abar42 = convert(T, -4.09071180555556) + abar43 = convert(T, 3.85561342592593) + abar51 = convert(T, 2.44433442449694) + abar52 = convert(T, -9.53630178503099) + abar53 = convert(T, 8.17612078274186) + abar54 = convert(T, -0.169259805186522) + abar61 = convert(T, 2.9748519087319) + abar62 = convert(T, -11.3596698113208) + abar63 = convert(T, 9.45765911709872) + abar64 = convert(T, 0.162068068588807) + abar65 = convert(T, -0.234909283098676) + abar71 = convert(T, 0.0933527131782946) + abar72 = convert(T, 0.0) + abar73 = convert(T, 0.461590022571385) + abar74 = convert(T, 0.833764581124072) + abar75 = convert(T, -0.609540650207085) + abar76 = convert(T, 0.220833333333333) + b1 = convert(T, 0.0933527131782946) + b2 = convert(T, 0.0) + b3 = convert(T, 0.319562323318651) + b4 = convert(T, 0.138960763520679) + b5 = convert(T, -0.0518758000176242) + b6 = convert(T, 0.0) + b7 = convert(T, 0.0) + bbar1 = convert(T, 0.0933527131782946) + bbar2 = convert(T, 0.0) + bbar3 = convert(T, 0.461590022571385) + bbar4 = convert(T, 0.833764581124072) + bbar5 = convert(T, -0.609540650207085) + bbar6 = convert(T, 0.220833333333333) + bbar7 = convert(T, 0.0) Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) end function Nystrom5ConstantCache(T::Type, T2::Type) - c1 = convert(T2, \\ ) - c2 = convert(T2, \\ ) - c3 = convert(T2, \\ ) - c4 = convert(T2, \\ ) - c5 = convert(T2, \\ ) - c6 = convert(T2, \\ ) - c7 = convert(T2, \\ ) - a21 = convert(T, \\ ) - a31 = convert(T, \\ ) - a32 = convert(T, \\ ) - a41 = convert(T, \\ ) - a42 = convert(T, \\ ) - a43 = convert(T, \\ ) - a51 = convert(T, \\ ) - a52 = convert(T, \\ ) - a53 = convert(T, \\ ) - a54 = convert(T, \\ ) - a61 = convert(T, \\ ) - a62 = convert(T, \\ ) - a63 = convert(T, \\ ) - a64 = convert(T, \\ ) - a65 = convert(T, \\ ) - a71 = convert(T, \\ ) - a72 = convert(T, \\ ) - a73 = convert(T, \\ ) - a74 = convert(T, \\ ) - a75 = convert(T, \\ ) - a76 = convert(T, \\ ) - abar21 = convert(T, \\ ) - abar31 = convert(T, \\ ) - abar32 = convert(T, \\ ) - abar41 = convert(T, \\ ) - abar42 = convert(T, \\ ) - abar43 = convert(T, \\ ) - abar51 = convert(T, \\ ) - abar52 = convert(T, \\ ) - abar53 = convert(T, \\ ) - abar54 = convert(T, \\ ) - abar61 = convert(T, \\ ) - abar62 = convert(T, \\ ) - abar63 = convert(T, \\ ) - abar64 = convert(T, \\ ) - abar65 = convert(T, \\ ) - abar71 = convert(T, \\ ) - abar72 = convert(T, \\ ) - abar73 = convert(T, \\ ) - abar74 = convert(T, \\ ) - abar75 = convert(T, \\ ) - abar76 = convert(T, \\ ) - b1 = convert(T, \\ ) - b2 = convert(T, \\ ) - b3 = convert(T, \\ ) - b4 = convert(T, \\ ) - b5 = convert(T, \\ ) - b6 = convert(T, \\ ) - b7 = convert(T, \\ ) - bbar1 = convert(T, \\ ) - bbar2 = convert(T, \\ ) - bbar3 = convert(T, \\ ) - bbar4 = convert(T, \\ ) - bbar5 = convert(T, \\ ) - bbar6 = convert(T, \\ ) - bbar7 = convert(T, \\ ) + c1 = convert(T2, 1 // 1) + c2 = convert(T2, 8 // 39) + c3 = convert(T2, 4 // 13) + c4 = convert(T2, 5 // 6) + c5 = convert(T2, 43 // 47) + c6 = convert(T2, 36463 // 36464) + c7 = convert(T2, 1 // 1) + a21 = convert(T, 32 // 1521) + a31 = convert(T, 4 // 169) + a32 = convert(T, 4 // 169) + a41 = convert(T, 175 // 5184) + a42 = convert(T, 0 // 1) + a43 = convert(T, 1625 // 5184) + a51 = convert(T, -342497279 // 5618900760) + a52 = convert(T, 6827067 // 46824173) + a53 = convert(T, 35048741 // 102161832) + a54 = convert(T, -2201514 // 234120865) + a61 = convert(T, -7079 // 52152) + a62 = convert(T, 767 // 2173) + a63 = convert(T, 14027 // 52152) + a64 = convert(T, 30 // 2173) + a65 = convert(T, 0 // 1) + a71 = convert(T, 4817 // 51600) + a72 = convert(T, 0 // 1) + a73 = convert(T, 388869 // 1216880) + a74 = convert(T, 3276 // 23575) + a75 = convert(T, -1142053 // 22015140) + a76 = convert(T, 0 // 1) + abar21 = convert(T, 8 // 39) + abar31 = convert(T, 1 // 13) + abar32 = convert(T, 3 // 13) + abar41 = convert(T, 7385 // 6912) + abar42 = convert(T, -9425 // 2304) + abar43 = convert(T, 13325 // 3456) + abar51 = convert(T, 223324757 // 91364240) + abar52 = convert(T, -174255393 // 18272848) + abar53 = convert(T, 382840094 // 46824173) + abar54 = convert(T, -39627252 // 234120865) + abar61 = convert(T, 108475 // 36464) + abar62 = convert(T, -9633 // 848) + abar63 = convert(T, 7624604 // 806183) + abar64 = convert(T, 8100 // 49979) + abar65 = convert(T, -4568212 // 19446707) + abar71 = convert(T, 4817 // 51600) + abar72 = convert(T, 0 // 1) + abar73 = convert(T, 1685099 // 3650640) + abar74 = convert(T, 19656 // 23575) + abar75 = convert(T, -53676491 // 88060560) + abar76 = convert(T, 53 // 240) + b1 = convert(T, 4817 // 51600) + b2 = convert(T, 0 // 1) + b3 = convert(T, 388869 // 1216880) + b4 = convert(T, 3276 // 23575) + b5 = convert(T, -1142053 // 22015140) + b6 = convert(T, 0 // 1) + b7 = convert(T, 0 // 1) + bbar1 = convert(T, 4817 // 51600) + bbar2 = convert(T, 0 // 1) + bbar3 = convert(T, 1685099 // 3650640) + bbar4 = convert(T, 19656 // 23575) + bbar5 = convert(T, -53676491 // 88060560) + bbar6 = convert(T, 53 // 240) + bbar7 = convert(T, 0 // 1) Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) end diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index b4a650eec5..b367388446 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -115,6 +115,9 @@ dts = 1 .// 2 .^ (9:-1:6) 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, Nystrom5(), dense_errors = true) +@test sim.𝒪est[:l2]≈5 rtol=1e-1 +@test sim.𝒪est[:L2]≈5 rtol=1e-1 sim = test_convergence(dts, prob, Nystrom4VelocityIndependent(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -287,6 +290,9 @@ dts = 1 .// 2 .^ (9:-1:6) 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, Nystrom5(), dense_errors = true) +@test sim.𝒪est[:l2]≈4 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 sim = test_convergence(dts, prob, Nystrom4VelocityIndependent(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 From ad7dd23785eb5f2f8c0c7b913f11d1e65c8b2c3c Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 26 May 2023 10:18:03 +0200 Subject: [PATCH 04/15] Used JuliaFormatter --- src/caches/rkn_caches.jl | 7 +-- src/perform_step/rkn_perform_step.jl | 66 ++++++++++++++------- src/tableaus/rkn_tableaus.jl | 86 +++++++++++++++------------- 3 files changed, 94 insertions(+), 65 deletions(-) diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 39c4b8d630..196aeb04ae 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -36,7 +36,8 @@ 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 Nystrom5Cache{uType, rateType, reducedRateType, TabType} <: OrdinaryDiffEqMutableCache +@cache struct Nystrom5Cache{uType, rateType, reducedRateType, TabType} <: + OrdinaryDiffEqMutableCache u::uType uprev::uType fsalfirst::rateType @@ -74,11 +75,9 @@ function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - Nystrom5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) + Nystrom5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end - - @cache struct Nystrom4VelocityIndependentCache{uType, rateType, reducedRateType} <: OrdinaryDiffEqMutableCache u::uType diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index 19f2cd53f6..52a1d72adc 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -67,7 +67,7 @@ end kdu = duprev + halfdt * k₁ k₂ = f.f1(kdu, ku, p, ttmp) - ku = uprev + halfdt * duprev + eighth_dtsq * k₁ + ku = uprev + halfdt * duprev + eighth_dtsq * k₁ kdu = duprev + halfdt * k₂ k₃ = f.f1(kdu, ku, p, ttmp) @@ -105,7 +105,7 @@ end @.. broadcast=false kdu=duprev + halfdt * k₁ f.f1(k₂, kdu, ku, p, ttmp) - @.. broadcast=false ku=uprev + halfdt * duprev + eighth_dtsq * k₁ + @.. broadcast=false ku=uprev + halfdt * duprev + eighth_dtsq * k₁ @.. broadcast=false kdu=duprev + halfdt * k₂ f.f1(k₃, kdu, ku, p, ttmp) @@ -122,7 +122,7 @@ end integrator.stats.nf2 += 1 end -@muladd function perform_step!(integrator, cache::Nystrom5ConstantCache, +@muladd function perform_step!(integrator, cache::Nystrom5ConstantCache, repeat_step = false) @unpack t, dt, f, p = integrator duprev, uprev = integrator.uprev.x @@ -145,12 +145,18 @@ end kdu = duprev + dt * (abar51 * k1 + abar52 * k2 + abar53 * k3 + abar54 * k4) k5 = f.f1(duprev, ku, p, t + dt * c5) - ku = uprev + dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) - kdu = duprev + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) + ku = uprev + + dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) + kdu = duprev + + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) k6 = f.f1(duprev, ku, p, t + dt * c6) - ku = uprev + dt * (c7 * duprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)) - kdu = duprev + dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + abar75 * k5 + abar76 * k6) + ku = uprev + + dt * (c7 * duprev + + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)) + kdu = duprev + + dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + abar75 * k5 + + abar76 * k6) k7 = f.f1(duprev, ku, p, t + dt * c7) u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # no b6, b7 @@ -168,9 +174,9 @@ end @unpack t, dt, f, p = integrator du, u = integrator.u.x duprev, uprev = integrator.uprev.x - @unpack tmp, fsalfirst, k2, k3, k4, k5, k6, k7, k = cache + @unpack tmp, fsalfirst, k2, k3, k4, k5, k6, k7, k = cache @unpack c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b3, b4, b5, b6, b7, bbar1, bbar3, bbar4, bbar5, bbar6, bbar7 = cache.tab #cache.tab? - + kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2] k1 = integrator.fsalfirst.x[1] @@ -181,25 +187,41 @@ end @.. 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 + a42 * k2 + a43 * k3)) + f.f1(k3, kdu, ku, p, t + dt * c3) + @.. broadcast=false ku=uprev + + dt * (c4 * duprev + dt * (a41 * k1 + a42 * k2 + a43 * k3)) @.. 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) - + @.. 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 ku=uprev + dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) - @.. broadcast=false kdu=duprev + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) - + @.. broadcast=false ku=uprev + + dt * (c6 * duprev + + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) + @.. broadcast=false kdu=duprev + + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + + abar65 * k5) + f.f1(k6, kdu, ku, p, t + dt * c6) - @.. broadcast=false ku=uprev + dt * (c7 * duprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)) - @.. broadcast=false kdu=duprev + dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + abar75 * k5 + abar76 * k6) - + @.. broadcast=false ku=uprev + + dt * (c7 * duprev + + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + + a76 * k6)) + @.. broadcast=false kdu=duprev + + dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + + abar75 * k5 + abar76 * k6) + f.f1(k7, kdu, ku, p, t + dt * c7) - @.. broadcast=false u=uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) - @.. broadcast=false du=duprev + dt * (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) + @.. broadcast=false u=uprev + + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) + @.. broadcast=false du=duprev + + dt * + (bbar1 * k1 + bbar3 * k3 + bbar4 * k4 + bbar5 * k5 + bbar6 * k6) f.f1(k.x[1], du, u, p, t + dt) f.f2(k.x[2], du, u, p, t + dt) diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 2c2dd2a7c1..9e23199c5f 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -93,42 +93,47 @@ function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFlo a74 = convert(T, 0.138960763520679) a75 = convert(T, -0.0518758000176242) a76 = convert(T, 0.0) - abar21 = convert(T, 0.205128205128205) - abar31 = convert(T, 0.0769230769230769) - abar32 = convert(T, 0.230769230769231) - abar41 = convert(T, 1.06843171296296) - abar42 = convert(T, -4.09071180555556) - abar43 = convert(T, 3.85561342592593) - abar51 = convert(T, 2.44433442449694) - abar52 = convert(T, -9.53630178503099) - abar53 = convert(T, 8.17612078274186) - abar54 = convert(T, -0.169259805186522) - abar61 = convert(T, 2.9748519087319) - abar62 = convert(T, -11.3596698113208) - abar63 = convert(T, 9.45765911709872) - abar64 = convert(T, 0.162068068588807) - abar65 = convert(T, -0.234909283098676) - abar71 = convert(T, 0.0933527131782946) - abar72 = convert(T, 0.0) - abar73 = convert(T, 0.461590022571385) - abar74 = convert(T, 0.833764581124072) - abar75 = convert(T, -0.609540650207085) - abar76 = convert(T, 0.220833333333333) - b1 = convert(T, 0.0933527131782946) - b2 = convert(T, 0.0) - b3 = convert(T, 0.319562323318651) - b4 = convert(T, 0.138960763520679) - b5 = convert(T, -0.0518758000176242) - b6 = convert(T, 0.0) - b7 = convert(T, 0.0) - bbar1 = convert(T, 0.0933527131782946) - bbar2 = convert(T, 0.0) - bbar3 = convert(T, 0.461590022571385) - bbar4 = convert(T, 0.833764581124072) - bbar5 = convert(T, -0.609540650207085) - bbar6 = convert(T, 0.220833333333333) - bbar7 = convert(T, 0.0) - Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) + abar21 = convert(T, 0.205128205128205) + abar31 = convert(T, 0.0769230769230769) + abar32 = convert(T, 0.230769230769231) + abar41 = convert(T, 1.06843171296296) + abar42 = convert(T, -4.09071180555556) + abar43 = convert(T, 3.85561342592593) + abar51 = convert(T, 2.44433442449694) + abar52 = convert(T, -9.53630178503099) + abar53 = convert(T, 8.17612078274186) + abar54 = convert(T, -0.169259805186522) + abar61 = convert(T, 2.9748519087319) + abar62 = convert(T, -11.3596698113208) + abar63 = convert(T, 9.45765911709872) + abar64 = convert(T, 0.162068068588807) + abar65 = convert(T, -0.234909283098676) + abar71 = convert(T, 0.0933527131782946) + abar72 = convert(T, 0.0) + abar73 = convert(T, 0.461590022571385) + abar74 = convert(T, 0.833764581124072) + abar75 = convert(T, -0.609540650207085) + abar76 = convert(T, 0.220833333333333) + b1 = convert(T, 0.0933527131782946) + b2 = convert(T, 0.0) + b3 = convert(T, 0.319562323318651) + b4 = convert(T, 0.138960763520679) + b5 = convert(T, -0.0518758000176242) + b6 = convert(T, 0.0) + b7 = convert(T, 0.0) + bbar1 = convert(T, 0.0933527131782946) + bbar2 = convert(T, 0.0) + bbar3 = convert(T, 0.461590022571385) + bbar4 = convert(T, 0.833764581124072) + bbar5 = convert(T, -0.609540650207085) + bbar6 = convert(T, 0.220833333333333) + bbar7 = convert(T, 0.0) + Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, + a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, + abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, + b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) end function Nystrom5ConstantCache(T::Type, T2::Type) @@ -195,11 +200,14 @@ function Nystrom5ConstantCache(T::Type, T2::Type) bbar5 = convert(T, -53676491 // 88060560) bbar6 = convert(T, 53 // 240) bbar7 = convert(T, 0 // 1) - Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) + Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, + a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, + abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, + b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) end -################################################################################# - struct IRKN3ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache bconst1::T bconst2::T From b7747fb8f0abb4434f1ca660121376fc796362f9 Mon Sep 17 00:00:00 2001 From: Henry Date: Fri, 26 May 2023 12:01:45 +0200 Subject: [PATCH 05/15] Cleaned up the Butcher Table for Nystrom5 und changed to perform_step function so that coefficients with the value of zero are excluded. --- src/perform_step/rkn_perform_step.jl | 26 ++++----- src/tableaus/rkn_tableaus.jl | 84 ++++++++++++++-------------- 2 files changed, 55 insertions(+), 55 deletions(-) diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index 52a1d72adc..14d8cf7bdc 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -126,7 +126,7 @@ end repeat_step = false) @unpack t, dt, f, p = integrator duprev, uprev = integrator.uprev.x - @unpack c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b3, b4, b5, b6, b7, bbar1, bbar3, bbar4, bbar5, bbar6, bbar7 = cache + @unpack c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6 = cache k1 = integrator.fsalfirst.x[1] ku = uprev + (c2 * dt * duprev) + (dt^2 * a21 * k1) @@ -137,7 +137,7 @@ end kdu = duprev + dt * (abar31 * k1 + abar32 * k2) k3 = f.f1(kdu, ku, p, t + dt * c3) - ku = uprev + dt * (c4 * duprev + dt * (a41 * k1 + a42 * k2 + a43 * k3)) + 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) @@ -146,17 +146,17 @@ end k5 = f.f1(duprev, ku, p, t + dt * c5) ku = uprev + - dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) + dt * (c6 * duprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4)) # a65 = 0 kdu = duprev + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) k6 = f.f1(duprev, ku, p, t + dt * c6) ku = uprev + dt * (c7 * duprev + - dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)) + dt * (a71 * k1 + a73 * k3 + a74 * k4 + a75 * k5)) # a72 = a76 = 0 kdu = duprev + - dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + abar75 * k5 + - abar76 * k6) + dt * (abar71 * k1 + abar73 * k3 + abar74 * k4 + abar75 * k5 + + abar76 * k6) # abar72 = 0 k7 = f.f1(duprev, ku, p, t + dt * c7) u = uprev + dt * (duprev + dt * (b1 * k1 + b3 * k3 + b4 * k4 + b5 * k5)) # no b6, b7 @@ -175,7 +175,7 @@ end du, u = integrator.u.x duprev, uprev = integrator.uprev.x @unpack tmp, fsalfirst, k2, k3, k4, k5, k6, k7, k = cache - @unpack c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar72, abar73, abar74, abar75, abar76, b1, b3, b4, b5, b6, b7, bbar1, bbar3, bbar4, bbar5, bbar6, bbar7 = cache.tab #cache.tab? + @unpack c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, abar71, abar73, abar74, abar75, abar76, b1, b3, b4, b5, bbar1, bbar3, bbar4, bbar5, bbar6 = cache.tab kdu, ku = integrator.cache.tmp.x[1], integrator.cache.tmp.x[2] k1 = integrator.fsalfirst.x[1] @@ -189,7 +189,7 @@ end f.f1(k3, kdu, ku, p, t + dt * c3) @.. broadcast=false ku=uprev + - dt * (c4 * duprev + dt * (a41 * k1 + a42 * k2 + a43 * k3)) + 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) @@ -202,7 +202,7 @@ end f.f1(k5, kdu, ku, p, t + dt * c5) @.. broadcast=false ku=uprev + dt * (c6 * duprev + - dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)) + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4)) # a65 = 0 @.. broadcast=false kdu=duprev + dt * (abar61 * k1 + abar62 * k2 + abar63 * k3 + abar64 * k4 + abar65 * k5) @@ -210,11 +210,11 @@ end f.f1(k6, kdu, ku, p, t + dt * c6) @.. broadcast=false ku=uprev + dt * (c7 * duprev + - dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + - a76 * k6)) + dt * (a71 * k1 + a73 * k3 + a74 * k4 + a75 * k5)) # a72 = a76 = 0 + @.. broadcast=false kdu=duprev + - dt * (abar71 * k1 + abar72 * k2 + abar73 * k3 + abar74 * k4 + - abar75 * k5 + abar76 * k6) + dt * (abar71 * k1 + abar73 * k3 + abar74 * k4 + + abar75 * k5 + abar76 * k6) # abar72 = 0 f.f1(k7, kdu, ku, p, t + dt * c7) @.. broadcast=false u=uprev + diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 9e23199c5f..ef3caad633 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -10,7 +10,7 @@ struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache a31::T a32::T a41::T - a42::T + #a42::T a43::T a51::T a52::T @@ -20,13 +20,13 @@ struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache a62::T a63::T a64::T - a65::T + #a65::T a71::T - a72::T + #a72::T a73::T a74::T a75::T - a76::T + #a76::T abar21::T abar31::T abar32::T @@ -43,40 +43,40 @@ struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache abar64::T abar65::T abar71::T - abar72::T + #abar72::T abar73::T abar74::T abar75::T abar76::T b1::T - b2::T + #b2::T b3::T b4::T b5::T - b6::T - b7::T + #b6::T + #b7::T bbar1::T - bbar2::T + #bbar2::T bbar3::T bbar4::T bbar5::T bbar6::T - bbar7::T + #bbar7::T end function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats}) c1 = convert(T2, 1.0) c2 = convert(T2, 0.2051282051282051) - c3 = convert(T2, 0.3076923076923077) + c3 = convert(T2, 0.307692307692308) c4 = convert(T2, 0.8333333333333334) c5 = convert(T2, 0.9148936170212766) - c6 = convert(T2, 0.9999725756910926) + c6 = convert(T2, 1.0) c7 = convert(T2, 1.0) a21 = convert(T, 0.0210387902695595) a31 = convert(T, 0.0236686390532544) a32 = convert(T, 0.0236686390532544) a41 = convert(T, 0.0337577160493827) - a42 = convert(T, 0.0) + #a42 = convert(T, 0.0) a43 = convert(T, 0.31346450617284) a51 = convert(T, -0.060954498687391) a52 = convert(T, 0.145802190676171) @@ -86,13 +86,13 @@ function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFlo a62 = convert(T, 0.352968246663599) a63 = convert(T, 0.268963798128547) a64 = convert(T, 0.0138057984353428) - a65 = convert(T, 0.0) + #a65 = convert(T, 0.0) a71 = convert(T, 0.0933527131782946) - a72 = convert(T, 0.0) + #a72 = convert(T, 0.0) a73 = convert(T, 0.319562323318651) a74 = convert(T, 0.138960763520679) a75 = convert(T, -0.0518758000176242) - a76 = convert(T, 0.0) + #a76 = convert(T, 0.0) abar21 = convert(T, 0.205128205128205) abar31 = convert(T, 0.0769230769230769) abar32 = convert(T, 0.230769230769231) @@ -109,31 +109,31 @@ function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFlo abar64 = convert(T, 0.162068068588807) abar65 = convert(T, -0.234909283098676) abar71 = convert(T, 0.0933527131782946) - abar72 = convert(T, 0.0) + #abar72 = convert(T, 0.0) abar73 = convert(T, 0.461590022571385) abar74 = convert(T, 0.833764581124072) abar75 = convert(T, -0.609540650207085) abar76 = convert(T, 0.220833333333333) b1 = convert(T, 0.0933527131782946) - b2 = convert(T, 0.0) + #b2 = convert(T, 0.0) b3 = convert(T, 0.319562323318651) b4 = convert(T, 0.138960763520679) b5 = convert(T, -0.0518758000176242) - b6 = convert(T, 0.0) - b7 = convert(T, 0.0) + #b6 = convert(T, 0.0) + #b7 = convert(T, 0.0) bbar1 = convert(T, 0.0933527131782946) bbar2 = convert(T, 0.0) bbar3 = convert(T, 0.461590022571385) bbar4 = convert(T, 0.833764581124072) bbar5 = convert(T, -0.609540650207085) bbar6 = convert(T, 0.220833333333333) - bbar7 = convert(T, 0.0) - Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, - a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, - a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + #bbar7 = convert(T, 0.0) + Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, + abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, - abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, - b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) + abar71, abar73, abar74, abar75, abar76, b1, b3, b4, + b5, bbar1, bbar3, bbar4, bbar5, bbar6) end function Nystrom5ConstantCache(T::Type, T2::Type) @@ -142,13 +142,13 @@ function Nystrom5ConstantCache(T::Type, T2::Type) c3 = convert(T2, 4 // 13) c4 = convert(T2, 5 // 6) c5 = convert(T2, 43 // 47) - c6 = convert(T2, 36463 // 36464) + c6 = convert(T2, 1 // 1) # 36463 // 36464 c7 = convert(T2, 1 // 1) a21 = convert(T, 32 // 1521) a31 = convert(T, 4 // 169) a32 = convert(T, 4 // 169) a41 = convert(T, 175 // 5184) - a42 = convert(T, 0 // 1) + #a42 = convert(T, 0 // 1) a43 = convert(T, 1625 // 5184) a51 = convert(T, -342497279 // 5618900760) a52 = convert(T, 6827067 // 46824173) @@ -158,13 +158,13 @@ function Nystrom5ConstantCache(T::Type, T2::Type) a62 = convert(T, 767 // 2173) a63 = convert(T, 14027 // 52152) a64 = convert(T, 30 // 2173) - a65 = convert(T, 0 // 1) + #a65 = convert(T, 0 // 1) a71 = convert(T, 4817 // 51600) - a72 = convert(T, 0 // 1) + #a72 = convert(T, 0 // 1) a73 = convert(T, 388869 // 1216880) a74 = convert(T, 3276 // 23575) a75 = convert(T, -1142053 // 22015140) - a76 = convert(T, 0 // 1) + #a76 = convert(T, 0 // 1) abar21 = convert(T, 8 // 39) abar31 = convert(T, 1 // 13) abar32 = convert(T, 3 // 13) @@ -181,31 +181,31 @@ function Nystrom5ConstantCache(T::Type, T2::Type) abar64 = convert(T, 8100 // 49979) abar65 = convert(T, -4568212 // 19446707) abar71 = convert(T, 4817 // 51600) - abar72 = convert(T, 0 // 1) + #abar72 = convert(T, 0 // 1) abar73 = convert(T, 1685099 // 3650640) abar74 = convert(T, 19656 // 23575) abar75 = convert(T, -53676491 // 88060560) abar76 = convert(T, 53 // 240) b1 = convert(T, 4817 // 51600) - b2 = convert(T, 0 // 1) + #b2 = convert(T, 0 // 1) b3 = convert(T, 388869 // 1216880) b4 = convert(T, 3276 // 23575) b5 = convert(T, -1142053 // 22015140) - b6 = convert(T, 0 // 1) - b7 = convert(T, 0 // 1) + #b6 = convert(T, 0 // 1) + #b7 = convert(T, 0 // 1) bbar1 = convert(T, 4817 // 51600) - bbar2 = convert(T, 0 // 1) + #bbar2 = convert(T, 0 // 1) bbar3 = convert(T, 1685099 // 3650640) bbar4 = convert(T, 19656 // 23575) bbar5 = convert(T, -53676491 // 88060560) bbar6 = convert(T, 53 // 240) - bbar7 = convert(T, 0 // 1) - Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a42, a43, a51, - a52, a53, a54, a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, - a76, abar21, abar31, abar32, abar41, abar42, abar43, abar51, + #bbar7 = convert(T, 0 // 1) + Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, + a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, + abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, - abar71, abar72, abar73, abar74, abar75, abar76, b1, b2, b3, b4, - b5, b6, b7, bbar1, bbar2, bbar3, bbar4, bbar5, bbar6, bbar7) + abar71, abar73, abar74, abar75, abar76, b1, b3, b4, + b5, bbar1, bbar3, bbar4, bbar5, bbar6) end struct IRKN3ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache From c7c1a892f4cf3fb92f4dfb517de159e5de1d8f28 Mon Sep 17 00:00:00 2001 From: Henry Date: Sun, 28 May 2023 08:55:06 +0200 Subject: [PATCH 06/15] Nystrom5 is now FineRKN5. Using rational tableau only with FineRKN5 --- src/OrdinaryDiffEq.jl | 2 +- src/alg_utils.jl | 2 +- src/algorithms.jl | 4 ++-- src/caches/rkn_caches.jl | 10 +++++----- src/perform_step/rkn_perform_step.jl | 8 ++++---- test/algconvergence/partitioned_methods_tests.jl | 4 ++-- 6 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index 76dc985fda..eda6ec1b61 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -423,7 +423,7 @@ export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, export SplitEuler -export Nystrom4, Nystrom5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, +export Nystrom4, FineRKN5, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7 export ROCK2, ROCK4, RKC, IRKC, ESERK4, ESERK5, SERK2 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 4514b2acb3..a9f7a811ee 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -487,7 +487,7 @@ alg_order(alg::SofSpa10) = 10 alg_order(alg::IRKN3) = 3 alg_order(alg::Nystrom4) = 4 -alg_order(alg::Nystrom4) = 5 +alg_order(alg::FineRKN5) = 5 alg_order(alg::Nystrom4VelocityIndependent) = 4 alg_order(alg::IRKN4) = 4 alg_order(alg::Nystrom5VelocityIndependent) = 5 diff --git a/src/algorithms.jl b/src/algorithms.jl index 2457568229..c0f994e1fc 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5247,7 +5247,7 @@ E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equatio struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end """ - Nystrom5 + FineRKN5 A 5th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. @@ -5266,7 +5266,7 @@ In case the ODE Problem is not dependent on the first derivative consider using publisher={Springer} } """ -struct Nystrom5 <: OrdinaryDiffEqPartitionedAlgorithm end +struct FineRKN5 <: OrdinaryDiffEqPartitionedAlgorithm end """ Nystrom4VelocityIdependent diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 196aeb04ae..704dea9c4e 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -36,7 +36,7 @@ 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 Nystrom5Cache{uType, rateType, reducedRateType, TabType} <: +@cache struct FineRKN5Cache{uType, rateType, reducedRateType, TabType} <: OrdinaryDiffEqMutableCache u::uType uprev::uType @@ -53,7 +53,7 @@ end tab::TabType end -function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, +function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} @@ -68,14 +68,14 @@ function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, k7 = zero(reduced_rate_prototype) k = zero(rate_prototype) tmp = zero(u) - Nystrom5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp, tab) + FineRKN5Cache(u, uprev, k1, k2, k3, k4, k5, k6, k7, k, tmp, tab) end -function alg_cache(alg::Nystrom5, u, rate_prototype, ::Type{uEltypeNoUnits}, +function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} - Nystrom5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) + FineRKN5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) end @cache struct Nystrom4VelocityIndependentCache{uType, rateType, reducedRateType} <: diff --git a/src/perform_step/rkn_perform_step.jl b/src/perform_step/rkn_perform_step.jl index 14d8cf7bdc..46aa5b9a2e 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/src/perform_step/rkn_perform_step.jl @@ -4,7 +4,7 @@ ## y₁ = y₀ + hy'₀ + h²∑b̄ᵢk'ᵢ ## y'₁ = y'₀ + h∑bᵢk'ᵢ -const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, Nystrom5ConstantCache, +const NystromCCDefaultInitialization = Union{Nystrom4ConstantCache, FineRKN5ConstantCache, Nystrom4VelocityIndependentConstantCache, Nystrom5VelocityIndependentConstantCache, IRKN3ConstantCache, IRKN4ConstantCache, @@ -25,7 +25,7 @@ function initialize!(integrator, cache::NystromCCDefaultInitialization) integrator.fsalfirst = ArrayPartition((kdu, ku)) end -const NystromDefaultInitialization = Union{Nystrom4Cache, Nystrom5Cache, +const NystromDefaultInitialization = Union{Nystrom4Cache, FineRKN5Cache, Nystrom4VelocityIndependentCache, Nystrom5VelocityIndependentCache, IRKN3Cache, IRKN4Cache, @@ -122,7 +122,7 @@ end integrator.stats.nf2 += 1 end -@muladd function perform_step!(integrator, cache::Nystrom5ConstantCache, +@muladd function perform_step!(integrator, cache::FineRKN5ConstantCache, repeat_step = false) @unpack t, dt, f, p = integrator duprev, uprev = integrator.uprev.x @@ -170,7 +170,7 @@ end integrator.k[2] = integrator.fsallast end -@muladd function perform_step!(integrator, cache::Nystrom5Cache, repeat_step = false) +@muladd function perform_step!(integrator, cache::FineRKN5Cache, repeat_step = false) @unpack t, dt, f, p = integrator du, u = integrator.u.x duprev, uprev = integrator.uprev.x diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index b367388446..3b2cc2a27e 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -115,7 +115,7 @@ dts = 1 .// 2 .^ (9:-1:6) 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, Nystrom5(), dense_errors = true) +sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true) @test sim.𝒪est[:l2]≈5 rtol=1e-1 @test sim.𝒪est[:L2]≈5 rtol=1e-1 sim = test_convergence(dts, prob, Nystrom4VelocityIndependent(), dense_errors = true) @@ -290,7 +290,7 @@ dts = 1 .// 2 .^ (9:-1:6) 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, Nystrom5(), dense_errors = true) +sim = test_convergence(dts, prob, FineRKN5(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 sim = test_convergence(dts, prob, Nystrom4VelocityIndependent(), dense_errors = true) From 13669a65f53ef1c09eb24ca1b74b3acf7f59cf4e Mon Sep 17 00:00:00 2001 From: Henry Date: Sun, 28 May 2023 09:06:57 +0200 Subject: [PATCH 07/15] Changed Nystrom5 to FineRKN5 --- src/tableaus/rkn_tableaus.jl | 77 ++---------------------------------- 1 file changed, 3 insertions(+), 74 deletions(-) diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index ef3caad633..3eb2c8339a 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -1,4 +1,4 @@ -struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache +struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache c1::T2 c2::T2 c3::T2 @@ -64,79 +64,8 @@ struct Nystrom5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache #bbar7::T end -function Nystrom5ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloats}) - c1 = convert(T2, 1.0) - c2 = convert(T2, 0.2051282051282051) - c3 = convert(T2, 0.307692307692308) - c4 = convert(T2, 0.8333333333333334) - c5 = convert(T2, 0.9148936170212766) - c6 = convert(T2, 1.0) - c7 = convert(T2, 1.0) - a21 = convert(T, 0.0210387902695595) - a31 = convert(T, 0.0236686390532544) - a32 = convert(T, 0.0236686390532544) - a41 = convert(T, 0.0337577160493827) - #a42 = convert(T, 0.0) - a43 = convert(T, 0.31346450617284) - a51 = convert(T, -0.060954498687391) - a52 = convert(T, 0.145802190676171) - a53 = convert(T, 0.34307079575472) - a54 = convert(T, -0.0094033225103623) - a61 = convert(T, -0.135737843227489) - a62 = convert(T, 0.352968246663599) - a63 = convert(T, 0.268963798128547) - a64 = convert(T, 0.0138057984353428) - #a65 = convert(T, 0.0) - a71 = convert(T, 0.0933527131782946) - #a72 = convert(T, 0.0) - a73 = convert(T, 0.319562323318651) - a74 = convert(T, 0.138960763520679) - a75 = convert(T, -0.0518758000176242) - #a76 = convert(T, 0.0) - abar21 = convert(T, 0.205128205128205) - abar31 = convert(T, 0.0769230769230769) - abar32 = convert(T, 0.230769230769231) - abar41 = convert(T, 1.06843171296296) - abar42 = convert(T, -4.09071180555556) - abar43 = convert(T, 3.85561342592593) - abar51 = convert(T, 2.44433442449694) - abar52 = convert(T, -9.53630178503099) - abar53 = convert(T, 8.17612078274186) - abar54 = convert(T, -0.169259805186522) - abar61 = convert(T, 2.9748519087319) - abar62 = convert(T, -11.3596698113208) - abar63 = convert(T, 9.45765911709872) - abar64 = convert(T, 0.162068068588807) - abar65 = convert(T, -0.234909283098676) - abar71 = convert(T, 0.0933527131782946) - #abar72 = convert(T, 0.0) - abar73 = convert(T, 0.461590022571385) - abar74 = convert(T, 0.833764581124072) - abar75 = convert(T, -0.609540650207085) - abar76 = convert(T, 0.220833333333333) - b1 = convert(T, 0.0933527131782946) - #b2 = convert(T, 0.0) - b3 = convert(T, 0.319562323318651) - b4 = convert(T, 0.138960763520679) - b5 = convert(T, -0.0518758000176242) - #b6 = convert(T, 0.0) - #b7 = convert(T, 0.0) - bbar1 = convert(T, 0.0933527131782946) - bbar2 = convert(T, 0.0) - bbar3 = convert(T, 0.461590022571385) - bbar4 = convert(T, 0.833764581124072) - bbar5 = convert(T, -0.609540650207085) - bbar6 = convert(T, 0.220833333333333) - #bbar7 = convert(T, 0.0) - Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, - a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, - abar21, abar31, abar32, abar41, abar42, abar43, abar51, - abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, - abar71, abar73, abar74, abar75, abar76, b1, b3, b4, - b5, bbar1, bbar3, bbar4, bbar5, bbar6) -end -function Nystrom5ConstantCache(T::Type, T2::Type) +function FineRKN5ConstantCache(T::Type, T2::Type) c1 = convert(T2, 1 // 1) c2 = convert(T2, 8 // 39) c3 = convert(T2, 4 // 13) @@ -200,7 +129,7 @@ function Nystrom5ConstantCache(T::Type, T2::Type) bbar5 = convert(T, -53676491 // 88060560) bbar6 = convert(T, 53 // 240) #bbar7 = convert(T, 0 // 1) - Nystrom5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, + FineRKN5ConstantCache(c1, c2, c3, c4, c5, c6, c7, a21, a31, a32, a41, a43, a51, a52, a53, a54, a61, a62, a63, a64, a71, a73, a74, a75, abar21, abar31, abar32, abar41, abar42, abar43, abar51, abar52, abar53, abar54, abar61, abar62, abar63, abar64, abar65, From c06063ecc31605ee9e79dccbdfe847a632ae64b7 Mon Sep 17 00:00:00 2001 From: Henry Date: Tue, 30 May 2023 13:34:18 +0200 Subject: [PATCH 08/15] Added convergence tests for --- src/caches/rkn_caches.jl | 2 +- test/algconvergence/partitioned_methods_tests.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 704dea9c4e..3f336efebf 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -58,7 +58,7 @@ function alg_cache(alg::FineRKN5, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} reduced_rate_prototype = rate_prototype.x[2] - tab = Nystrom5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) + tab = FineRKN5ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) k1 = zero(rate_prototype) k2 = zero(reduced_rate_prototype) k3 = zero(reduced_rate_prototype) diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 3b2cc2a27e..8273f3791f 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -115,9 +115,6 @@ dts = 1 .// 2 .^ (9:-1:6) 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, FineRKN5(), dense_errors = true) -@test sim.𝒪est[:l2]≈5 rtol=1e-1 -@test sim.𝒪est[:L2]≈5 rtol=1e-1 sim = test_convergence(dts, prob, Nystrom4VelocityIndependent(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -131,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, FineRKN5(), dense_errors = true) +@test sim.𝒪est[:l2]≈6 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 dts = 1.0 ./ 2.0 .^ (2:-1:-2) sim = test_convergence(dts, prob, SofSpa10(), dense_errors = true) From b187d0ded8815c410589fd38a2b613cbe4c7c206 Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Tue, 30 May 2023 13:47:58 +0000 Subject: [PATCH 09/15] Added missing out of olace convergence test for ```FineRKN5```. --- src/caches/rkn_caches.jl | 1 - src/tableaus/rkn_tableaus.jl | 1 - test/algconvergence/partitioned_methods_tests.jl | 6 +++--- 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/caches/rkn_caches.jl b/src/caches/rkn_caches.jl index 3f336efebf..8780fef10a 100644 --- a/src/caches/rkn_caches.jl +++ b/src/caches/rkn_caches.jl @@ -41,7 +41,6 @@ end u::uType uprev::uType fsalfirst::rateType - #k1::rateType k2::reducedRateType k3::reducedRateType k4::reducedRateType diff --git a/src/tableaus/rkn_tableaus.jl b/src/tableaus/rkn_tableaus.jl index 3eb2c8339a..3ded6acaee 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/src/tableaus/rkn_tableaus.jl @@ -64,7 +64,6 @@ struct FineRKN5ConstantCache{T, T2} <: OrdinaryDiffEqConstantCache #bbar7::T end - function FineRKN5ConstantCache(T::Type, T2::Type) c1 = convert(T2, 1 // 1) c2 = convert(T2, 8 // 39) diff --git a/test/algconvergence/partitioned_methods_tests.jl b/test/algconvergence/partitioned_methods_tests.jl index 8273f3791f..2b99f25fca 100644 --- a/test/algconvergence/partitioned_methods_tests.jl +++ b/test/algconvergence/partitioned_methods_tests.jl @@ -290,9 +290,6 @@ dts = 1 .// 2 .^ (9:-1:6) 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, FineRKN5(), dense_errors = true) -@test sim.𝒪est[:l2]≈4 rtol=1e-1 -@test sim.𝒪est[:L2]≈4 rtol=1e-1 sim = test_convergence(dts, prob, Nystrom4VelocityIndependent(), dense_errors = true) @test sim.𝒪est[:l2]≈4 rtol=1e-1 @test sim.𝒪est[:L2]≈4 rtol=1e-1 @@ -306,6 +303,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, FineRKN5(), dense_errors = true) +@test sim.𝒪est[:l2]≈6 rtol=1e-1 +@test sim.𝒪est[:L2]≈4 rtol=1e-1 dts = 1.0 ./ 2.0 .^ (2:-1:-2) sim = test_convergence(dts, prob, SofSpa10(), dense_errors = true) From d132868a7b303fae881c1f2be49d03976ac8296d Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 1 Jun 2023 12:44:34 +0200 Subject: [PATCH 10/15] Update src/algorithms.jl Co-authored-by: Hendrik Ranocha --- src/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index c0f994e1fc..c022b82859 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5247,7 +5247,7 @@ E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equatio struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end """ - FineRKN5 + FineRKN5() A 5th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. From 22b2b7cb91e553a37889a2318dc9df210431833a Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 1 Jun 2023 12:44:48 +0200 Subject: [PATCH 11/15] Update src/algorithms.jl Co-authored-by: Hendrik Ranocha --- src/algorithms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index c022b82859..528e2ab2a9 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5251,7 +5251,7 @@ struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end A 5th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. -In case the ODE Problem is not dependent on the first derivative consider using +In case the ODE does not depend on the first derivative, consider using [`Nystrom5VelocityIndependent`](@ref) to increase performance. ## References From 4a27a6ee20f322fec7ba94bf7b8ded5f1d5bfd6c Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 1 Jun 2023 12:45:04 +0200 Subject: [PATCH 12/15] Update src/algorithms.jl Co-authored-by: Hendrik Ranocha --- src/algorithms.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 528e2ab2a9..490c78928f 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5249,7 +5249,8 @@ struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end """ FineRKN5() -A 5th order explicit Runge-Kutta-Nyström method which can be applied directly on second order ODEs. Can only be used with fixed time steps. +A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. +Can only be used with fixed time steps. In case the ODE does not depend on the first derivative, consider using [`Nystrom5VelocityIndependent`](@ref) to increase performance. From 6e4a43135e6afbd62faa59b93e10f1f410f7147e Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 1 Jun 2023 12:45:20 +0200 Subject: [PATCH 13/15] Update src/algorithms.jl Co-authored-by: Hendrik Ranocha --- src/algorithms.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index 490c78928f..8a827bd316 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5256,8 +5256,9 @@ In case the ODE does not depend on the first derivative, consider using [`Nystrom5VelocityIndependent`](@ref) to increase performance. ## References +``` @article{fine1987low, - title={Low order practical Runge-Kutta-Nystr{\"o}m methods}, + title={Low order practical {R}unge-{K}utta-{N}ystr{\"o}m methods}, author={Fine, Jerry Michael}, journal={Computing}, volume={38}, @@ -5266,6 +5267,7 @@ In case the ODE does not depend on the first derivative, consider using year={1987}, publisher={Springer} } +``` """ struct FineRKN5 <: OrdinaryDiffEqPartitionedAlgorithm end From 23c8bf4ffb596db2d8b46e3b96d89e798f6aa6a5 Mon Sep 17 00:00:00 2001 From: Henry Langner <128313572+HenryLangner@users.noreply.github.com> Date: Thu, 1 Jun 2023 10:47:45 +0000 Subject: [PATCH 14/15] Added ```FineRKN5``` to ```docs/src/dynamical/nystrom.m```. --- docs/src/dynamical/nystrom.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/dynamical/nystrom.md b/docs/src/dynamical/nystrom.md index 2c63ece155..288d5c92c9 100644 --- a/docs/src/dynamical/nystrom.md +++ b/docs/src/dynamical/nystrom.md @@ -5,6 +5,7 @@ IRKN3 Nystrom4 Nystrom4VelocityIndependent IRKN4 +FineRKN5 Nystrom5VelocityIndependent DPRKN6 DPRKN8 From 419d07a285ce0ad1ac9a49e1c8d8528f1e6dd372 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 7 Jun 2023 14:02:57 +0200 Subject: [PATCH 15/15] Update src/algorithms.jl --- src/algorithms.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index effe835050..2143adf150 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -5257,8 +5257,8 @@ struct Nystrom4 <: OrdinaryDiffEqPartitionedAlgorithm end A 5th order explicit Runge-Kutta-Nyström method which can be applied directly to second order ODEs. Can only be used with fixed time steps. -In case the ODE does not depend on the first derivative, consider using -[`Nystrom5VelocityIndependent`](@ref) to increase performance. +This method requires that the acceleration equation needs to be independent of the velocity +otherwise it results in convergence order loss. ## References ```