Add a method to compute the length of stability regions

Author: YingboMa

REQUIRE

@@ -3,6 +3,7 @@ RecipesBase 0.1.0
 DiffEqBase 3.0.0
 RecursiveArrayTools 0.4.0
 DiffEqPDEBase 0.4.0
+NLsolve 0.14.1
 DiffEqMonteCarlo
 DiffEqNoiseProcess
 Juno

src/tableau_info.jl

@@ -1,3 +1,4 @@
+using NLsolve
 """
 `Base.length(tab::ODERKTableau)`
 
@@ -15,3 +16,16 @@ r(z) = \\frac{\\det(I-zA+zeb^T)}{\\det(I-zA)}
 ```
 """
 stability_region(z,tab::ODERKTableau) = det(eye(tab.stages)- z*tab.A + z*ones(tab.stages)*tab.α')/det(eye(tab.stages)-z*tab.A)
+
+"""
+`stability_region(tab::ODERKTableau)`
+
+Calculates the length of the stability region in the real axis.
+"""
+function stability_region(tab::ODERKTableau; initial_guess=-3.0)
+  residual! = function (resid, x)
+    resid[1] = abs(stability_region(x[1], tab)) - 1
+  end
+  sol = nlsolve(residual!, [initial_guess])
+  sol.zero[1]
+end

test/stability_region_test.jl

@@ -0,0 +1,5 @@
+using DiffEqDevTools, Base.Test
+
+@test stability_region(constructDormandPrince6(), initial_guess=-3.5) + 3.95413 < 1e-3
+@test stability_region(constructTsitourasPapakostas6(), initial_guess=-3.5) + 3.95413 < 1e-3
+@test stability_region(constructRadauIIA5(), initial_guess=12.) - 11.84 < 1e-2