From 86da8b3b41af1985221ef697b4930860ace371bc 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/lib/SimpleNonlinearSolve/src/bracketing/itp.jl b/lib/SimpleNonlinearSolve/src/bracketing/itp.jl index 2972d1c5c..4a542a958 100644 --- a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl +++ b/lib/SimpleNonlinearSolve/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 2e63a1346e8209e91528ca1a6ab7dee3903965ed 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 --- lib/SimpleNonlinearSolve/src/bracketing/itp.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl b/lib/SimpleNonlinearSolve/src/bracketing/itp.jl index 4a542a958..e75642c5b 100644 --- a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl +++ b/lib/SimpleNonlinearSolve/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 ec1e309598547a20fc6006551cd416dbad78d01e 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 --- lib/SimpleNonlinearSolve/src/bracketing/itp.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl b/lib/SimpleNonlinearSolve/src/bracketing/itp.jl index e75642c5b..4cc9e9b0c 100644 --- a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl +++ b/lib/SimpleNonlinearSolve/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 145d5f221c5b4c43c24436b79390fe1b1c881620 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 --- lib/SimpleNonlinearSolve/src/bracketing/itp.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl b/lib/SimpleNonlinearSolve/src/bracketing/itp.jl index 4cc9e9b0c..4c92645ba 100644 --- a/lib/SimpleNonlinearSolve/src/bracketing/itp.jl +++ b/lib/SimpleNonlinearSolve/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