Merge pull request #149 from SciML/hr/imaginary_stability_interval

implement `imaginary_stability_interval`

Author: ChrisRackauckas

Project.toml

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

src/DiffEqDevTools.jl

@@ -45,7 +45,7 @@ export test_convergence, analyticless_test_convergence, appxtrue!, appxtrue
 export get_sample_errors
 
 #Tab Functions
-export stability_region, residual_order_condition, check_tableau
+export stability_region, residual_order_condition, check_tableau, imaginary_stability_interval
 
 #Tableaus
 export deduce_Butcher_tableau

src/tableau_info.jl

@@ -27,6 +27,7 @@ end
 `stability_region(tab::ODERKTableau; 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...)
     residual! = function (resid, x)
@@ -36,6 +37,24 @@ function stability_region(tab::ODERKTableau; initial_guess = -3.0, kw...)
     sol.zero[1]
 end
 
+"""
+    imaginary_stability_interval(tab::ODERKTableau;
+                                 initial_guess = length(tab) - 1)
+
+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(tab::ODERKTableau;
+                                      initial_guess = length(tab) - one(eltype(tab.A)),
+                                      kw...)
+    residual! = function (resid, x)
+        resid[1] = abs(stability_region(im * x[1], tab)) - 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

@@ -3,3 +3,9 @@ using DiffEqDevTools, Test
 @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(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