Merge pull request #150 from SciML/hr/stability_region_for_methods

allow computing the stability function directly from an algorithm

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "DiffEqDevTools"
 uuid = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "2.46.0"
+version = "2.47.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

src/DiffEqDevTools.jl

@@ -1,5 +1,6 @@
 module DiffEqDevTools
 
+using DiffEqBase: AbstractODEAlgorithm
 using DiffEqBase, RecipesBase, RecursiveArrayTools, DiffEqNoiseProcess, StructArrays
 using NLsolve, LinearAlgebra, RootedTrees
 

src/ode_tableaus.jl

@@ -1,5 +1,5 @@
 """
-    deduce_Butcher_tableau(erk, T=Float64)
+    deduce_Butcher_tableau(erk, T = Float64)
 
 Deduce and return the Butcher coefficients `A, b, c` by solving some
 specific ordinary differential equations using the explicit Runge-Kutta

src/tableau_info.jl

@@ -6,9 +6,12 @@ Defines the length of a Runge-Kutta method to be the number of stages.
 Base.length(tab::ODERKTableau) = tab.stages
 
 """
-`stability_region(z,tab::ODERKTableau)`
+`stability_region(z, tab::ODERKTableau; embedded = false)`
 
 Calculates the stability function from the tableau at `z`. Stable if <1.
+If `embedded = true`, the stability function is calculated for the embedded
+method. Otherwise, the stability function is calculated for the main method
+(default).
 
 ```math
 r(z) = 1 + z bᵀ(I - zA)⁻¹ e
@@ -24,14 +27,49 @@ function stability_region(z, tab::ODERKTableau; embedded = false)
 end
 
 """
-`stability_region(tab::ODERKTableau; initial_guess=-3.0)`
+    stability_region(z, alg::AbstractODEAlgorithm)
+
+Calculates the stability function from the algorithm `alg` at `z`.
+The stability region of a possible embedded method cannot be calculated
+using this method.
+
+If you use an implicit method, you may run into convergence issues when
+the value of `z` is outside of the stability region, e.g.,
+
+```julia-repl
+julia> typemin(Float64)
+-Inf
+
+julia> stability_region(typemin(Float64), ImplicitEuler())
+┌ Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt.
+
+julia> nextfloat(typemin(Float64))
+-1.7976931348623157e308
+
+julia> stability_region(nextfloat(typemin(Float64)), ImplicitEuler())
+0.0
+```
+"""
+function stability_region(z, alg::AbstractODEAlgorithm)
+    u0 = one(z)
+    tspan = (zero(real(z)), one(real(z)))
+    ode = ODEProblem{false}((u, p, t) -> z * u, u0, tspan)
+    integrator = init(ode, alg; adaptive = false, dt = 1)
+    step!(integrator)
+    return integrator.u[end]
+end
+
+"""
+`stability_region(tab_or_alg::Union{ODERKTableau, AbstractODEAlgorithm};
+                  initial_guess=-3.0)`
 
 Calculates the length of the stability region in the real axis.
 See also [`imaginary_stability_interval`](@ref).
 """
-function stability_region(tab::ODERKTableau; initial_guess = -3.0, kw...)
+function stability_region(tab_or_alg::Union{ODERKTableau, AbstractODEAlgorithm};
+                          initial_guess = -3.0, kw...)
     residual! = function (resid, x)
-        resid[1] = abs(stability_region(x[1], tab)) - 1
+        resid[1] = abs(stability_region(x[1], tab_or_alg)) - 1
     end
     sol = nlsolve(residual!, [initial_guess]; kw...)
     sol.zero[1]
@@ -55,6 +93,24 @@ function imaginary_stability_interval(tab::ODERKTableau;
     sol.zero[1]
 end
 
+"""
+    imaginary_stability_interval(alg::ODERKTableau;
+                                 initial_guess = 20.0)
+
+Calculates the length of the imaginary stability interval, i.e.,
+the size of the stability region on the imaginary axis.
+See also [`stability_region`](@ref).
+"""
+function imaginary_stability_interval(alg::AbstractODEAlgorithm;
+                                      initial_guess = 20.0,
+                                      kw...)
+    residual! = function (resid, x)
+        resid[1] = abs(stability_region(im * x[1], alg)) - 1
+    end
+    sol = nlsolve(residual!, [initial_guess]; kw...)
+    sol.zero[1]
+end
+
 function RootedTrees.residual_order_condition(tab::ODERKTableau, order::Int,
         reducer = nothing, mapper = x -> x^2;
         embedded = false)

test/stability_region_test.jl

@@ -1,11 +1,21 @@
 using DiffEqDevTools, Test
+using OrdinaryDiffEq
 
 @test stability_region(constructDormandPrince6(), initial_guess = -3.5)≈-3.95413 rtol=1e-3
 @test stability_region(constructTsitourasPapakostas6(), initial_guess = -3.5)≈-3.95413 rtol=1e-3
 @test stability_region(constructRadauIIA5(), initial_guess = 12.0)≈11.84 rtol=1e-2
 
+@test @inferred(stability_region(constructTsitouras5()))≈@inferred(stability_region(Tsit5()))
+@test @inferred(stability_region(constructSSPRK104()))≈@inferred(stability_region(SSPRK104()))
+@test @inferred(stability_region(constructImplicitEuler()))≈@inferred(stability_region(ImplicitEuler()))
+@test @inferred(abs(stability_region(nextfloat(typemin(Float64)), ImplicitEuler()))) < eps(Float64)
+
 @test @inferred(imaginary_stability_interval(constructSSPRK33()))≈sqrt(3)
 @test @inferred(imaginary_stability_interval(constructSSPRK33(Float32)))≈sqrt(3.0f0)
 @test @inferred(imaginary_stability_interval(constructKutta3()))≈sqrt(3)
 @test @inferred(imaginary_stability_interval(constructKutta3(Float32)))≈sqrt(3.0f0)
 @test @inferred(imaginary_stability_interval(constructRK4()))≈2.8284271247
+
+@test @inferred(imaginary_stability_interval(SSPRK33()))≈sqrt(3)
+@test @inferred(imaginary_stability_interval(SSPRK33(), initial_guess = 5.0f0))≈sqrt(3.0f0)
+@test @inferred(imaginary_stability_interval(RK4()))≈2.8284271247