From 319e95b6ec8313554f51433fe5163da03f3292bf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Sun, 23 Jul 2023 14:32:29 +0200 Subject: [PATCH 1/4] ITP improvements --- lib/SimpleNonlinearSolve/src/itp.jl | 22 ++++++++++++--------- lib/SimpleNonlinearSolve/test/basictests.jl | 15 +++++++------- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/itp.jl b/lib/SimpleNonlinearSolve/src/itp.jl index 648208b80..b762eaea9 100644 --- a/lib/SimpleNonlinearSolve/src/itp.jl +++ b/lib/SimpleNonlinearSolve/src/itp.jl @@ -78,7 +78,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, k1 = alg.k1 k2 = alg.k2 n0 = alg.n0 - n_h = ceil(log2((right - left) / (2 * ϵ))) + n_h = ceil(log2(abs(right - left) / (2 * ϵ))) mid = (left + right) / 2 x_f = (fr * left - fl * right) / (fr - fl) xt = left @@ -89,9 +89,9 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, ϵ_s = ϵ * 2^(n_h + n0) i = 0 #iteration while i <= maxiters - #mid = (left + right) / 2 - r = ϵ_s - ((right - left) / 2) - δ = k1 * ((right - left)^k2) + span = abs(right - left) + r = ϵ_s - (span / 2) + δ = k1 * (span^k2) ## Interpolation step ## x_f = (fr * left - fl * right) / (fr - fl) @@ -112,6 +112,9 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, end ## Update ## + tmin, tmax = minmax(left, right) + xp >= tmax && (xp = prevfloat(tmax)) + xp <= tmin && (xp = nextfloat(tmin)) yp = f(xp) yps = yp * sign(fr) if yps > 0 @@ -121,16 +124,17 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, left = xp fl = yp else - left = xp - right = xp + return SciMLBase.build_solution(prob, alg, xp, yps; + retcode = ReturnCode.Success, left = left, + right = right) end i += 1 mid = (left + right) / 2 ϵ_s /= 2 - if (right - left < 2 * ϵ) - return SciMLBase.build_solution(prob, alg, mid, f(mid); - retcode = ReturnCode.Success, left = left, + if nextfloat_tdir(left, prob.tspan...) == right + return SciMLBase.build_solution(prob, alg, left, fl; + retcode = ReturnCode.FloatingPointLimit, left = left, right = right) end end diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl index ee3170ffb..678cf3a6a 100644 --- a/lib/SimpleNonlinearSolve/test/basictests.jl +++ b/lib/SimpleNonlinearSolve/test/basictests.jl @@ -540,18 +540,19 @@ for alg in (SimpleNewtonRaphson(), SimpleTrustRegion()) @test abs.(sol.u) ≈ sqrt.(p) end -# Flipped signs test +# Flipped signs & reversed tsoan test for bracketing algorithms f1(u, p) = u * u - p f2(u, p) = p - u * u -for Alg in (Alefeld, Bisection, Falsi, Brent, ITP, Ridder) - alg = Alg() +for alg in (Alefeld(), Bisection(), Falsi(), Brent(), ITP(), Ridder()) for p ∈ 1:4 inp1 = IntervalNonlinearProblem(f1, (1.0, 2.0), p) inp2 = IntervalNonlinearProblem(f2, (1.0, 2.0), p) - sol = solve(inp1, alg) - @test abs.(sol.u) ≈ sqrt.(p) - sol = solve(inp2, alg) - @test abs.(sol.u) ≈ sqrt.(p) + inp3 = IntervalNonlinearProblem(f1, (2.0, 1.0), p) + inp4 = IntervalNonlinearProblem(f2, (2.0, 1.0), p) + @test abs.(solve(inp1, alg).u) ≈ sqrt.(p) + @test abs.(solve(inp2, alg).u) ≈ sqrt.(p) + @test abs.(solve(inp3, alg).u) ≈ sqrt.(p) + @test abs.(solve(inp4, alg).u) ≈ sqrt.(p) end end \ No newline at end of file From 8e8ef86e67b9f0ce264fdfe5389cc407970834f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Sun, 23 Jul 2023 14:39:15 +0200 Subject: [PATCH 2/4] Use safer formula for falsi point in ITP --- lib/SimpleNonlinearSolve/src/itp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/itp.jl b/lib/SimpleNonlinearSolve/src/itp.jl index b762eaea9..2d8c0a69f 100644 --- a/lib/SimpleNonlinearSolve/src/itp.jl +++ b/lib/SimpleNonlinearSolve/src/itp.jl @@ -80,7 +80,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, n0 = alg.n0 n_h = ceil(log2(abs(right - left) / (2 * ϵ))) mid = (left + right) / 2 - x_f = (fr * left - fl * right) / (fr - fl) + x_f = left + (right - left) * (fl/(fl - fr)) xt = left xp = left r = zero(left) #minmax radius @@ -94,7 +94,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, δ = k1 * (span^k2) ## Interpolation step ## - x_f = (fr * left - fl * right) / (fr - fl) + x_f = left + (right - left) * (fl/(fl - fr)) ## Truncation step ## σ = sign(mid - x_f) From 65f24fe17c1e3bb8106fb06c6c6bc5935165c71e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Sun, 23 Jul 2023 14:45:08 +0200 Subject: [PATCH 3/4] =?UTF-8?q?Refine=20ITP=20return=20solution=20on=20hit?= =?UTF-8?q?ting=20f(=C2=B7)=20zero?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- lib/SimpleNonlinearSolve/src/itp.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/itp.jl b/lib/SimpleNonlinearSolve/src/itp.jl index 2d8c0a69f..1f1efbafe 100644 --- a/lib/SimpleNonlinearSolve/src/itp.jl +++ b/lib/SimpleNonlinearSolve/src/itp.jl @@ -125,8 +125,8 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, fl = yp else return SciMLBase.build_solution(prob, alg, xp, yps; - retcode = ReturnCode.Success, left = left, - right = right) + retcode = ReturnCode.Success, left = xp, + right = xp) end i += 1 mid = (left + right) / 2 From fbd4f68614d02689c080a2af173ed863da660a87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Gonz=C3=A1lez?= Date: Sun, 23 Jul 2023 14:45:45 +0200 Subject: [PATCH 4/4] typo --- lib/SimpleNonlinearSolve/test/basictests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl index 95f8d58e2..6ab9e56d8 100644 --- a/lib/SimpleNonlinearSolve/test/basictests.jl +++ b/lib/SimpleNonlinearSolve/test/basictests.jl @@ -540,7 +540,7 @@ for alg in (SimpleNewtonRaphson(), SimpleTrustRegion()) @test abs.(sol.u) ≈ sqrt.(p) end -# Flipped signs & reversed tsoan test for bracketing algorithms +# Flipped signs & reversed tspan test for bracketing algorithms f1(u, p) = u * u - p f2(u, p) = p - u * u