Merge pull request #62 from JuliaReach/schillic/46

Revise code for spectral bounds

Author: schillic

src/error_bounds.jl

@@ -13,10 +13,10 @@
 # =========================================================================================================================
 
 # --- Error bounds using a priori estimate from [1] ---
-
-# See Theorem 4.2 in [1]. This is a bound based on an a priori estimate
-# of the norm of the exact solution x(t).
 # These bounds use the supremum norm (p = Inf).
+
+# see Theorem 4.2 in [1]; this bound is based on an a priori estimate
+# of the norm of the exact solution x(t)
 function error_bound_apriori(α, F₁, F₂; N)
     nF₂ = opnorm(F₂, Inf)
     μF₁ = logarithmic_norm(F₁, Inf)
@@ -26,7 +26,7 @@ function error_bound_apriori(α, F₁, F₂; N)
     return ε
 end
 
-# See Theorem 4.2 in [1]
+# see Theorem 4.2 in [1]
 function convergence_radius_apriori(α, F₁, F₂)
     nF₂ = opnorm(F₂, Inf)
     μF₁ = logarithmic_norm(F₁, Inf)
@@ -41,7 +41,7 @@ end
 
 # --- Error bounds using power series method from [1] ---
 
-# See Theorem 4.3 in [1], which uses the power series method.
+# see Theorem 4.3 in [1], which uses the power series method
 function error_bound_pseries(x₀, F₁, F₂; N)
     nx₀ = norm(x₀, Inf)
     nF₁ = opnorm(F₁, Inf)
@@ -52,7 +52,7 @@ function error_bound_pseries(x₀, F₁, F₂; N)
     return ε
 end
 
-# See Theorem 4.3 in [1].
+# see Theorem 4.3 in [1]
 function convergence_radius_pseries(x₀, F₁, F₂)
     nx₀ = norm(x₀, Inf)
     nF₁ = opnorm(F₁, Inf)
@@ -64,47 +64,37 @@ function convergence_radius_pseries(x₀, F₁, F₂)
 end
 
 # --- Error bounds using spectral abscissa from [2] ---
+# These bounds use the spectral norm (p = 2).
 
 # compute eigenvalues and sort them by increasing real part
-function _error_bound_specabs_Re_λ₁(F₁; check=true)
+function _error_bound_specabs_Re_λ₁(F₁)
     λ = eigvals(F₁; sortby=real)
-    λ₁ = last(λ)
-    Re_λ₁ = real(λ₁)
-    if check && Re_λ₁ > 0
-        throw(ArgumentError("expected Re(λ₁) ≤ 0, got $Re_λ₁"))
-    end
-    return Re_λ₁
+    return real(last(λ))
 end
 
-# See Definition (2.2) in [2]. These bounds use the spectral norm (p = 2)
-function _error_bound_specabs_R(x₀, F₂, Re_λ₁)
+# see Corollary 1 in [2]
+function error_bound_specabs(x₀, F₁, F₂; N, check=true)
+    Re_λ₁ = _error_bound_specabs_Re_λ₁(F₁)
+    if check && Re_λ₁ >= 0
+        throw(ArgumentError("expected Re(λ₁) < 0, got $Re_λ₁"))
+    end
+
+    # Equation (2.2)
     nx₀ = norm(x₀, 2)
     nF₂ = opnorm(F₂, 2)
-    R = nx₀ * nF₂ / abs(Re_λ₁)
-    return R
-end
-
-# See Corollary 1 in [2]
-function error_bound_specabs(x₀, F₁, F₂; N, check=true)
-    Re_λ₁ = _error_bound_specabs_Re_λ₁(F₁; check=check)
-    R = _error_bound_specabs_R(x₀, F₂, Re_λ₁)
+    R = nx₀ * nF₂ / -Re_λ₁
     if check && R >= 1
         throw(ArgumentError("expected R < 1, got R = $R; try scaling the ODE"))
     end
 
-    nx₀ = norm(x₀, 2)
-    if iszero(Re_λ₁)
-        nF₂ = opnorm(F₂, 2)
-        ε = t -> nx₀ * (nx₀ * nF₂ * t)^N
-    else
-        ε = t -> nx₀ * R^N * (1 - exp(Re_λ₁ * t))^N
-    end
+    # Equation (4.29)
+    ε = t -> nx₀ * R^N * (1 - exp(Re_λ₁ * t))^N
     return ε
 end
 
-# See Corollary 1 in [2]
-function convergence_radius_specabs(x₀, F₁, F₂; check=true)
-    Re_λ₁ = _error_bound_specabs_Re_λ₁(F₁; check=check)
+# see Corollary 1 in [2]
+function convergence_radius_specabs(x₀, F₁, F₂)
+    Re_λ₁ = _error_bound_specabs_Re_λ₁(F₁)
 
     if Re_λ₁ < 0
         T = Inf

test/runtests.jl

@@ -156,7 +156,6 @@ end
     # Re_λ₁ > 0
     F1_pos = [0.0 0; 0 1]
     @test_throws ArgumentError convergence_radius_specabs(x0, F1_pos, F2)
-    @test_throws ArgumentError convergence_radius_specabs(x0, F1_pos, F2; check=false)
 
     F1_neg = [-2.0 -1; -1 -2]
     F1 = [0.0 1; -1 -2]