diff --git a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl b/lib/SimpleNonlinearSolve/src/bracketing/itp.jl index 2972d1c5c..4c92645ba 100644 --- a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl +++ b/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