Merge pull request #6 from YingboMa/master

Add a method to compute the length of the stability region

Author: ChrisRackauckas

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; initial_guess=-3.0)`
+
+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/runtests.jl

@@ -7,3 +7,4 @@ using Base.Test
 @time @testset "Analyticless Convergence Tests" begin include("analyticless_convergence_tests.jl") end
 @time @testset "ODE Tableau Convergence Tests" begin include("ode_tableau_convergence_tests.jl") end ## Windows 32-bit fails on Butcher62 convergence test
 @time @testset "Analyticless Stochastic WP" begin include("analyticless_stochastic_wp.jl") end
+@time @testset "Stability Region Tests" begin include("stability_region_test.jl") end

test/stability_region_test.jl

@@ -0,0 +1,5 @@
+using DiffEqDevTools, Base.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.) ≈ 11.84 rtol=1e-2