diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index bc7b5a436d..97989c7464 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -22,6 +22,7 @@ jobs: - AlgConvergence_II - AlgConvergence_III - Downstream + - Symplectic - Extrapolation - StabilizedRK - StabilizedIRK diff --git a/lib/OrdinaryDiffEqRKN/Project.toml b/lib/OrdinaryDiffEqRKN/Project.toml new file mode 100644 index 0000000000..05d567ea4f --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/Project.toml @@ -0,0 +1,23 @@ +name = "OrdinaryDiffEqRKN" +uuid = "af6ede74-add8-4cfd-b1df-9a4dbb109d7a" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898" +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" + +[compat] +julia = "1.10" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl new file mode 100644 index 0000000000..7ddfbc2c8e --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl @@ -0,0 +1,30 @@ +module OrdinaryDiffEqRKN + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, alg_extrapolates, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqAdaptivePartitionedAlgorithm, + OrdinaryDiffEqPartitionedAlgorithm, + OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, _ode_interpolant, + trivial_limiter!, _ode_interpolant!, _ode_addsteps! +using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools +using DiffEqBase: @def, @tight_loop_macros + +include("algorithms.jl") +include("alg_utils.jl") +include("rkn_tableaus.jl") +include("rkn_caches.jl") +include("interp_func.jl") +include("interpolants.jl") +include("rkn_perform_step.jl") + +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, + Nystrom5VelocityIndependent, + IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, + RKN4 + +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/alg_utils.jl b/lib/OrdinaryDiffEqRKN/src/alg_utils.jl new file mode 100644 index 0000000000..8f22b97e9c --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/alg_utils.jl @@ -0,0 +1,20 @@ +alg_extrapolates(alg::IRKN3) = true +alg_extrapolates(alg::IRKN4) = true + +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 +alg_order(alg::Nystrom5VelocityIndependent) = 5 +alg_order(alg::DPRKN4) = 4 +alg_order(alg::DPRKN5) = 5 +alg_order(alg::DPRKN6) = 6 +alg_order(alg::DPRKN6FM) = 6 +alg_order(alg::DPRKN8) = 8 +alg_order(alg::DPRKN12) = 12 +alg_order(alg::ERKN4) = 4 +alg_order(alg::ERKN5) = 5 +alg_order(alg::ERKN7) = 7 +alg_order(alg::RKN4) = 4 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/algorithms.jl b/lib/OrdinaryDiffEqRKN/src/algorithms.jl new file mode 100644 index 0000000000..f1cf14fd93 --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/algorithms.jl @@ -0,0 +1,341 @@ +""" + IRKN3 + +Improved Runge-Kutta-Nyström method of order three, which minimizes the amount of evaluated functions in each step. Fixed time steps only. + +Second order ODE should not depend on the first derivative. + +## References + +@article{rabiei2012numerical, +title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, +author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, +publisher={Citeseer} +} +""" +struct IRKN3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + Nystrom4 + +A 4th 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 +[`Nystrom4VelocityIndependent`](@ref) to increase performance. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +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() + +A 5th 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 FineRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + Nystrom4VelocityIdependent + +A 4th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). + +More efficient then [`Nystrom4`](@ref) on velocity independent problems, since less evaluations are needed. + +Fixed time steps only. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom4VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + IRKN4 + +Improves Runge-Kutta-Nyström method of order four, which minimizes the amount of evaluated functions in each step. Fixed time steps only. + +Second order ODE should not be dependent on the first derivative. + +Recommended for smooth problems with expensive functions to evaluate. + +## References + +@article{rabiei2012numerical, +title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, +author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, +publisher={Citeseer} +} +""" +struct IRKN4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + Nystrom5VelocityIndependent + +A 5th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). +Fixed time steps only. + +## References + +E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. +Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, +Springer-Verlag. +""" +struct Nystrom5VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end + +""" + DPRKN4 + +4th order explicit Runge-Kutta-Nyström methods. The second order ODE should not depend on the first derivative. + +## References + +@article{Dormand1987FamiliesOR, +title={Families of Runge-Kutta-Nystrom Formulae}, +author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, +journal={Ima Journal of Numerical Analysis}, +year={1987}, +volume={7}, +pages={235-250} +} +""" +struct DPRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN5 + +5th order explicit Runge-Kutta-Nyström mehod. The second order ODE should not depend on the first derivative. + +## References + +@article{Bettis1973ARN, +title={A Runge-Kutta Nystrom algorithm}, +author={Dale G. Bettis}, +journal={Celestial mechanics}, +year={1973}, +volume={8}, +pages={229-233}, +publisher={Springer} +} +""" +struct DPRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN6 + +6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. Free 6th order interpolant. + +## References + +@article{dormand1987runge, +title={Runge-kutta-nystrom triples}, +author={Dormand, JR and Prince, PJ}, +journal={Computers \\& Mathematics with Applications}, +volume={13}, +number={12}, +pages={937--949}, +year={1987}, +publisher={Elsevier} +} +""" +struct DPRKN6 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN6FM + +6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Compared to [`DPRKN6`](@ref), this method has smaller truncation error coefficients which leads to performance gain +when only the main solution points are considered. + +## References + +@article{Dormand1987FamiliesOR, +title={Families of Runge-Kutta-Nystrom Formulae}, +author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, +journal={Ima Journal of Numerical Analysis}, +year={1987}, +volume={7}, +pages={235-250} +} +""" +struct DPRKN6FM <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN8 + +8th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Not as efficient as [`DPRKN12`](@ref) when high accuracy is needed, however this solver is competitive with +[`DPRKN6`](@ref) at lax tolerances and, depending on the problem, might be a good option between performance and accuracy. + +## References + +@article{dormand1987high, +title={High-order embedded Runge-Kutta-Nystrom formulae}, +author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, +journal={IMA Journal of Numerical Analysis}, +volume={7}, +number={4}, +pages={423--430}, +year={1987}, +publisher={Oxford University Press} +} +""" +struct DPRKN8 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + DPRKN12 + +12th order explicit Rugne-Kutta-Nyström method. The second order ODE should not depend on the first derivative. + +Most efficient when high accuracy is needed. + +## References + +@article{dormand1987high, +title={High-order embedded Runge-Kutta-Nystrom formulae}, +author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, +journal={IMA Journal of Numerical Analysis}, +volume={7}, +number={4}, +pages={423--430}, +year={1987}, +publisher={Oxford University Press} +} +""" +struct DPRKN12 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN4 + +Embedded 4(3) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic problems. + +## References + +@article{demba2017embedded, +title={An Embedded 4 (3) Pair of Explicit Trigonometrically-Fitted Runge-Kutta-Nystr{\"o}m Method for Solving Periodic Initial Value Problems}, +author={Demba, MA and Senu, N and Ismail, F}, +journal={Applied Mathematical Sciences}, +volume={11}, +number={17}, +pages={819--838}, +year={2017} +} +""" +struct ERKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN5 + +Embedded 5(4) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic problems. + +## References + +@article{demba20165, +title={A 5 (4) Embedded Pair of Explicit Trigonometrically-Fitted Runge--Kutta--Nystr{\"o}m Methods for the Numerical Solution of Oscillatory Initial Value Problems}, +author={Demba, Musa A and Senu, Norazak and Ismail, Fudziah}, +journal={Mathematical and Computational Applications}, +volume={21}, +number={4}, +pages={46}, +year={2016}, +publisher={Multidisciplinary Digital Publishing Institute} +} +""" +struct ERKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" + ERKN7 + +Embedded pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. + +The second order ODE should not depend on the first derivative. + +Uses adaptive step size control. This method is extra efficient on periodic Problems. + +## References + +@article{SimosOnHO, +title={On high order Runge-Kutta-Nystr{\"o}m pairs}, +author={Theodore E. Simos and Ch. Tsitouras}, +journal={J. Comput. Appl. Math.}, +volume={400}, +pages={113753} +} +""" +struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end + +""" +3 stage fourth order Runge-Kutta Nystrom method to solve second order linear inhomogeneous IVPs. + +Does not include an adaptive method. Solves for for d-dimensional differential systems of second order linear inhomogeneous equations. + +!!! warn + + This method is only fourth order for these systems, the method is second order otherwise! + +## References + +@article{MONTIJANO2024115533, +title = {Explicit Runge–Kutta–Nyström methods for the numerical solution of second order linear inhomogeneous IVPs}, +author = {J.I. Montijano and L. Rández and M. Calvo}, +journal = {Journal of Computational and Applied Mathematics}, +volume = {438}, +pages = {115533}, +year = {2024}, +} +""" +struct RKN4 <: OrdinaryDiffEqAlgorithm end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/interp_func.jl b/lib/OrdinaryDiffEqRKN/src/interp_func.jl new file mode 100644 index 0000000000..f5b2c03a2b --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/interp_func.jl @@ -0,0 +1,7 @@ +function DiffEqBase.interp_summary(::Type{cacheType}, + dense::Bool) where { + cacheType <: + Union{DPRKN6ConstantCache, + DPRKN6Cache}} +dense ? "specialized 6th order interpolation" : "1st order linear" +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/src/interpolants.jl b/lib/OrdinaryDiffEqRKN/src/interpolants.jl new file mode 100644 index 0000000000..4e919f9f47 --- /dev/null +++ b/lib/OrdinaryDiffEqRKN/src/interpolants.jl @@ -0,0 +1,141 @@ +@def dprkn6unpack begin + if cache isa OrdinaryDiffEqMutableCache + @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache.tab + else + @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache + end +end + +@def dprkn6pre0 begin + @dprkn6unpack + b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) + b3Θ = Θ * @evalpoly(Θ, r31, r32, r33, r34) + b4Θ = Θ * @evalpoly(Θ, r41, r42, r43, r44) + b5Θ = Θ * @evalpoly(Θ, r51, r52, r53, r54) + b6Θ = Θ * @evalpoly(Θ, r61, r62, r63, r64) + + bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14) + bp3Θ = Θ * @evalpoly(Θ, rp31, rp32, rp33, rp34) + bp4Θ = Θ * @evalpoly(Θ, rp41, rp42, rp43, rp44) + bp5Θ = Θ * @evalpoly(Θ, rp51, rp52, rp53, rp54) + bp6Θ = Θ * @evalpoly(Θ, rp61, rp62, rp63, rp64) + + kk1, kk2, kk3 = k + k1, k2 = kk1.x + k3, k4 = kk2.x + k5, k6 = kk3.x + + duprev, uprev = y₀.x + dtsq = dt^2 +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +return ArrayPartition( + duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6), + uprev + + dt * Θ * + (duprev + + dt * Θ * (b1Θ * k1 + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +return ArrayPartition( + duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]), + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs]))) +end + +@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Number, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +halfsize = length(y₀) ÷ 2 +if idxs <= halfsize + duprev[idxs] + + dt * Θ * + (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + + bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]) +else + idxs = idxs - halfsize + uprev[idxs] + + dt * Θ * + (duprev[idxs] + + dt * Θ * + (b1Θ * k1[idxs] + + b3Θ * k3[idxs] + + b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs])) +end +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, + idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +@inbounds @.. broadcast=false out.x[2]=uprev + + dt * Θ * + (duprev + + dt * Θ * + (b1Θ * k1 + + b3Θ * k3 + + b4Θ * k4 + b5Θ * k5 + b6Θ * k6)) +@inbounds @.. broadcast=false out.x[1]=duprev + + dt * Θ * + (bp1Θ * k1 + bp3Θ * k3 + + bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6) +#for i in eachindex(out.x[1]) +# out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + +# b3Θ*k3[i] + +# b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) +# out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + +# bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) +#end +out +end + +@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, + cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, + T::Type{Val{0}}, differential_vars::Nothing) +@dprkn6pre0 +halfsize = length(y₀) ÷ 2 +isfirsthalf = idxs .<= halfsize +secondhalf = idxs .> halfsize +firstidxs = idxs[isfirsthalf] +secondidxs_shifted = idxs[secondhalf] +secondidxs = secondidxs_shifted .- halfsize + +@views @.. broadcast=false out[secondhalf]=uprev[secondidxs] + + dt * Θ * + (duprev[secondidxs] + + dt * Θ * + (b1Θ * k1[secondidxs] + + b3Θ * k3[secondidxs] + + b4Θ * k4[secondidxs] + + b5Θ * k5[secondidxs] + + b6Θ * k6[secondidxs])) +@views @.. broadcast=false out[isfirsthalf]=duprev[firstidxs] + + dt * Θ * + (bp1Θ * k1[firstidxs] + + bp3Θ * k3[firstidxs] + + bp4Θ * k4[firstidxs] + + bp5Θ * k5[firstidxs] + + bp6Θ * k6[firstidxs]) +out +end \ No newline at end of file diff --git a/src/caches/rkn_caches.jl b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl similarity index 99% rename from src/caches/rkn_caches.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_caches.jl index 02014deef1..9119c135f3 100644 --- a/src/caches/rkn_caches.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_caches.jl @@ -712,4 +712,4 @@ function alg_cache(alg::RKN4, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} RKN4ConstantCache() -end +end \ No newline at end of file diff --git a/src/perform_step/rkn_perform_step.jl b/lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl similarity index 99% rename from src/perform_step/rkn_perform_step.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl index 70077f9268..1bacc80b93 100644 --- a/src/perform_step/rkn_perform_step.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_perform_step.jl @@ -1907,4 +1907,4 @@ end integrator.stats.nf += 2 integrator.stats.nf2 += 1 -end +end \ No newline at end of file diff --git a/src/tableaus/rkn_tableaus.jl b/lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl similarity index 99% rename from src/tableaus/rkn_tableaus.jl rename to lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl index c9b50e05ff..63dc672144 100644 --- a/src/tableaus/rkn_tableaus.jl +++ b/lib/OrdinaryDiffEqRKN/src/rkn_tableaus.jl @@ -2469,4 +2469,4 @@ function DPRKN12ConstantCache(T::Type{<:CompiledFloats}, T2::Type{<:CompiledFloa convert(T, 0.015662325288859434), convert(T, -0.029155189630077964), convert(T, 0.025)) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqRKN/test/runtests.jl b/lib/OrdinaryDiffEqRKN/test/runtests.jl new file mode 100644 index 0000000000..e69de29bb2 diff --git a/lib/OrdinaryDiffEqSymplecticRK/Project.toml b/lib/OrdinaryDiffEqSymplecticRK/Project.toml new file mode 100644 index 0000000000..d25f5cfd73 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/Project.toml @@ -0,0 +1,19 @@ +name = "OrdinaryDiffEqSymplecticRK" +uuid = "fa646aed-7ef9-47eb-84c4-9443fc8cbfa8" +authors = ["ParamThakkar123 "] +version = "1.0.0" + +[deps] +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" + +[extras] +DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +julia = "1.10" + +[targets] +test = ["DiffEqDevTools", "Random", "SafeTestsets", "Test"] \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl new file mode 100644 index 0000000000..7fc1144e47 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl @@ -0,0 +1,26 @@ +module OrdinaryDiffEqSymplecticRK + +import OrdinaryDiffEq: alg_order, calculate_residuals!, + initialize!, perform_step!, @unpack, unwrap_alg, + calculate_residuals, + OrdinaryDiffEqAlgorithm, + OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache, + OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEqPartitionedAlgorithm, CompiledFloats, uses_uprev, + alg_cache, _vec, _reshape, @cache, isfsal, full_cache, + constvalue, _unwrap_val, du_alias_or_new, + explicit_rk_docstring, trivial_limiter!, + _ode_interpolant!, _ode_addsteps! +using DiffEqBase, FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools + + +include("algorithms.jl") +include("alg_utils.jl") +include("symplectic_caches.jl") +include("symplectic_tableaus.jl") +include("symplectic_perform_step.jl") + +export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, + McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, + CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 + +end diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl b/lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl new file mode 100644 index 0000000000..ee6e0bee37 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/src/alg_utils.jl @@ -0,0 +1,17 @@ +alg_order(alg::SymplecticEuler) = 1 +alg_order(alg::VelocityVerlet) = 2 +alg_order(alg::VerletLeapfrog) = 2 +alg_order(alg::PseudoVerletLeapfrog) = 2 +alg_order(alg::McAte2) = 2 +alg_order(alg::Ruth3) = 3 +alg_order(alg::McAte3) = 3 +alg_order(alg::McAte4) = 4 +alg_order(alg::CandyRoz4) = 4 +alg_order(alg::CalvoSanz4) = 4 +alg_order(alg::McAte42) = 4 +alg_order(alg::McAte5) = 5 +alg_order(alg::Yoshida6) = 6 +alg_order(alg::KahanLi6) = 6 +alg_order(alg::McAte8) = 8 +alg_order(alg::KahanLi8) = 8 +alg_order(alg::SofSpa10) = 10 \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl b/lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl new file mode 100644 index 0000000000..22f0f9258b --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/src/algorithms.jl @@ -0,0 +1,224 @@ +struct SymplecticEuler <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{verlet1967computer, +title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, +author={Verlet, Loup}, +journal={Physical review}, +volume={159}, +number={1}, +pages={98}, +year={1967}, +publisher={APS} +} +""" + +struct VelocityVerlet <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{verlet1967computer, +title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, +author={Verlet, Loup}, +journal={Physical review}, +volume={159}, +number={1}, +pages={98}, +year={1967}, +publisher={APS} +} +""" + +struct VerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{verlet1967computer, +title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, +author={Verlet, Loup}, +journal={Physical review}, +volume={159}, +number={1}, +pages={98}, +year={1967}, +publisher={APS} +} +""" + +struct PseudoVerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte2 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{ruth1983canonical, +title={A canonical integration technique}, +author={Ruth, Ronald D}, +journal={IEEE Trans. Nucl. Sci.}, +volume={30}, +number={CERN-LEP-TH-83-14}, +pages={2669--2671}, +year={1983} +} +""" + +struct Ruth3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte3 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{candy1991symplectic, +title={A symplectic integration algorithm for separable Hamiltonian functions}, +author={Candy, J and Rozmus, W}, +journal={Journal of Computational Physics}, +volume={92}, +number={1}, +pages={230--256}, +year={1991}, +publisher={Elsevier} +} +""" + +struct CandyRoz4 <: OrdinaryDiffEqPartitionedAlgorithm end + +struct McAte4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{sanz1993symplectic, +title={Symplectic numerical methods for Hamiltonian problems}, +author={Sanz-Serna, Jes{\'u}s Maria and Calvo, Mari-Paz}, +journal={International Journal of Modern Physics C}, +volume={4}, +number={02}, +pages={385--392}, +year={1993}, +publisher={World Scientific} +} +""" + +struct CalvoSanz4 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte42 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1992accuracy, +title={The accuracy of symplectic integrators}, +author={McLachlan, Robert I and Atela, Pau}, +journal={Nonlinearity}, +volume={5}, +number={2}, +pages={541}, +year={1992}, +publisher={IOP Publishing} +} +""" + +struct McAte5 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{yoshida1990construction, +title={Construction of higher order symplectic integrators}, +author={Yoshida, Haruo}, +journal={Physics letters A}, +volume={150}, +number={5-7}, +pages={262--268}, +year={1990}, +publisher={Elsevier} +} +""" + +struct Yoshida6 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{kahan1997composition, +title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, +author={Kahan, William and Li, Ren-Cang}, +journal={Mathematics of computation}, +volume={66}, +number={219}, +pages={1089--1099}, +year={1997} +} +""" + +struct KahanLi6 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{mclachlan1995numerical, +title={On the numerical integration of ordinary differential equations by symmetric composition methods}, +author={McLachlan, Robert I}, +journal={SIAM Journal on Scientific Computing}, +volume={16}, +number={1}, +pages={151--168}, +year={1995}, +publisher={SIAM} +} +""" + +struct McAte8 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{kahan1997composition, +title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, +author={Kahan, William and Li, Ren-Cang}, +journal={Mathematics of computation}, +volume={66}, +number={219}, +pages={1089--1099}, +year={1997} +} +""" + +struct KahanLi8 <: OrdinaryDiffEqPartitionedAlgorithm end + +""" +@article{sofroniou2005derivation, +title={Derivation of symmetric composition constants for symmetric integrators}, +author={Sofroniou, Mark and Spaletta, Giulia}, +journal={Optimization Methods and Software}, +volume={20}, +number={4-5}, +pages={597--613}, +year={2005}, +publisher={Taylor \\& Francis} +} +""" +struct SofSpa10 <: OrdinaryDiffEqPartitionedAlgorithm end \ No newline at end of file diff --git a/src/caches/symplectic_caches.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl similarity index 99% rename from src/caches/symplectic_caches.jl rename to lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl index 71c55340dc..483fa6996b 100644 --- a/src/caches/symplectic_caches.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_caches.jl @@ -416,4 +416,4 @@ function alg_cache(alg::SofSpa10, u, rate_prototype, ::Type{uEltypeNoUnits}, dt, reltol, p, calck, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} SofSpa10ConstantCache(constvalue(uBottomEltypeNoUnits), constvalue(tTypeNoUnits)) -end +end \ No newline at end of file diff --git a/src/perform_step/symplectic_perform_step.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl similarity index 99% rename from src/perform_step/symplectic_perform_step.jl rename to lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl index 413f1882a2..f309e50098 100644 --- a/src/perform_step/symplectic_perform_step.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_perform_step.jl @@ -2022,4 +2022,4 @@ end integrator.stats.nf += 36 integrator.stats.nf2 += 36 store_symp_state!(integrator, cache, kdu, ku) -end +end \ No newline at end of file diff --git a/src/tableaus/symplectic_tableaus.jl b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_tableaus.jl similarity index 99% rename from src/tableaus/symplectic_tableaus.jl rename to lib/OrdinaryDiffEqSymplecticRK/src/symplectic_tableaus.jl index 98676979c3..c7bee07b0c 100644 --- a/src/tableaus/symplectic_tableaus.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/src/symplectic_tableaus.jl @@ -804,4 +804,4 @@ function SofSpa10ConstantCache(T::Type, T2::Type) b19, b20, b21, b22, b23, b24, b25, b26, b27, b28, b29, b30, b31, b32, b33, b34, b35, b36) -end +end \ No newline at end of file diff --git a/lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl b/lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl new file mode 100644 index 0000000000..231d0ee343 --- /dev/null +++ b/lib/OrdinaryDiffEqSymplecticRK/test/runtests.jl @@ -0,0 +1,3 @@ +using SafeTestsets + +@time @safetestset "Synplectic Tests" include("symplectic_tests.jl") \ No newline at end of file diff --git a/test/algconvergence/symplectic_tests.jl b/lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl similarity index 99% rename from test/algconvergence/symplectic_tests.jl rename to lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl index 489f86a9c5..d47e6312d9 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl @@ -99,4 +99,4 @@ ref = solve( else @test maximum(ref[4, :] - sol[4, :]) < 3e-3 end -end +end \ No newline at end of file diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl index a688dacd6a..e1a0319b57 100644 --- a/src/OrdinaryDiffEq.jl +++ b/src/OrdinaryDiffEq.jl @@ -154,9 +154,7 @@ include("caches/firk_caches.jl") include("caches/kencarp_kvaerno_caches.jl") include("caches/linear_caches.jl") include("caches/linear_nonlinear_caches.jl") -include("caches/symplectic_caches.jl") include("caches/rosenbrock_caches.jl") -include("caches/rkn_caches.jl") include("caches/adams_bashforth_moulton_caches.jl") include("caches/nordsieck_caches.jl") include("caches/bdf_caches.jl") @@ -167,12 +165,10 @@ include("caches/qprk_caches.jl") include("tableaus/low_order_rk_tableaus.jl") include("tableaus/high_order_rk_tableaus.jl") -include("tableaus/symplectic_tableaus.jl") include("tableaus/verner_tableaus.jl") include("tableaus/rosenbrock_tableaus.jl") include("tableaus/sdirk_tableaus.jl") include("tableaus/firk_tableaus.jl") -include("tableaus/rkn_tableaus.jl") include("tableaus/qprk_tableaus.jl") include("integrators/type.jl") @@ -184,8 +180,6 @@ include("initialize_dae.jl") include("wrappers.jl") include("perform_step/fixed_timestep_perform_step.jl") -include("perform_step/symplectic_perform_step.jl") -include("perform_step/rkn_perform_step.jl") include("perform_step/split_perform_step.jl") include("perform_step/linear_perform_step.jl") include("perform_step/exponential_rk_perform_step.jl") @@ -266,6 +260,19 @@ include("../lib/OrdinaryDiffEqFeagin/src/OrdinaryDiffEqFeagin.jl") using ..OrdinaryDiffEqFeagin export Feagin10, Feagin12, Feagin14 +include("../lib/OrdinaryDiffEqSymplecticRK/src/OrdinaryDiffEqSymplecticRK.jl") +using ..OrdinaryDiffEqSymplecticRK +export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, + McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, + CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 + +include("../lib/OrdinaryDiffEqRKN/src/OrdinaryDiffEqRKN.jl") +using ..OrdinaryDiffEqRKN +export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, + Nystrom5VelocityIndependent, + IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, + RKN4 + import PrecompileTools PrecompileTools.@compile_workload begin @@ -426,19 +433,10 @@ export LawsonEuler, NorsettEuler, ETD1, ETDRK2, ETDRK3, ETDRK4, HochOst4, Exp4, EPIRK4s3B, EPIRK5s3, EXPRB53s3, EPIRK5P1, EPIRK5P2, ETD2, Exprb32, Exprb43 -export SymplecticEuler, VelocityVerlet, VerletLeapfrog, PseudoVerletLeapfrog, - McAte2, Ruth3, McAte3, CandyRoz4, McAte4, McAte42, McAte5, - CalvoSanz4, Yoshida6, KahanLi6, McAte8, KahanLi8, SofSpa10 - export SHLDDRK52, SHLDDRK_2N export SplitEuler -export Nystrom4, FineRKN4, FineRKN5, Nystrom4VelocityIndependent, - Nystrom5VelocityIndependent, - IRKN3, IRKN4, DPRKN4, DPRKN5, DPRKN6, DPRKN6FM, DPRKN8, DPRKN12, ERKN4, ERKN5, ERKN7, - RKN4 - export AB3, AB4, AB5, ABM32, ABM43, ABM54 export VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5 diff --git a/src/alg_utils.jl b/src/alg_utils.jl index 4467a10eb1..ff78c80801 100644 --- a/src/alg_utils.jl +++ b/src/alg_utils.jl @@ -355,8 +355,6 @@ alg_extrapolates(alg::DImplicitEuler) = true alg_extrapolates(alg::DABDF2) = true alg_extrapolates(alg::Trapezoid) = true alg_extrapolates(alg::SDIRK22) = true -alg_extrapolates(alg::IRKN4) = true -alg_extrapolates(alg::IRKN3) = true alg_extrapolates(alg::ABDF2) = true alg_extrapolates(alg::SBDF) = true alg_extrapolates(alg::MEBDF2) = true @@ -421,42 +419,6 @@ alg_order(alg::KuttaPRK2p5) = 5 alg_order(alg::RKO65) = 5 alg_order(alg::FRK65) = 6 -alg_order(alg::SymplecticEuler) = 1 -alg_order(alg::VelocityVerlet) = 2 -alg_order(alg::VerletLeapfrog) = 2 -alg_order(alg::PseudoVerletLeapfrog) = 2 -alg_order(alg::McAte2) = 2 -alg_order(alg::Ruth3) = 3 -alg_order(alg::McAte3) = 3 -alg_order(alg::McAte4) = 4 -alg_order(alg::CandyRoz4) = 4 -alg_order(alg::CalvoSanz4) = 4 -alg_order(alg::McAte42) = 4 -alg_order(alg::McAte5) = 5 -alg_order(alg::Yoshida6) = 6 -alg_order(alg::KahanLi6) = 6 -alg_order(alg::McAte8) = 8 -alg_order(alg::KahanLi8) = 8 -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 -alg_order(alg::Nystrom5VelocityIndependent) = 5 -alg_order(alg::DPRKN4) = 4 -alg_order(alg::DPRKN5) = 5 -alg_order(alg::DPRKN6) = 6 -alg_order(alg::DPRKN6FM) = 6 -alg_order(alg::DPRKN8) = 8 -alg_order(alg::DPRKN12) = 12 -alg_order(alg::ERKN4) = 4 -alg_order(alg::ERKN5) = 5 -alg_order(alg::ERKN7) = 7 -alg_order(alg::RKN4) = 4 - alg_order(alg::Midpoint) = 2 alg_order(alg::RK4) = 4 diff --git a/src/algorithms.jl b/src/algorithms.jl index 973ebb43d4..7b7ae3ce6e 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -307,564 +307,6 @@ end ################################################################################ -# Symplectic methods - -struct SymplecticEuler <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{verlet1967computer, -title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, -author={Verlet, Loup}, -journal={Physical review}, -volume={159}, -number={1}, -pages={98}, -year={1967}, -publisher={APS} -} -""" -struct VelocityVerlet <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{verlet1967computer, -title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, -author={Verlet, Loup}, -journal={Physical review}, -volume={159}, -number={1}, -pages={98}, -year={1967}, -publisher={APS} -} -""" -struct VerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{verlet1967computer, -title={Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules}, -author={Verlet, Loup}, -journal={Physical review}, -volume={159}, -number={1}, -pages={98}, -year={1967}, -publisher={APS} -} -""" -struct PseudoVerletLeapfrog <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte2 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{ruth1983canonical, -title={A canonical integration technique}, -author={Ruth, Ronald D}, -journal={IEEE Trans. Nucl. Sci.}, -volume={30}, -number={CERN-LEP-TH-83-14}, -pages={2669--2671}, -year={1983} -} -""" -struct Ruth3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{candy1991symplectic, -title={A symplectic integration algorithm for separable Hamiltonian functions}, -author={Candy, J and Rozmus, W}, -journal={Journal of Computational Physics}, -volume={92}, -number={1}, -pages={230--256}, -year={1991}, -publisher={Elsevier} -} -""" -struct CandyRoz4 <: OrdinaryDiffEqPartitionedAlgorithm end -struct McAte4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{sanz1993symplectic, -title={Symplectic numerical methods for Hamiltonian problems}, -author={Sanz-Serna, Jes{\'u}s Maria and Calvo, Mari-Paz}, -journal={International Journal of Modern Physics C}, -volume={4}, -number={02}, -pages={385--392}, -year={1993}, -publisher={World Scientific} -} -""" -struct CalvoSanz4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte42 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1992accuracy, -title={The accuracy of symplectic integrators}, -author={McLachlan, Robert I and Atela, Pau}, -journal={Nonlinearity}, -volume={5}, -number={2}, -pages={541}, -year={1992}, -publisher={IOP Publishing} -} -""" -struct McAte5 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{yoshida1990construction, -title={Construction of higher order symplectic integrators}, -author={Yoshida, Haruo}, -journal={Physics letters A}, -volume={150}, -number={5-7}, -pages={262--268}, -year={1990}, -publisher={Elsevier} -} -""" -struct Yoshida6 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{kahan1997composition, -title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, -author={Kahan, William and Li, Ren-Cang}, -journal={Mathematics of computation}, -volume={66}, -number={219}, -pages={1089--1099}, -year={1997} -} -""" -struct KahanLi6 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{mclachlan1995numerical, -title={On the numerical integration of ordinary differential equations by symmetric composition methods}, -author={McLachlan, Robert I}, -journal={SIAM Journal on Scientific Computing}, -volume={16}, -number={1}, -pages={151--168}, -year={1995}, -publisher={SIAM} -} -""" -struct McAte8 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{kahan1997composition, -title={Composition constants for raising the orders of unconventional schemes for ordinary differential equations}, -author={Kahan, William and Li, Ren-Cang}, -journal={Mathematics of computation}, -volume={66}, -number={219}, -pages={1089--1099}, -year={1997} -} -""" -struct KahanLi8 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" -@article{sofroniou2005derivation, -title={Derivation of symmetric composition constants for symmetric integrators}, -author={Sofroniou, Mark and Spaletta, Giulia}, -journal={Optimization Methods and Software}, -volume={20}, -number={4-5}, -pages={597--613}, -year={2005}, -publisher={Taylor \\& Francis} -} -""" -struct SofSpa10 <: OrdinaryDiffEqPartitionedAlgorithm end - -# Nyström methods - -""" - IRKN3 - -Improved Runge-Kutta-Nyström method of order three, which minimizes the amount of evaluated functions in each step. Fixed time steps only. - -Second order ODE should not depend on the first derivative. - -## References - -@article{rabiei2012numerical, -title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, -author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, -publisher={Citeseer} -} -""" -struct IRKN3 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - Nystrom4 - -A 4th 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 -[`Nystrom4VelocityIndependent`](@ref) to increase performance. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -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() - -A 5th 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 FineRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - Nystrom4VelocityIdependent - -A 4th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). - -More efficient then [`Nystrom4`](@ref) on velocity independent problems, since less evaluations are needed. - -Fixed time steps only. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom4VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - IRKN4 - -Improves Runge-Kutta-Nyström method of order four, which minimizes the amount of evaluated functions in each step. Fixed time steps only. - -Second order ODE should not be dependent on the first derivative. - -Recommended for smooth problems with expensive functions to evaluate. - -## References - -@article{rabiei2012numerical, -title={Numerical Solution of Second-Order Ordinary Differential Equations by Improved Runge-Kutta Nystrom Method}, -author={Rabiei, Faranak and Ismail, Fudziah and Norazak, S and Emadi, Saeid}, -publisher={Citeseer} -} -""" -struct IRKN4 <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - Nystrom5VelocityIndependent - -A 5th order explicit Runkge-Kutta-Nyström method. Used directly on second order ODEs, where the acceleration is independent from velocity (ODE Problem is not dependent on the first derivative). -Fixed time steps only. - -## References - -E. Hairer, S.P. Norsett, G. Wanner, (1993) Solving Ordinary Differential Equations I. -Nonstiff Problems. 2nd Edition. Springer Series in Computational Mathematics, -Springer-Verlag. -""" -struct Nystrom5VelocityIndependent <: OrdinaryDiffEqPartitionedAlgorithm end - -""" - DPRKN4 - -4th order explicit Runge-Kutta-Nyström methods. The second order ODE should not depend on the first derivative. - -## References - -@article{Dormand1987FamiliesOR, -title={Families of Runge-Kutta-Nystrom Formulae}, -author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, -journal={Ima Journal of Numerical Analysis}, -year={1987}, -volume={7}, -pages={235-250} -} -""" -struct DPRKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN5 - -5th order explicit Runge-Kutta-Nyström mehod. The second order ODE should not depend on the first derivative. - -## References - -@article{Bettis1973ARN, -title={A Runge-Kutta Nystrom algorithm}, -author={Dale G. Bettis}, -journal={Celestial mechanics}, -year={1973}, -volume={8}, -pages={229-233}, -publisher={Springer} -} -""" -struct DPRKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN6 - -6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. Free 6th order interpolant. - -## References - -@article{dormand1987runge, -title={Runge-kutta-nystrom triples}, -author={Dormand, JR and Prince, PJ}, -journal={Computers \\& Mathematics with Applications}, -volume={13}, -number={12}, -pages={937--949}, -year={1987}, -publisher={Elsevier} -} -""" -struct DPRKN6 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN6FM - -6th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Compared to [`DPRKN6`](@ref), this method has smaller truncation error coefficients which leads to performance gain -when only the main solution points are considered. - -## References - -@article{Dormand1987FamiliesOR, -title={Families of Runge-Kutta-Nystrom Formulae}, -author={J. R. Dormand and Moawwad E. A. El-Mikkawy and P. J. Prince}, -journal={Ima Journal of Numerical Analysis}, -year={1987}, -volume={7}, -pages={235-250} -} -""" -struct DPRKN6FM <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN8 - -8th order explicit Runge-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Not as efficient as [`DPRKN12`](@ref) when high accuracy is needed, however this solver is competitive with -[`DPRKN6`](@ref) at lax tolerances and, depending on the problem, might be a good option between performance and accuracy. - -## References - -@article{dormand1987high, -title={High-order embedded Runge-Kutta-Nystrom formulae}, -author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, -journal={IMA Journal of Numerical Analysis}, -volume={7}, -number={4}, -pages={423--430}, -year={1987}, -publisher={Oxford University Press} -} -""" -struct DPRKN8 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - DPRKN12 - -12th order explicit Rugne-Kutta-Nyström method. The second order ODE should not depend on the first derivative. - -Most efficient when high accuracy is needed. - -## References - -@article{dormand1987high, -title={High-order embedded Runge-Kutta-Nystrom formulae}, -author={Dormand, JR and El-Mikkawy, MEA and Prince, PJ}, -journal={IMA Journal of Numerical Analysis}, -volume={7}, -number={4}, -pages={423--430}, -year={1987}, -publisher={Oxford University Press} -} -""" -struct DPRKN12 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN4 - -Embedded 4(3) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic problems. - -## References - -@article{demba2017embedded, -title={An Embedded 4 (3) Pair of Explicit Trigonometrically-Fitted Runge-Kutta-Nystr{\"o}m Method for Solving Periodic Initial Value Problems}, -author={Demba, MA and Senu, N and Ismail, F}, -journal={Applied Mathematical Sciences}, -volume={11}, -number={17}, -pages={819--838}, -year={2017} -} -""" -struct ERKN4 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN5 - -Embedded 5(4) pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic problems. - -## References - -@article{demba20165, -title={A 5 (4) Embedded Pair of Explicit Trigonometrically-Fitted Runge--Kutta--Nystr{\"o}m Methods for the Numerical Solution of Oscillatory Initial Value Problems}, -author={Demba, Musa A and Senu, Norazak and Ismail, Fudziah}, -journal={Mathematical and Computational Applications}, -volume={21}, -number={4}, -pages={46}, -year={2016}, -publisher={Multidisciplinary Digital Publishing Institute} -} -""" -struct ERKN5 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" - ERKN7 - -Embedded pair of explicit Runge-Kutta-Nyström methods. Integrates the periodic properties of the harmonic oscillator exactly. - -The second order ODE should not depend on the first derivative. - -Uses adaptive step size control. This method is extra efficient on periodic Problems. - -## References - -@article{SimosOnHO, -title={On high order Runge-Kutta-Nystr{\"o}m pairs}, -author={Theodore E. Simos and Ch. Tsitouras}, -journal={J. Comput. Appl. Math.}, -volume={400}, -pages={113753} -} -""" -struct ERKN7 <: OrdinaryDiffEqAdaptivePartitionedAlgorithm end - -""" -3 stage fourth order Runge-Kutta Nystrom method to solve second order linear inhomogeneous IVPs. - -Does not include an adaptive method. Solves for for d-dimensional differential systems of second order linear inhomogeneous equations. - -!!! warn - - This method is only fourth order for these systems, the method is second order otherwise! - -## References - -@article{MONTIJANO2024115533, -title = {Explicit Runge–Kutta–Nyström methods for the numerical solution of second order linear inhomogeneous IVPs}, -author = {J.I. Montijano and L. Rández and M. Calvo}, -journal = {Journal of Computational and Applied Mathematics}, -volume = {438}, -pages = {115533}, -year = {2024}, -} -""" -struct RKN4 <: OrdinaryDiffEqAlgorithm end - -################################################################################ - # Adams Bashforth and Adams moulton methods """ @@ -2847,7 +2289,7 @@ function PDIRK44(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{tru end ### Algorithm Groups -const MultistepAlgorithms = Union{IRKN3, IRKN4, +const MultistepAlgorithms = Union{ ABDF2, AB3, AB4, AB5, ABM32, ABM43, ABM54} diff --git a/src/dense/interpolants.jl b/src/dense/interpolants.jl index 7356b8865d..ba99d9f4aa 100644 --- a/src/dense/interpolants.jl +++ b/src/dense/interpolants.jl @@ -2725,145 +2725,4 @@ end end """ -""" -@def dprkn6unpack begin - if cache isa OrdinaryDiffEqMutableCache - @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache.tab - else - @unpack r14, r13, r12, r11, r10, r34, r33, r32, r31, r44, r43, r42, r41, r54, r53, r52, r51, r64, r63, r62, r61, rp14, rp13, rp12, rp11, rp10, rp34, rp33, rp32, rp31, rp44, rp43, rp42, rp41, rp54, rp53, rp52, rp51, rp64, rp63, rp62, rp61 = cache - end -end - -@def dprkn6pre0 begin - @dprkn6unpack - b1Θ = @evalpoly(Θ, r10, r11, r12, r13, r14) - b3Θ = Θ * @evalpoly(Θ, r31, r32, r33, r34) - b4Θ = Θ * @evalpoly(Θ, r41, r42, r43, r44) - b5Θ = Θ * @evalpoly(Θ, r51, r52, r53, r54) - b6Θ = Θ * @evalpoly(Θ, r61, r62, r63, r64) - - bp1Θ = @evalpoly(Θ, rp10, rp11, rp12, rp13, rp14) - bp3Θ = Θ * @evalpoly(Θ, rp31, rp32, rp33, rp34) - bp4Θ = Θ * @evalpoly(Θ, rp41, rp42, rp43, rp44) - bp5Θ = Θ * @evalpoly(Θ, rp51, rp52, rp53, rp54) - bp6Θ = Θ * @evalpoly(Θ, rp61, rp62, rp63, rp64) - - kk1, kk2, kk3 = k - k1, k2 = kk1.x - k3, k4 = kk2.x - k5, k6 = kk3.x - - duprev, uprev = y₀.x - dtsq = dt^2 -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - return ArrayPartition( - duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6), - uprev + - dt * Θ * - (duprev + - dt * Θ * (b1Θ * k1 + b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6))) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - return ArrayPartition( - duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]), - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs]))) -end - -@muladd function _ode_interpolant(Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs::Number, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - halfsize = length(y₀) ÷ 2 - if idxs <= halfsize - duprev[idxs] + - dt * Θ * - (bp1Θ * k1[idxs] + bp3Θ * k3[idxs] + - bp4Θ * k4[idxs] + bp5Θ * k5[idxs] + bp6Θ * k6[idxs]) - else - idxs = idxs - halfsize - uprev[idxs] + - dt * Θ * - (duprev[idxs] + - dt * Θ * - (b1Θ * k1[idxs] + - b3Θ * k3[idxs] + - b4Θ * k4[idxs] + b5Θ * k5[idxs] + b6Θ * k6[idxs])) - end -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, - idxs::Nothing, T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - @inbounds @.. broadcast=false out.x[2]=uprev + - dt * Θ * - (duprev + - dt * Θ * - (b1Θ * k1 + - b3Θ * k3 + - b4Θ * k4 + b5Θ * k5 + b6Θ * k6)) - @inbounds @.. broadcast=false out.x[1]=duprev + - dt * Θ * - (bp1Θ * k1 + bp3Θ * k3 + - bp4Θ * k4 + bp5Θ * k5 + bp6Θ * k6) - #for i in eachindex(out.x[1]) - # out.x[2][i] = uprev[i] + dt*Θ*(duprev[i] + dt*Θ*(b1Θ*k1[i] + - # b3Θ*k3[i] + - # b4Θ*k4[i] + b5Θ*k5[i] + b6Θ*k6[i])) - # out.x[1][i] = duprev[i] + dt*Θ*(bp1Θ*k1[i] + bp3Θ*k3[i] + - # bp4Θ*k4[i] + bp5Θ*k5[i] + bp6Θ*k6[i]) - #end - out -end - -@muladd function _ode_interpolant!(out, Θ, dt, y₀, y₁, k, - cache::Union{DPRKN6ConstantCache, DPRKN6Cache}, idxs, - T::Type{Val{0}}, differential_vars::Nothing) - @dprkn6pre0 - halfsize = length(y₀) ÷ 2 - isfirsthalf = idxs .<= halfsize - secondhalf = idxs .> halfsize - firstidxs = idxs[isfirsthalf] - secondidxs_shifted = idxs[secondhalf] - secondidxs = secondidxs_shifted .- halfsize - - @views @.. broadcast=false out[secondhalf]=uprev[secondidxs] + - dt * Θ * - (duprev[secondidxs] + - dt * Θ * - (b1Θ * k1[secondidxs] + - b3Θ * k3[secondidxs] + - b4Θ * k4[secondidxs] + - b5Θ * k5[secondidxs] + - b6Θ * k6[secondidxs])) - @views @.. broadcast=false out[isfirsthalf]=duprev[firstidxs] + - dt * Θ * - (bp1Θ * k1[firstidxs] + - bp3Θ * k3[firstidxs] + - bp4Θ * k4[firstidxs] + - bp5Θ * k5[firstidxs] + - bp6Θ * k6[firstidxs]) - out -end +""" \ No newline at end of file diff --git a/src/interp_func.jl b/src/interp_func.jl index 14c95b7aba..4265219b7d 100644 --- a/src/interp_func.jl +++ b/src/interp_func.jl @@ -148,13 +148,6 @@ function DiffEqBase.interp_summary(::Type{cacheType}, caches = fieldtype(cacheType, :caches) join([DiffEqBase.interp_summary(ct, dense) for ct in fieldtypes(caches)], ", ") end -function DiffEqBase.interp_summary(::Type{cacheType}, - dense::Bool) where { - cacheType <: - Union{DPRKN6ConstantCache, - DPRKN6Cache}} - dense ? "specialized 6th order interpolation" : "1st order linear" -end function (interp::InterpolationData)(tvals, idxs, deriv, p, continuity::Symbol = :left) ode_interpolation(tvals, interp, idxs, deriv, p, continuity) diff --git a/test/runtests.jl b/test/runtests.jl index 14aa5a0703..b5a0fd333e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,12 @@ function activate_feagin() Pkg.instantiate() end +function activate_symplectic_rk() + Pkg.activate("../lib/OrdinaryDiffEqSymplecticRK") + Pkg.develop(PackageSpec(path = dirname(@__DIR__))) + Pkg.instantiate() +end + #Start Test Script @time begin @@ -199,10 +205,13 @@ end @time @safetestset "FIRK Tests" include("algconvergence/ode_firk_tests.jl") @time @safetestset "Linear-Nonlinear Methods Tests" include("algconvergence/linear_nonlinear_convergence_tests.jl") @time @safetestset "Linear-Nonlinear Krylov Methods Tests" include("algconvergence/linear_nonlinear_krylov_tests.jl") - @time @safetestset "Symplectic Tests" include("algconvergence/symplectic_tests.jl") @time @safetestset "Quadruple precision Runge-Kutta Tests" include("algconvergence/ode_quadruple_precision_tests.jl") end + if !is_APPVEYOR && GROUP == "Symplectic" + @time @safetestset "Symplectic Tests" include("../lib/OrdinaryDiffEqSymplecticRK/test/symplectic_tests.jl") + end + if !is_APPVEYOR && GROUP == "Extrapolation" @time @safetestset "Extrapolation Tests" include("../lib/OrdinaryDiffEqExtrapolation/test/runtests.jl") end