From 5b0ec527e2c74dd778ba31924556a4905fb18780 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Fri, 1 Sep 2023 02:03:41 +0200 Subject: [PATCH 01/19] tolerance fix bisection --- lib/SimpleNonlinearSolve/src/bisection.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index 673f77067..6e39afda9 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -20,12 +20,12 @@ function Bisection(; exact_left = false, exact_right = false) end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args...; - maxiters = 1000, + maxiters = 1000, abstol = nothing, kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - + atol = abstol !== nothing ? abstol : eps(1.0)^(3 // 5) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, From 12267cd4cb751ac612abe0e1e7c3a19f1fb16865 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Fri, 1 Sep 2023 02:28:42 +0200 Subject: [PATCH 02/19] right exact solution --- lib/SimpleNonlinearSolve/src/bisection.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index 6e39afda9..58f31f184 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -31,6 +31,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... retcode = ReturnCode.ExactSolutionLeft, left = left, right = right) end + if iszero(fr) + return SciMLBase.build_solution(prob, alg, right, fr; + retcode = ReturnCode.ExactSolutionRight, left = left, + right = right) + end i = 1 if !iszero(fr) From 0d54cc553a886ff12b6ff61e150fad995feda58f Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Tue, 5 Sep 2023 00:37:54 +0200 Subject: [PATCH 03/19] tol for Bisection() --- lib/SimpleNonlinearSolve/src/bisection.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index 58f31f184..1cb0e5d54 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -46,6 +46,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) + if abs((right - left) / 2) < atol + return SciMLBase.build_solution(prob, alg, mid, fm; + retcode = ReturnCode.Success, + left = left, right = right) + end if iszero(fm) right = mid break @@ -68,6 +73,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) + if abs((right - left) / 2) < atol + return SciMLBase.build_solution(prob, alg, mid, fm; + retcode = ReturnCode.Success, + left = left, right = right) + end if iszero(fm) right = mid fr = fm From a76441b54c2f302eaa2ec319c7981b2a61f6e57b Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Wed, 6 Sep 2023 00:09:42 +0200 Subject: [PATCH 04/19] tests for tol changed default tol to machine epsilon --- lib/SimpleNonlinearSolve/src/bisection.jl | 2 +- lib/SimpleNonlinearSolve/test/basictests.jl | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index 1cb0e5d54..17836caf4 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -25,7 +25,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : eps(1.0)^(3 // 5) + atol = abstol !== nothing ? abstol : eps(1.0) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl index 6ab9e56d8..a1e2c5669 100644 --- a/lib/SimpleNonlinearSolve/test/basictests.jl +++ b/lib/SimpleNonlinearSolve/test/basictests.jl @@ -373,6 +373,15 @@ probB = IntervalNonlinearProblem(f, tspan) sol = solve(probB, ITP()) @test sol.u ≈ sqrt(2.0) +# Tolerance tests for Interval methods + +tols = [0.1, 0.01, 0.001, 0.0001, 1e-5, 1e-6, 1e-7] + +for abstol in tols + sol = solve(probB, Bisection()) + @test abs(sol.u - sqrt(2)) < abstol +end + # Garuntee Tests for Bisection f = function (u, p) if u < 2.0 From 1dfb5dc0edf8467601668fa8996fef1e7cd79717 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Thu, 7 Sep 2023 21:01:04 +0200 Subject: [PATCH 05/19] tol fix for brent, ridder and falsi --- lib/SimpleNonlinearSolve/src/brent.jl | 13 ++++++++++++- lib/SimpleNonlinearSolve/src/falsi.jl | 17 ++++++++++++++++- lib/SimpleNonlinearSolve/src/ridder.jl | 17 ++++++++++++++++- 3 files changed, 44 insertions(+), 3 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/brent.jl b/lib/SimpleNonlinearSolve/src/brent.jl index 1cedad134..2618211bd 100644 --- a/lib/SimpleNonlinearSolve/src/brent.jl +++ b/lib/SimpleNonlinearSolve/src/brent.jl @@ -7,12 +7,13 @@ A non-allocating Brent method struct Brent <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; - maxiters = 1000, + maxiters = 1000, abstol = nothing, kwargs...) f = Base.Fix2(prob.f, prob.p) a, b = prob.tspan fa, fb = f(a), f(b) ϵ = eps(convert(typeof(fa), 1.0)) + atol = abstol !== nothing ? abstol : eps(1.0) if iszero(fa) return SciMLBase.build_solution(prob, alg, a, fa; @@ -60,6 +61,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; cond = false end fs = f(s) + if abs((b - a) / 2) < atol + return SciMLBase.build_solution(prob, alg, s, fs; + retcode = ReturnCode.Success, + left = a, right = b) + end if iszero(fs) if b < a a = b @@ -99,6 +105,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; left = a, right = b) end fc = f(c) + if abs((b - a) / 2) < atol + return SciMLBase.build_solution(prob, alg, c, fc; + retcode = ReturnCode.Success, + left = a, right = b) + end if iszero(fc) b = c fb = fc diff --git a/lib/SimpleNonlinearSolve/src/falsi.jl b/lib/SimpleNonlinearSolve/src/falsi.jl index cce11a811..b08a2a1d1 100644 --- a/lib/SimpleNonlinearSolve/src/falsi.jl +++ b/lib/SimpleNonlinearSolve/src/falsi.jl @@ -4,16 +4,21 @@ struct Falsi <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; - maxiters = 1000, + maxiters = 1000, abstol = nothing, kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) + atol = abstol !== nothing ? abstol : eps(1.0) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, right = right) + elseif iszero(fr) + return SciMLBase.build_solution(prob, alg, right, fr; + retcode = ReturnCode.ExactSolutionRight, left = left, + right = right) end i = 1 @@ -32,6 +37,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; break end fm = f(mid) + if abs((right - left) / 2) < atol + return SciMLBase.build_solution(prob, alg, mid, fm; + retcode = ReturnCode.Success, + left = left, right = right) + end if iszero(fm) right = mid break @@ -54,6 +64,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) + if abs((right - left) / 2) < atol + return SciMLBase.build_solution(prob, alg, mid, fm; + retcode = ReturnCode.Success, + left = left, right = right) + end if iszero(fm) right = mid fr = fm diff --git a/lib/SimpleNonlinearSolve/src/ridder.jl b/lib/SimpleNonlinearSolve/src/ridder.jl index 62b5a931a..099891ddd 100644 --- a/lib/SimpleNonlinearSolve/src/ridder.jl +++ b/lib/SimpleNonlinearSolve/src/ridder.jl @@ -7,16 +7,21 @@ A non-allocating ridder method struct Ridder <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; - maxiters = 1000, + maxiters = 1000, abstol = nothing, kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) + atol = abstol !== nothing ? abstol : eps(1.0) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, right = right) + elseif iszero(fr) + return SciMLBase.build_solution(prob, alg, right, fr; + retcode = ReturnCode.ExactSolutionRight, left = left, + right = right) end xo = oftype(left, Inf) @@ -37,6 +42,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; x = mid + (mid - left) * sign(fl - fr) * fm / s fx = f(x) xo = x + if abs((right - left) / 2) < atol + return SciMLBase.build_solution(prob, alg, mid, fm; + retcode = ReturnCode.Success, + left = left, right = right) + end if iszero(fx) right = x fr = fx @@ -66,6 +76,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) + if abs((right - left) / 2) < atol + return SciMLBase.build_solution(prob, alg, mid, fm; + retcode = ReturnCode.Success, + left = left, right = right) + end if iszero(fm) right = mid fr = fm From c0d76ac6296046490993e03f933aabd9e4738af4 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Thu, 7 Sep 2023 21:01:31 +0200 Subject: [PATCH 06/19] tol tests - both upper and lower bounds --- lib/SimpleNonlinearSolve/test/basictests.jl | 26 +++++++++++++++++---- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl index a1e2c5669..23253240a 100644 --- a/lib/SimpleNonlinearSolve/test/basictests.jl +++ b/lib/SimpleNonlinearSolve/test/basictests.jl @@ -374,12 +374,28 @@ sol = solve(probB, ITP()) @test sol.u ≈ sqrt(2.0) # Tolerance tests for Interval methods - +f, tspan = (u, p) -> u .* u .- 2.0, (1.0, 10.0) +probB = IntervalNonlinearProblem(f, tspan) tols = [0.1, 0.01, 0.001, 0.0001, 1e-5, 1e-6, 1e-7] - -for abstol in tols - sol = solve(probB, Bisection()) - @test abs(sol.u - sqrt(2)) < abstol +ϵ = eps(1.0) #least possible tol for all methods + +for atol in tols + sol = solve(probB, Bisection(), abstol = atol) + @test abs(sol.u - sqrt(2)) < atol + @test abs(sol.u - sqrt(2)) > ϵ #test that the solution is not calculated upto max precision + sol = solve(probB, Falsi(), abstol = atol) + @test abs(sol.u - sqrt(2)) < atol + @test abs(sol.u - sqrt(2)) > ϵ +end + +tols = [0.1] # Ridder and Brent converge rapidly so as we lower tolerance below 0.01, it converges with max precision to the solution +for atol in tols + sol = solve(probB, Ridder(), abstol = atol) + @test abs(sol.u - sqrt(2)) < atol + @test abs(sol.u - sqrt(2)) > ϵ + sol = solve(probB, Brent(), abstol = atol) + @test abs(sol.u - sqrt(2)) < atol + @test abs(sol.u - sqrt(2)) > ϵ end # Garuntee Tests for Bisection From 10bb570b40a631e9e038b40a621f1165b2915fd7 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Thu, 7 Sep 2023 21:03:36 +0200 Subject: [PATCH 07/19] right check for brent --- lib/SimpleNonlinearSolve/src/brent.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/lib/SimpleNonlinearSolve/src/brent.jl b/lib/SimpleNonlinearSolve/src/brent.jl index 2618211bd..a6706248a 100644 --- a/lib/SimpleNonlinearSolve/src/brent.jl +++ b/lib/SimpleNonlinearSolve/src/brent.jl @@ -19,6 +19,10 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; return SciMLBase.build_solution(prob, alg, a, fa; retcode = ReturnCode.ExactSolutionLeft, left = a, right = b) + elseif iszero(fb) + return SciMLBase.build_solution(prob, alg, b, fb; + retcode = ReturnCode.ExactSolutionRight, left = a, + right = b) end if abs(fa) < abs(fb) c = b From 62e0845bbe44b8a2c11215157cd7e65d7da8ce58 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:09:32 -0400 Subject: [PATCH 08/19] Update src/ridder.jl --- lib/SimpleNonlinearSolve/src/ridder.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/ridder.jl b/lib/SimpleNonlinearSolve/src/ridder.jl index 099891ddd..190a37153 100644 --- a/lib/SimpleNonlinearSolve/src/ridder.jl +++ b/lib/SimpleNonlinearSolve/src/ridder.jl @@ -12,7 +12,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : eps(1.0) + atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; From 703d984c6de1b4830e697e0f9d223252d117e18f Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:09:38 -0400 Subject: [PATCH 09/19] Update src/falsi.jl --- lib/SimpleNonlinearSolve/src/falsi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/falsi.jl b/lib/SimpleNonlinearSolve/src/falsi.jl index b08a2a1d1..549a52c13 100644 --- a/lib/SimpleNonlinearSolve/src/falsi.jl +++ b/lib/SimpleNonlinearSolve/src/falsi.jl @@ -9,7 +9,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : eps(1.0) + atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; From c378fad4c50d686951c45712af7c6a18a7cc7681 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:09:47 -0400 Subject: [PATCH 10/19] Update src/brent.jl --- lib/SimpleNonlinearSolve/src/brent.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/brent.jl b/lib/SimpleNonlinearSolve/src/brent.jl index a6706248a..1ddc0f348 100644 --- a/lib/SimpleNonlinearSolve/src/brent.jl +++ b/lib/SimpleNonlinearSolve/src/brent.jl @@ -13,7 +13,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; a, b = prob.tspan fa, fb = f(a), f(b) ϵ = eps(convert(typeof(fa), 1.0)) - atol = abstol !== nothing ? abstol : eps(1.0) + atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) if iszero(fa) return SciMLBase.build_solution(prob, alg, a, fa; From 11b7d02b0c9d4f890cce360c8dff542817536040 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:09:50 -0400 Subject: [PATCH 11/19] Update src/bisection.jl --- lib/SimpleNonlinearSolve/src/bisection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index 17836caf4..d0c134c04 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -25,7 +25,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : eps(1.0) + atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, From f5f61047cbb8cbc5e16562971d03d1c7ccdf585a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:10:34 -0400 Subject: [PATCH 12/19] Update src/bisection.jl --- lib/SimpleNonlinearSolve/src/bisection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index d0c134c04..f1b69bd7b 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -25,7 +25,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) + atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, From 0a94645e7670ac580b9a6c70a62b801e21e4472c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:11:07 -0400 Subject: [PATCH 13/19] Update src/ridder.jl --- lib/SimpleNonlinearSolve/src/ridder.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/ridder.jl b/lib/SimpleNonlinearSolve/src/ridder.jl index 190a37153..6e6e9464e 100644 --- a/lib/SimpleNonlinearSolve/src/ridder.jl +++ b/lib/SimpleNonlinearSolve/src/ridder.jl @@ -12,7 +12,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) + atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; From f1564c92973d1449a72d9db26ded8c99a9706f3a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:11:10 -0400 Subject: [PATCH 14/19] Update src/falsi.jl --- lib/SimpleNonlinearSolve/src/falsi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/falsi.jl b/lib/SimpleNonlinearSolve/src/falsi.jl index 549a52c13..52506f871 100644 --- a/lib/SimpleNonlinearSolve/src/falsi.jl +++ b/lib/SimpleNonlinearSolve/src/falsi.jl @@ -9,7 +9,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) + atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; From 4f89b5060f02ca61c55f57c9613c3c3c1decbeb5 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Sat, 9 Sep 2023 05:11:19 -0400 Subject: [PATCH 15/19] Update src/brent.jl --- lib/SimpleNonlinearSolve/src/brent.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lib/SimpleNonlinearSolve/src/brent.jl b/lib/SimpleNonlinearSolve/src/brent.jl index 1ddc0f348..11fcd5f34 100644 --- a/lib/SimpleNonlinearSolve/src/brent.jl +++ b/lib/SimpleNonlinearSolve/src/brent.jl @@ -13,7 +13,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; a, b = prob.tspan fa, fb = f(a), f(b) ϵ = eps(convert(typeof(fa), 1.0)) - atol = abstol !== nothing ? abstol : minimum(eps(left), eps(right)) + atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fa) return SciMLBase.build_solution(prob, alg, a, fa; From de98fe98a1dd57512d6748957d7d9001f9e4a7f7 Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Mon, 11 Sep 2023 01:10:58 +0200 Subject: [PATCH 16/19] tol fix for itp --- lib/SimpleNonlinearSolve/src/brent.jl | 2 +- lib/SimpleNonlinearSolve/src/itp.jl | 10 ++++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/brent.jl b/lib/SimpleNonlinearSolve/src/brent.jl index 11fcd5f34..df2cd1d7d 100644 --- a/lib/SimpleNonlinearSolve/src/brent.jl +++ b/lib/SimpleNonlinearSolve/src/brent.jl @@ -13,7 +13,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; a, b = prob.tspan fa, fb = f(a), f(b) ϵ = eps(convert(typeof(fa), 1.0)) - atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) + atol = abstol !== nothing ? abstol : min(eps(a), eps(b)) if iszero(fa) return SciMLBase.build_solution(prob, alg, a, fa; diff --git a/lib/SimpleNonlinearSolve/src/itp.jl b/lib/SimpleNonlinearSolve/src/itp.jl index fa390aa45..fa3bebe26 100644 --- a/lib/SimpleNonlinearSolve/src/itp.jl +++ b/lib/SimpleNonlinearSolve/src/itp.jl @@ -59,12 +59,12 @@ struct ITP{T} <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, - args...; abstol = 1.0e-15, + args...; abstol = nothing, maxiters = 1000, kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan # a and b fl, fr = f(left), f(right) - ϵ = abstol + ϵ = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, @@ -111,6 +111,12 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, xp = mid - (σ * r) end + if abs((left - right) / 2) < ϵ + return SciMLBase.build_solution(prob, alg, mid, f(mid); + retcode = ReturnCode.Success, + left = left, right = right) + end + ## Update ## tmin, tmax = minmax(left, right) xp >= tmax && (xp = prevfloat(tmax)) From 8a32ee67c571029ba6eeb6c106ce7c78d3e8663c Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Mon, 11 Sep 2023 01:11:07 +0200 Subject: [PATCH 17/19] test for itp --- lib/SimpleNonlinearSolve/test/basictests.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl index 23253240a..34468be58 100644 --- a/lib/SimpleNonlinearSolve/test/basictests.jl +++ b/lib/SimpleNonlinearSolve/test/basictests.jl @@ -386,6 +386,9 @@ for atol in tols sol = solve(probB, Falsi(), abstol = atol) @test abs(sol.u - sqrt(2)) < atol @test abs(sol.u - sqrt(2)) > ϵ + sol = solve(probB, ITP(), abstol = atol) + @test abs(sol.u - sqrt(2)) < atol + @test abs(sol.u - sqrt(2)) > ϵ end tols = [0.1] # Ridder and Brent converge rapidly so as we lower tolerance below 0.01, it converges with max precision to the solution From 4402777519f9c1dd6469a09d8d923493a6824efd Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Tue, 12 Sep 2023 00:06:25 +0200 Subject: [PATCH 18/19] minor arg fix --- lib/SimpleNonlinearSolve/src/bisection.jl | 8 ++++---- lib/SimpleNonlinearSolve/src/brent.jl | 7 +++---- lib/SimpleNonlinearSolve/src/falsi.jl | 7 +++---- lib/SimpleNonlinearSolve/src/itp.jl | 4 ++-- lib/SimpleNonlinearSolve/src/ridder.jl | 7 +++---- 5 files changed, 15 insertions(+), 18 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index f1b69bd7b..d583bd35d 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -20,12 +20,12 @@ function Bisection(; exact_left = false, exact_right = false) end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args...; - maxiters = 1000, abstol = nothing, + maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])), kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) + #atol = abstol if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, @@ -46,7 +46,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) - if abs((right - left) / 2) < atol + if abs((right - left) / 2) < abstol return SciMLBase.build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success, left = left, right = right) @@ -73,7 +73,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) - if abs((right - left) / 2) < atol + if abs((right - left) / 2) < abstol return SciMLBase.build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success, left = left, right = right) diff --git a/lib/SimpleNonlinearSolve/src/brent.jl b/lib/SimpleNonlinearSolve/src/brent.jl index df2cd1d7d..47e5495f0 100644 --- a/lib/SimpleNonlinearSolve/src/brent.jl +++ b/lib/SimpleNonlinearSolve/src/brent.jl @@ -7,13 +7,12 @@ A non-allocating Brent method struct Brent <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; - maxiters = 1000, abstol = nothing, + maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])), kwargs...) f = Base.Fix2(prob.f, prob.p) a, b = prob.tspan fa, fb = f(a), f(b) ϵ = eps(convert(typeof(fa), 1.0)) - atol = abstol !== nothing ? abstol : min(eps(a), eps(b)) if iszero(fa) return SciMLBase.build_solution(prob, alg, a, fa; @@ -65,7 +64,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; cond = false end fs = f(s) - if abs((b - a) / 2) < atol + if abs((b - a) / 2) < abstol return SciMLBase.build_solution(prob, alg, s, fs; retcode = ReturnCode.Success, left = a, right = b) @@ -109,7 +108,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...; left = a, right = b) end fc = f(c) - if abs((b - a) / 2) < atol + if abs((b - a) / 2) < abstol return SciMLBase.build_solution(prob, alg, c, fc; retcode = ReturnCode.Success, left = a, right = b) diff --git a/lib/SimpleNonlinearSolve/src/falsi.jl b/lib/SimpleNonlinearSolve/src/falsi.jl index 52506f871..de1079beb 100644 --- a/lib/SimpleNonlinearSolve/src/falsi.jl +++ b/lib/SimpleNonlinearSolve/src/falsi.jl @@ -4,12 +4,11 @@ struct Falsi <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; - maxiters = 1000, abstol = nothing, + maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])), kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; @@ -37,7 +36,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; break end fm = f(mid) - if abs((right - left) / 2) < atol + if abs((right - left) / 2) < abstol return SciMLBase.build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success, left = left, right = right) @@ -64,7 +63,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...; retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) - if abs((right - left) / 2) < atol + if abs((right - left) / 2) < abstol return SciMLBase.build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success, left = left, right = right) diff --git a/lib/SimpleNonlinearSolve/src/itp.jl b/lib/SimpleNonlinearSolve/src/itp.jl index fa3bebe26..c5d87bbbe 100644 --- a/lib/SimpleNonlinearSolve/src/itp.jl +++ b/lib/SimpleNonlinearSolve/src/itp.jl @@ -59,12 +59,12 @@ struct ITP{T} <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, - args...; abstol = nothing, + args...; abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])), maxiters = 1000, kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan # a and b fl, fr = f(left), f(right) - ϵ = abstol !== nothing ? abstol : min(eps(left), eps(right)) + ϵ = abstol if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, diff --git a/lib/SimpleNonlinearSolve/src/ridder.jl b/lib/SimpleNonlinearSolve/src/ridder.jl index 6e6e9464e..ce95a178a 100644 --- a/lib/SimpleNonlinearSolve/src/ridder.jl +++ b/lib/SimpleNonlinearSolve/src/ridder.jl @@ -7,12 +7,11 @@ A non-allocating ridder method struct Ridder <: AbstractBracketingAlgorithm end function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; - maxiters = 1000, abstol = nothing, + maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])), kwargs...) f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - atol = abstol !== nothing ? abstol : min(eps(left), eps(right)) if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; @@ -42,7 +41,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; x = mid + (mid - left) * sign(fl - fr) * fm / s fx = f(x) xo = x - if abs((right - left) / 2) < atol + if abs((right - left) / 2) < abstol return SciMLBase.build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success, left = left, right = right) @@ -76,7 +75,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...; retcode = ReturnCode.FloatingPointLimit, left = left, right = right) fm = f(mid) - if abs((right - left) / 2) < atol + if abs((right - left) / 2) < abstol return SciMLBase.build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success, left = left, right = right) From 35247e6daa819071fbd884cc96fa8e9e9b67e38e Mon Sep 17 00:00:00 2001 From: Yash Raj Singh Date: Tue, 12 Sep 2023 00:10:45 +0200 Subject: [PATCH 19/19] format --- lib/SimpleNonlinearSolve/src/bisection.jl | 1 - lib/SimpleNonlinearSolve/src/itp.jl | 6 +++--- 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/lib/SimpleNonlinearSolve/src/bisection.jl b/lib/SimpleNonlinearSolve/src/bisection.jl index d583bd35d..24db4adcc 100644 --- a/lib/SimpleNonlinearSolve/src/bisection.jl +++ b/lib/SimpleNonlinearSolve/src/bisection.jl @@ -25,7 +25,6 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args... f = Base.Fix2(prob.f, prob.p) left, right = prob.tspan fl, fr = f(left), f(right) - #atol = abstol if iszero(fl) return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left = left, diff --git a/lib/SimpleNonlinearSolve/src/itp.jl b/lib/SimpleNonlinearSolve/src/itp.jl index c5d87bbbe..f6688381c 100644 --- a/lib/SimpleNonlinearSolve/src/itp.jl +++ b/lib/SimpleNonlinearSolve/src/itp.jl @@ -113,10 +113,10 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, if abs((left - right) / 2) < ϵ return SciMLBase.build_solution(prob, alg, mid, f(mid); - retcode = ReturnCode.Success, - left = left, right = right) + retcode = ReturnCode.Success, + left = left, right = right) end - + ## Update ## tmin, tmax = minmax(left, right) xp >= tmax && (xp = prevfloat(tmax))