Merge pull request #155 from DaniGlez/main

Default ITP parameters update

lib/SimpleNonlinearSolve/src/bracketing/itp.jl

@@ -15,15 +15,20 @@ I. F. D. Oliveira and R. H. C. Takahashi.
 The following keyword parameters are accepted.
 
   - `n₀::Int = 10`, the 'slack'. Must not be negative. When n₀ = 0 the worst-case is
-    identical to that of bisection, but increacing n₀ provides greater oppotunity for
+    identical to that of bisection, but increasing n₀ provides greater opportunity for
     superlinearity.
-  - `κ₁::Float64 = 0.007`. Must not be negative. The recomended value is `0.2/(x₂ - x₁)`.
+  - `scaled_κ₁::Float64 = 0.2`. Must not be negative. The recommended value is `0.2`.
     Lower values produce tighter asymptotic behaviour, while higher values improve the
     steady-state behaviour when truncation is not helpful.
-  - `κ₂::Real = 1.5`. Must lie in [1, 1+ϕ ≈ 2.62). Higher values allow for a greater
+  - `κ₂::Real = 2`. Must lie in [1, 1+ϕ ≈ 2.62). Higher values allow for a greater
     convergence rate, but also make the method more succeptable to worst-case performance.
-    In practice, κ=1,2 seems to work well due to the computational simplicity, as κ₂ is used
-    as an exponent in the method.
+    In practice, κ₂=1, 2 seems to work well due to the computational simplicity, as κ₂ is
+    used as an exponent in the method.
+
+### Computation of κ₁
+
+In the current implementation, we compute κ₁ = scaled_κ₁·|Δx₀|^(1 - κ₂); this allows κ₁ to
+adapt to the length of the interval and keep the proposed steps proportional to Δx.
 
 ### Worst Case Performance
 
@@ -35,19 +40,19 @@ n½ + `n₀` iterations, where n½ is the number of iterations using bisection
 If `f` is twice differentiable and the root is simple, then with `n₀` > 0 the convergence
 rate is √`κ₂`.
 """
-struct ITP{T} <: AbstractBracketingAlgorithm
-    k1::T
-    k2::T
+struct ITP{T₁, T₂} <: AbstractBracketingAlgorithm
+    scaled_k1::T₁
+    k2::T₂
     n0::Int
-    function ITP(; k1::Real = 0.007, k2::Real = 1.5, n0::Int = 10)
-        k1 < 0 && error("Hyper-parameter κ₁ should not be negative")
+    function ITP(;
+            scaled_k1::T₁ = 0.2, k2::T₂ = 2, n0::Int = 10) where {T₁ <: Real, T₂ <: Real}
+        scaled_k1 < 0 && error("Hyper-parameter κ₁ should not be negative")
         n0 < 0 && error("Hyper-parameter n₀ should not be negative")
         if k2 < 1 || k2 > (1.5 + sqrt(5) / 2)
             throw(ArgumentError("Hyper-parameter κ₂ should be between 1 and 1 + ϕ where \
                                  ϕ ≈ 1.618... is the golden ratio"))
         end
-        T = promote_type(eltype(k1), eltype(k2))
-        return new{T}(k1, k2, n0)
+        return new{T₁, T₂}(scaled_k1, k2, n0)
     end
 end
 
@@ -72,8 +77,8 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, args...;
     end
     ϵ = abstol
     #defining variables/cache
-    k1 = alg.k1
     k2 = alg.k2
+    k1 = alg.scaled_k1 * abs(right - left)^(1 - k2)
     n0 = alg.n0
     n_h = ceil(log2(abs(right - left) / (2 * ϵ)))
     mid = (left + right) / 2
@@ -88,7 +93,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, args...;
     while i <= maxiters
         span = abs(right - left)
         r = ϵ_s - (span / 2)
-        δ = k1 * (span^k2)
+        δ = k1 * ((k2 == 2) ? span^2 : (span^k2))
 
         ## Interpolation step ##
         x_f = left + (right - left) * (fl / (fl - fr))
@@ -119,10 +124,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, args...;
         xp <= tmin && (xp = nextfloat(tmin))
         yp = f(xp)
         yps = yp * sign(fr)
-        if yps > 0
+        T0 = zero(yps)
+        if yps > T0
             right = xp
             fr = yp
-        elseif yps < 0
+        elseif yps < T0
             left = xp
             fl = yp
         else