Added Rosenbrock methods

lib/OrdinaryDiffEqBDF/Project.toml

@@ -0,0 +1,22 @@
+name = "OrdinaryDiffEqBDF"
+uuid = "34d55129-80d4-4dd1-bb14-71372a5ab94f"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "1.0.0"
+
+[deps]
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+
+[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"]

lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl

@@ -0,0 +1,25 @@
+module OrdinaryDiffEqBDF
+
+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("bdf_caches.jl")
+include("bdf_perform_step.jl")
+
+export SBDF, ABDF2, QNDF1, QNDF2, QNDF, MEBDF2
+
+end

lib/OrdinaryDiffEqBDF/src/alg_utils.jl

@@ -0,0 +1,23 @@
+alg_extrapolates(alg::ABDF2) = true
+alg_extrapolates(alg::SBDF) = true
+alg_extrapolates(alg::MEBDF2) = true
+
+alg_order(alg::ABDF2) = 2
+alg_order(alg::SBDF) = alg.order
+alg_order(alg::QNDF1) = 1
+alg_order(alg::MEBDF2) = 2
+alg_order(alg::FBDF) = 1 #dummy value
+alg_order(alg::QNDF2) = 2
+alg_order(alg::QNDF) = 1 #dummy value
+
+get_current_alg_order(alg::QNDF, cache) = cache.order
+get_current_alg_order(alg::FBDF, cache) = cache.order
+
+get_current_adaptive_order(alg::FBDF, cache) = cache.order
+
+qsteady_max_default(alg::QNDF) = 2 // 1
+qsteady_max_default(alg::QNDF2) = 2 // 1
+qsteady_max_default(alg::QNDF1) = 2 // 1
+
+qsteady_min_default(alg::FBDF) = 9 // 10
+qsteady_max_default(alg::FBDF) = 2 // 1

lib/OrdinaryDiffEqBDF/src/algorithms.jl

@@ -0,0 +1,223 @@
+"""
+Uri M. Ascher, Steven J. Ruuth, Brian T. R. Wetton. Implicit-Explicit Methods for Time-
+Dependent Partial Differential Equations. 1995 Society for Industrial and Applied Mathematics
+Journal on Numerical Analysis, 32(3), pp 797-823, 1995. doi: https://doi.org/10.1137/0732037
+"""
+struct SBDF{CS, AD, F, F2, P, FDT, ST, CJ, K, T} <:
+       OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    κ::K
+    tol::T
+    extrapolant::Symbol
+    order::Int
+    ark::Bool
+end
+
+function SBDF(order; chunk_size = Val{0}(), autodiff = Val{true}(),
+        standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        tol = nothing,
+        extrapolant = :linear, ark = false)
+    SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(tol)}(linsolve,
+        nlsolve,
+        precs,
+        κ,
+        tol,
+        extrapolant,
+        order,
+        ark)
+end
+
+# All keyword form needed for remake
+function SBDF(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+        concrete_jac = nothing, diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        tol = nothing,
+        extrapolant = :linear,
+        order, ark = false)
+    SBDF{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(tol)}(linsolve,
+        nlsolve,
+        precs,
+        κ,
+        tol,
+        extrapolant,
+        order,
+        ark)
+end
+
+"""
+E. Alberdi Celayaa, J. J. Anza Aguirrezabalab, P. Chatzipantelidisc. Implementation of
+an Adaptive BDF2 Formula and Comparison with The MATLAB Ode15s. Procedia Computer Science,
+29, pp 1014-1026, 2014. doi: https://doi.org/10.1016/j.procs.2014.05.091
+
+ABDF2: Multistep Method
+An adaptive order 2 L-stable fixed leading coefficient multistep BDF method.
+"""
+struct ABDF2{CS, AD, F, F2, P, FDT, ST, CJ, K, T, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    κ::K
+    tol::T
+    smooth_est::Bool
+    extrapolant::Symbol
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+function ABDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(),
+        concrete_jac = nothing, diff_type = Val{:forward},
+        κ = nothing, tol = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
+        nlsolve = NLNewton(),
+        smooth_est = true, extrapolant = :linear,
+        controller = :Standard, step_limiter! = trivial_limiter!)
+    ABDF2{
+        _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(κ), typeof(tol), typeof(step_limiter!)}(linsolve, nlsolve, precs, κ, tol,
+        smooth_est, extrapolant, controller, step_limiter!)
+end
+
+"""
+QNDF1: Multistep Method
+An adaptive order 1 quasi-constant timestep L-stable numerical differentiation function (NDF) method.
+Optional parameter kappa defaults to Shampine's accuracy-optimal -0.1850.
+
+See also `QNDF`.
+"""
+struct QNDF1{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    extrapolant::Symbol
+    kappa::κType
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function QNDF1(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+        concrete_jac = nothing, diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        extrapolant = :linear, kappa = -0.1850,
+        controller = :Standard, step_limiter! = trivial_limiter!)
+    QNDF1{
+        _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(kappa), typeof(step_limiter!)}(linsolve,
+        nlsolve,
+        precs,
+        extrapolant,
+        kappa,
+        controller,
+        step_limiter!)
+end
+
+"""
+QNDF2: Multistep Method
+An adaptive order 2 quasi-constant timestep L-stable numerical differentiation function (NDF) method.
+
+See also `QNDF`.
+"""
+struct QNDF2{CS, AD, F, F2, P, FDT, ST, CJ, κType, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    extrapolant::Symbol
+    kappa::κType
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function QNDF2(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(),
+        concrete_jac = nothing, diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        extrapolant = :linear, kappa = -1 // 9,
+        controller = :Standard, step_limiter! = trivial_limiter!)
+    QNDF2{
+        _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve), typeof(nlsolve),
+        typeof(precs), diff_type, _unwrap_val(standardtag), _unwrap_val(concrete_jac),
+        typeof(kappa), typeof(step_limiter!)}(linsolve,
+        nlsolve,
+        precs,
+        extrapolant,
+        kappa,
+        controller,
+        step_limiter!)
+end
+
+"""
+QNDF: Multistep Method
+An adaptive order quasi-constant timestep NDF method.
+Utilizes Shampine's accuracy-optimal kappa values as defaults (has a keyword argument for a tuple of kappa coefficients).
+
+@article{shampine1997matlab,
+title={The matlab ode suite},
+author={Shampine, Lawrence F and Reichelt, Mark W},
+journal={SIAM journal on scientific computing},
+volume={18},
+number={1},
+pages={1--22},
+year={1997},
+publisher={SIAM}
+}
+"""
+struct QNDF{MO, CS, AD, F, F2, P, FDT, ST, CJ, K, T, κType, StepLimiter} <:
+       OrdinaryDiffEqNewtonAdaptiveAlgorithm{CS, AD, FDT, ST, CJ}
+    max_order::Val{MO}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    κ::K
+    tol::T
+    extrapolant::Symbol
+    kappa::κType
+    controller::Symbol
+    step_limiter!::StepLimiter
+end
+
+function QNDF(; max_order::Val{MO} = Val{5}(), chunk_size = Val{0}(),
+        autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing,
+        diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(), κ = nothing,
+        tol = nothing,
+        extrapolant = :linear, kappa = promote(-0.1850, -1 // 9, -0.0823, -0.0415, 0),
+        controller = :Standard, step_limiter! = trivial_limiter!) where {MO}
+    QNDF{MO, _unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac),
+        typeof(κ), typeof(tol), typeof(kappa), typeof(step_limiter!)}(
+        max_order, linsolve, nlsolve, precs, κ, tol,
+        extrapolant, kappa, controller, step_limiter!)
+end
+
+"""
+MEBDF2: Multistep Method
+The second order Modified Extended BDF method, which has improved stability properties over the standard BDF.
+Fixed timestep only.
+"""
+struct MEBDF2{CS, AD, F, F2, P, FDT, ST, CJ} <:
+       OrdinaryDiffEqNewtonAlgorithm{CS, AD, FDT, ST, CJ}
+    linsolve::F
+    nlsolve::F2
+    precs::P
+    extrapolant::Symbol
+end
+function MEBDF2(; chunk_size = Val{0}(), autodiff = true, standardtag = Val{true}(),
+        concrete_jac = nothing, diff_type = Val{:forward},
+        linsolve = nothing, precs = DEFAULT_PRECS, nlsolve = NLNewton(),
+        extrapolant = :constant)
+    MEBDF2{_unwrap_val(chunk_size), _unwrap_val(autodiff), typeof(linsolve),
+        typeof(nlsolve), typeof(precs), diff_type, _unwrap_val(standardtag),
+        _unwrap_val(concrete_jac)}(linsolve,
+        nlsolve,
+        precs,
+        extrapolant)
+end

src/perform_step/bdf_perform_step.jl → lib/OrdinaryDiffEqBDF/src/bdf_perform_step.jl

@@ -1342,4 +1342,4 @@ function perform_step!(integrator, cache::FBDFCache{max_order},
 
     f(integrator.fsallast, u, p, tdt)
     integrator.stats.nf += 1
-end
+end

lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl

@@ -1,8 +1,8 @@
 module OrdinaryDiffEqDefault
 
-using OrdinaryDiffEq: Vern7, Vern8, Vern9, Vern6, Tsit5, Rosenbrock23, Rodas5P, FBDF,
-    alg_stability_size, beta2_default, beta1_default, AutoSwitchCache, ODEIntegrator,
-    CompositeAlgorithm, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, AutoAlgSwitch
+using OrdinaryDiffEq: Vern7, Vern8, Vern9, Vern6, Tsit5, FBDF,
+                      alg_stability_size, beta2_default, beta1_default, AutoSwitchCache, ODEIntegrator,
+                      CompositeAlgorithm, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, AutoAlgSwitch
 import OrdinaryDiffEq: is_mass_matrix_alg, default_autoswitch
 import LinearSolve
 using LinearAlgebra: I, isdiag

lib/OrdinaryDiffEqExplicitRK/Project.toml

@@ -0,0 +1,24 @@
+name = "OrdinaryDiffEqExplicitRK"
+uuid = "a6f9a3f1-1c5b-4457-b8d8-bad5275751bf"
+authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
+version = "0.1.0"
+
+[deps]
+FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+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"]

lib/OrdinaryDiffEqExplicitRK/src/OrdinaryDiffEqExplicitRK.jl

@@ -0,0 +1,22 @@
+module OrdinaryDiffEqExplicitRK
+
+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("explicit_rk_caches.jl")
+include("explicit_rk_perform_step.jl")
+
+end

lib/OrdinaryDiffEqExplicitRK/src/alg_utils.jl

@@ -0,0 +1,7 @@
+isfsal(tab::DiffEqBase.ExplicitRKTableau) = tab.fsal
+
+alg_order(alg::ExplicitRK) = alg.tableau.order
+
+alg_adaptive_order(alg::ExplicitRK) = alg.tableau.adaptiveorder
+
+alg_stability_size(alg::ExplicitRK) = alg.tableau.stability_size

lib/OrdinaryDiffEqExplicitRK/src/algorithms.jl

@@ -0,0 +1,6 @@
+struct ExplicitRK{TabType} <: OrdinaryDiffEqAdaptiveAlgorithm
+    tableau::TabType
+end
+ExplicitRK(; tableau = ODE_DEFAULT_TABLEAU) = ExplicitRK(tableau)
+
+TruncatedStacktraces.@truncate_stacktrace ExplicitRK