From d1a434ec1106ef9c4dbd4d7351c609b9f8ad7bb0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Sat, 13 Jul 2024 14:16:02 +0200 Subject: [PATCH 1/4] New default values --- src/bracketing/itp.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/bracketing/itp.jl b/src/bracketing/itp.jl index 2972d1c5..4a542a95 100644 --- a/src/bracketing/itp.jl +++ b/src/bracketing/itp.jl @@ -17,13 +17,13 @@ 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 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 recomended 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. ### Worst Case Performance @@ -35,19 +35,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,7 +72,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, args...; end ϵ = abstol #defining variables/cache - k1 = alg.k1 + k1 = alg.scaled_k1 / abs(right - left) k2 = alg.k2 n0 = alg.n0 n_h = ceil(log2(abs(right - left) / (2 * ϵ))) @@ -88,7 +88,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)) From 4298697d74b5c2d5ef3ed425b1f2c2a5c4db87de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Fri, 9 Aug 2024 16:37:08 +0200 Subject: [PATCH 2/4] Docs, scaling factor --- src/bracketing/itp.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/bracketing/itp.jl b/src/bracketing/itp.jl index 4a542a95..e75642c5 100644 --- a/src/bracketing/itp.jl +++ b/src/bracketing/itp.jl @@ -15,16 +15,22 @@ 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. - - `scaled_κ₁::Float64 = 0.2`. Must not be negative. The recomended value is `0.2`. + - `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 = 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 + 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 dimension of the problem in order to keep the proposed initial step +proportional to Δx₀. + ### Worst Case Performance n½ + `n₀` iterations, where n½ is the number of iterations using bisection @@ -72,8 +78,8 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, args...; end ϵ = abstol #defining variables/cache - k1 = alg.scaled_k1 / abs(right - left) 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 From 53017a9e3e3d7c28411e5ff90453d8d789fc2f36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Fri, 9 Aug 2024 17:24:46 +0200 Subject: [PATCH 3/4] Docstring update --- src/bracketing/itp.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/bracketing/itp.jl b/src/bracketing/itp.jl index e75642c5..4cc9e9b0 100644 --- a/src/bracketing/itp.jl +++ b/src/bracketing/itp.jl @@ -28,8 +28,7 @@ The following keyword parameters are accepted. ### Computation of κ₁ In the current implementation, we compute κ₁ = scaled_κ₁·|Δx₀|^(1 - κ₂); this allows κ₁ to -adapt to the dimension of the problem in order to keep the proposed initial step -proportional to Δx₀. +adapt to the length of the interval and keep the proposed steps proportional to Δx. ### Worst Case Performance From c035945ee924b4c57320d3b2b0ec91ba0be421ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Fri, 9 Aug 2024 18:09:38 +0200 Subject: [PATCH 4/4] Fix ITP for generic output numbers --- src/bracketing/itp.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/bracketing/itp.jl b/src/bracketing/itp.jl index 4cc9e9b0..4c92645b 100644 --- a/src/bracketing/itp.jl +++ b/src/bracketing/itp.jl @@ -124,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