Merge pull request #70 from SciML/myb/stability_ply

Optimize stability polynomial evaluation and add order star plots

Project.toml

@@ -1,7 +1,7 @@
 name = "DiffEqDevTools"
 uuid = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "2.23.0"
+version = "2.24.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"

src/plotrecipes.jl

@@ -53,12 +53,21 @@ end
   errors,times
 end
 
-@recipe function f(tab::ODERKTableau;dx=1/100,dy=1/100,xlim=[-6,1],ylim=[-5,5])
-  x = xlim[1]:dx:xlim[2]
-  y = ylim[1]:dy:ylim[2]
-  f = (u,v)-> abs(stability_region(u+v*im,tab))<1
+@recipe function f(tab::ODERKTableau;dx=1/100,dy=1/100,order_star=false,embedded=false)
+  xlims = get(plotattributes, :xlims, (-6,1))
+  ylims = get(plotattributes, :ylims, (-5,5))
+  x = xlims[1]:dx:xlims[2]
+  y = ylims[1]:dy:ylims[2]
+
+  if order_star
+    f = (u,v)-> abs(stability_region(u+v*im,tab; embedded=embedded)/exp(u+v*im)) < 1
+  else
+    f = (u,v)-> abs(stability_region(u+v*im,tab; embedded=embedded)) < 1
+  end
   seriestype --> :contour
   fill --> true
-  cbar --> false
+  colorbar --> false
+  seriescolor --> :grays
+  aspect_ratio --> 1
   x,y,f
 end

src/tableau_info.jl

@@ -12,10 +12,17 @@ Base.length(tab::ODERKTableau) = tab.stages
 Calculates the stability function from the tableau at `z`. Stable if <1.
 
 ```math
-r(z) = \\frac{\\det(I-zA+zeb^T)}{\\det(I-zA)}
+r(z) = 1 + z bᵀ(I - zA)⁻¹ 𝟙
 ```
+where 𝟙 denotes a vector of ones.
 """
-stability_region(z,tab::ODERKTableau) = det(Matrix{Float64}(I,tab.stages,tab.stages)- z*tab.A + z*ones(tab.stages)*tab.α')/det(Matrix{Float64}(I,tab.stages,tab.stages)-z*tab.A)
+function stability_region(z,tab::ODERKTableau; embedded=false)
+  A, c = tab.A, tab.c
+  b = embedded ? tab.αEEst : tab.α
+  𝟙 = ones(eltype(A), length(b))
+  stages = (I - z*A) \ 𝟙
+  1 + z * (transpose(b) * stages)
+end
 
 """
 `stability_region(tab::ODERKTableau; initial_guess=-3.0)`