diff --git a/lib/BracketingNonlinearSolve/Project.toml b/lib/BracketingNonlinearSolve/Project.toml index 682873daa..e4b92fd1d 100644 --- a/lib/BracketingNonlinearSolve/Project.toml +++ b/lib/BracketingNonlinearSolve/Project.toml @@ -1,7 +1,7 @@ name = "BracketingNonlinearSolve" uuid = "70df07ce-3d50-431d-a3e7-ca6ddb60ac1e" authors = ["Avik Pal and contributors"] -version = "1.6.1" +version = "1.6.2" [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" diff --git a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl index fd9398832..6063d21fb 100644 --- a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl +++ b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl @@ -10,6 +10,7 @@ using SciMLBase: SciMLBase, IntervalNonlinearProblem, ReturnCode abstract type AbstractBracketingAlgorithm <: AbstractNonlinearSolveAlgorithm end +include("utils.jl") include("common.jl") include("alefeld.jl") diff --git a/lib/BracketingNonlinearSolve/src/alefeld.jl b/lib/BracketingNonlinearSolve/src/alefeld.jl index 86807986a..06aef3fb8 100644 --- a/lib/BracketingNonlinearSolve/src/alefeld.jl +++ b/lib/BracketingNonlinearSolve/src/alefeld.jl @@ -18,15 +18,11 @@ function SciMLBase.__solve( fc = f(c) if a == c || b == c - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit, left = a, right = b - ) + return build_bracketing_solution(prob, alg, c, fc, a, b, ReturnCode.FloatingPointLimit) end if iszero(fc) - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.Success, left = a, right = b - ) + return build_exact_solution(prob, alg, c, fc, ReturnCode.Success) end a, b, d = Impl.bracket(f, a, b, c) @@ -46,14 +42,11 @@ function SciMLBase.__solve( ē, fc = d, f(c) if a == c || b == c - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit, - left = a, right = b) + return build_bracketing_solution(prob, alg, c, fc, a, b, ReturnCode.FloatingPointLimit) end if iszero(fc) - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.Success, left = a, right = b) + return build_exact_solution(prob, alg, c, fc, ReturnCode.Success) end ā, b̄, d̄ = Impl.bracket(f, a, b, c) @@ -71,16 +64,11 @@ function SciMLBase.__solve( fc = f(c) if ā == c || b̄ == c - return SciMLBase.build_solution( - prob, alg, c, fc; - retcode = ReturnCode.FloatingPointLimit, left = ā, right = b̄ - ) + return build_bracketing_solution(prob, alg, c, fc, ā, b̄, ReturnCode.FloatingPointLimit) end if iszero(fc) - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.Success, left = ā, right = b̄ - ) + return build_exact_solution(prob, alg, c, fc, ReturnCode.Success) end ā, b̄, d̄ = Impl.bracket(f, ā, b̄, c) @@ -94,16 +82,11 @@ function SciMLBase.__solve( fc = f(c) if ā == c || b̄ == c - return SciMLBase.build_solution( - prob, alg, c, fc; - retcode = ReturnCode.FloatingPointLimit, left = ā, right = b̄ - ) + return build_bracketing_solution(prob, alg, c, fc, ā, b̄, ReturnCode.FloatingPointLimit) end if iszero(fc) - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.Success, left = ā, right = b̄ - ) + return build_exact_solution(prob, alg, c, fc, ReturnCode.Success) end ā, b̄, d = Impl.bracket(f, ā, b̄, c) @@ -117,15 +100,10 @@ function SciMLBase.__solve( fc = f(c) if ā == c || b̄ == c - return SciMLBase.build_solution( - prob, alg, c, fc; - retcode = ReturnCode.FloatingPointLimit, left = ā, right = b̄ - ) + return build_bracketing_solution(prob, alg, c, fc, ā, b̄, ReturnCode.FloatingPointLimit) end if iszero(fc) - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.Success, left = ā, right = b̄ - ) + return build_exact_solution(prob, alg, c, fc, ReturnCode.Success) end a, b, d = Impl.bracket(f, ā, b̄, c) end @@ -140,7 +118,5 @@ function SciMLBase.__solve( fc = f(c) # Return solution when run out of max iteration - return SciMLBase.build_solution( - prob, alg, c, fc; retcode = ReturnCode.MaxIters, left = a, right = b - ) + return build_bracketing_solution(prob, alg, c, fc, a, b, ReturnCode.MaxIters) end diff --git a/lib/BracketingNonlinearSolve/src/bisection.jl b/lib/BracketingNonlinearSolve/src/bisection.jl index 153772d27..e3f74e662 100644 --- a/lib/BracketingNonlinearSolve/src/bisection.jl +++ b/lib/BracketingNonlinearSolve/src/bisection.jl @@ -33,46 +33,41 @@ function SciMLBase.__solve( left, abstol, promote_type(eltype(left), eltype(right))) if iszero(fl) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right - ) + return build_exact_solution(prob, alg, left, fl, ReturnCode.ExactSolutionLeft) end if iszero(fr) - return SciMLBase.build_solution( - prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right - ) + return build_exact_solution(prob, alg, right, fr, ReturnCode.ExactSolutionRight) end if sign(fl) == sign(fr) @SciMLMessage("The interval is not an enclosing interval, opposite signs at the \ boundaries are required.", verbose, :non_enclosing_interval) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.InitialFailure) end + return internal_bisection(f, left, right, fl, fr, abstol, maxiters, prob, alg) +end + +# Bisection main loop is implemented in separate function so that it can be reused +# as a fallback solver in other solvers +function internal_bisection(f::F, left, right, fl, fr, abstol, maxiters, prob, alg) where {F} i = 1 while i ≤ maxiters mid = (left + right) / 2 if mid == left || mid == right - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.FloatingPointLimit, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.FloatingPointLimit) end fm = f(mid) - if abs((right - left) / 2) < abstol - return SciMLBase.build_solution( - prob, alg, mid, fm; retcode = ReturnCode.Success, left, right - ) + if iszero(fm) + return build_exact_solution(prob, alg, mid, fm, ReturnCode.Success) end - if iszero(fm) - right = mid - break + if abs((right - left) / 2) < abstol + return build_bracketing_solution(prob, alg, mid, fm, left, right, ReturnCode.Success) end if sign(fl) == sign(fm) @@ -86,14 +81,5 @@ function SciMLBase.__solve( i += 1 end - sol, i, left, right, - fl, fr = Impl.bisection( - left, right, fl, fr, f, abstol, maxiters - i, prob, alg - ) - - sol !== nothing && return sol - - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right - ) -end + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.MaxIters) +end \ No newline at end of file diff --git a/lib/BracketingNonlinearSolve/src/brent.jl b/lib/BracketingNonlinearSolve/src/brent.jl index cce3428ad..e07b09903 100644 --- a/lib/BracketingNonlinearSolve/src/brent.jl +++ b/lib/BracketingNonlinearSolve/src/brent.jl @@ -31,24 +31,18 @@ function SciMLBase.__solve( ) if iszero(fl) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right - ) + return build_exact_solution(prob, alg, left, fl, ReturnCode.ExactSolutionLeft) end if iszero(fr) - return SciMLBase.build_solution( - prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right - ) + return build_exact_solution(prob, alg, right, fr, ReturnCode.ExactSolutionRight) end if sign(fl) == sign(fr) @SciMLMessage("The interval is not an enclosing interval, opposite signs at the \ boundaries are required.", verbose, :non_enclosing_interval) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.InitialFailure) end if abs(fl) < abs(fr) @@ -83,10 +77,7 @@ function SciMLBase.__solve( # Bisection method s = (left + right) / 2 if s == left || s == right - return SciMLBase.build_solution( - prob, alg, left, fl; - retcode = ReturnCode.FloatingPointLimit, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.FloatingPointLimit) end cond = true else @@ -94,20 +85,12 @@ function SciMLBase.__solve( end fs = f(s) - if abs((right - left) / 2) < abstol - return SciMLBase.build_solution( - prob, alg, s, fs; retcode = ReturnCode.Success, left, right - ) + if iszero(fs) + return build_exact_solution(prob, alg, s, fs, ReturnCode.Success) end - if iszero(fs) - if right < left - left = right - fl = fr - end - right = s - fr = fs - break + if abs((right - left) / 2) < abstol + return build_bracketing_solution(prob, alg, s, fs, left, right, ReturnCode.Success) end if fl * fs < 0 @@ -128,14 +111,5 @@ function SciMLBase.__solve( i += 1 end - sol, i, left, right, - fl, fr = Impl.bisection( - left, right, fl, fr, f, abstol, maxiters - i, prob, alg - ) - - sol !== nothing && return sol - - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.MaxIters) end diff --git a/lib/BracketingNonlinearSolve/src/common.jl b/lib/BracketingNonlinearSolve/src/common.jl index 0127869d8..d7728155a 100644 --- a/lib/BracketingNonlinearSolve/src/common.jl +++ b/lib/BracketingNonlinearSolve/src/common.jl @@ -1,42 +1,5 @@ module Impl -using SciMLBase: SciMLBase, ReturnCode - -function bisection(left, right, fl, fr, f::F, abstol, maxiters, prob, alg) where {F} - i = 1 - sol = nothing - while i ≤ maxiters - mid = (left + right) / 2 - - if mid == left || mid == right - sol = SciMLBase.build_solution( - prob, alg, left, fl; left, right, retcode = ReturnCode.FloatingPointLimit - ) - break - end - - fm = f(mid) - if abs((right - left) / 2) < abstol - sol = SciMLBase.build_solution( - prob, alg, mid, fm; left, right, retcode = ReturnCode.Success - ) - break - end - - if iszero(fm) - right = mid - fr = fm - else - left = mid - fl = fm - end - - i += 1 - end - - return sol, i, left, right, fl, fr -end - prevfloat_tdir(x, x0, x1) = ifelse(x1 > x0, prevfloat(x), nextfloat(x)) nextfloat_tdir(x, x0, x1) = ifelse(x1 > x0, nextfloat(x), prevfloat(x)) max_tdir(a, b, x0, x1) = ifelse(x1 > x0, max(a, b), min(a, b)) diff --git a/lib/BracketingNonlinearSolve/src/falsi.jl b/lib/BracketingNonlinearSolve/src/falsi.jl index edd57cbc6..f0d09c95f 100644 --- a/lib/BracketingNonlinearSolve/src/falsi.jl +++ b/lib/BracketingNonlinearSolve/src/falsi.jl @@ -30,32 +30,24 @@ function SciMLBase.__solve( left, abstol, promote_type(eltype(left), eltype(right))) if iszero(fl) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right - ) + return build_exact_solution(prob, alg, left, fl, ReturnCode.ExactSolutionLeft) end if iszero(fr) - return SciMLBase.build_solution( - prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right - ) + return build_exact_solution(prob, alg, right, fr, ReturnCode.ExactSolutionRight) end if sign(fl) == sign(fr) @SciMLMessage("The interval is not an enclosing interval, opposite signs at the \ boundaries are required.", verbose, :non_enclosing_interval) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.InitialFailure) end i = 1 while i ≤ maxiters if Impl.nextfloat_tdir(left, l, r) == right - return SciMLBase.build_solution( - prob, alg, left, fl; left, right, retcode = ReturnCode.FloatingPointLimit - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.FloatingPointLimit) end mid = (fr * left - fl * right) / (fr - fl) @@ -66,10 +58,12 @@ function SciMLBase.__solve( (mid == left || mid == right) && break fm = f(mid) + if iszero(fm) + return build_exact_solution(prob, alg, mid, fm, ReturnCode.Success) + end + if abs((right - left) / 2) < abstol - return SciMLBase.build_solution( - prob, alg, mid, fm; left, right, retcode = ReturnCode.Success - ) + return build_bracketing_solution(prob, alg, mid, fm, left, right, ReturnCode.Success) end if abs(fm) < abstol @@ -86,14 +80,6 @@ function SciMLBase.__solve( i += 1 end - sol, i, left, right, - fl, fr = Impl.bisection( - left, right, fl, fr, f, abstol, maxiters - i, prob, alg - ) - - sol !== nothing && return sol - - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right - ) + # Fallback to bisection solver + return internal_bisection(f, left, right, fl, fr, abstol, maxiters - i, prob, alg) end diff --git a/lib/BracketingNonlinearSolve/src/itp.jl b/lib/BracketingNonlinearSolve/src/itp.jl index 2eb840749..b39b25675 100644 --- a/lib/BracketingNonlinearSolve/src/itp.jl +++ b/lib/BracketingNonlinearSolve/src/itp.jl @@ -71,24 +71,18 @@ function SciMLBase.__solve( ) if iszero(fl) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right - ) + return build_exact_solution(prob, alg, left, fl, ReturnCode.ExactSolutionLeft) end if iszero(fr) - return SciMLBase.build_solution( - prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right - ) + return build_exact_solution(prob, alg, right, fr, ReturnCode.ExactSolutionRight) end if sign(fl) == sign(fr) @SciMLMessage("The interval is not an enclosing interval, opposite signs at the \ boundaries are required.", verbose, :non_enclosing_interval) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.InitialFailure) end ϵ = abstol @@ -120,9 +114,7 @@ function SciMLBase.__solve( xp = ifelse(abs(xt - mid) ≤ r, xt, mid - copysign(r, diff)) # Projection Step if span < 2ϵ - return SciMLBase.build_solution( - prob, alg, xt, f(xt); retcode = ReturnCode.Success, left, right - ) + return build_bracketing_solution(prob, alg, xt, f(xt), left, right, ReturnCode.Success) end yp = f(xp) yps = yp * sign(fr) @@ -131,22 +123,16 @@ function SciMLBase.__solve( elseif yps < T0 left, fl = xp, yp else - return SciMLBase.build_solution( - prob, alg, xp, yps; retcode = ReturnCode.Success, left, right - ) + return build_exact_solution(prob, alg, xp, yps, ReturnCode.Success) end i += 1 ϵ_s /= 2 if nextfloat(left) == right - return SciMLBase.build_solution( - prob, alg, right, fr; retcode = ReturnCode.FloatingPointLimit, left, right - ) + return build_bracketing_solution(prob, alg, right, fr, left, right, ReturnCode.FloatingPointLimit) end end - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.MaxIters) end diff --git a/lib/BracketingNonlinearSolve/src/ridder.jl b/lib/BracketingNonlinearSolve/src/ridder.jl index e0686be95..2785ab35b 100644 --- a/lib/BracketingNonlinearSolve/src/ridder.jl +++ b/lib/BracketingNonlinearSolve/src/ridder.jl @@ -20,57 +20,46 @@ function SciMLBase.__solve( ) if iszero(fl) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft, left, right - ) + return build_exact_solution(prob, alg, left, fl, ReturnCode.ExactSolutionLeft) end if iszero(fr) - return SciMLBase.build_solution( - prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight, left, right - ) + return build_exact_solution(prob, alg, right, fr, ReturnCode.ExactSolutionRight) end if sign(fl) == sign(fr) @SciMLMessage("The interval is not an enclosing interval, opposite signs at the \ boundaries are required.", verbose, :non_enclosing_interval) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.InitialFailure, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.InitialFailure) end - xo = oftype(left, Inf) i = 1 while i ≤ maxiters mid = (left + right) / 2 if mid == left || mid == right - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.FloatingPointLimit, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.FloatingPointLimit) end fm = f(mid) + if iszero(fm) + return build_exact_solution(prob, alg, mid, fm, ReturnCode.Success) + end + s = sqrt(fm^2 - fl * fr) if iszero(s) - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.Failure, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.Failure) end x = mid + (mid - left) * sign(fl - fm) * fm / s fx = f(x) - xo = x - if abs((right - left) / 2) < abstol - return SciMLBase.build_solution( - prob, alg, mid, fm; retcode = ReturnCode.Success, left, right - ) + if iszero(fx) + return build_exact_solution(prob, alg, x, fx, ReturnCode.Success) end - if iszero(fx) - right, fr = x, fx - break + if abs((right - left) / 2) < abstol + return build_bracketing_solution(prob, alg, mid, fm, left, right, ReturnCode.Success) end if sign(fx) != sign(fm) @@ -90,14 +79,5 @@ function SciMLBase.__solve( i += 1 end - sol, i, left, right, - fl, fr = Impl.bisection( - left, right, fl, fr, f, abstol, maxiters - i, prob, alg - ) - - sol !== nothing && return sol - - return SciMLBase.build_solution( - prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right - ) + return build_bracketing_solution(prob, alg, left, fl, left, right, ReturnCode.MaxIters) end diff --git a/lib/BracketingNonlinearSolve/src/utils.jl b/lib/BracketingNonlinearSolve/src/utils.jl new file mode 100644 index 000000000..637dd0e3e --- /dev/null +++ b/lib/BracketingNonlinearSolve/src/utils.jl @@ -0,0 +1,24 @@ +""" +Builder for exact bracketing problems solution. +""" +function build_exact_solution(prob, alg, u, resid, retcode) + SciMLBase.build_solution( + prob, alg, u, resid; retcode = retcode, left = u, right = u + ) +end + +""" +Builder for bracketing problems solution. +Ensure that left/right are in the same order as tspan for consistency. +""" +function build_bracketing_solution(prob, alg, u, resid, bound1, bound2, retcode) + if xor(bound1 < bound2, prob.tspan[1] < prob.tspan[2]) + SciMLBase.build_solution( + prob, alg, u, resid; retcode = retcode, left = bound2, right = bound1 + ) + else + SciMLBase.build_solution( + prob, alg, u, resid; retcode = retcode, left = bound1, right = bound2 + ) + end +end \ No newline at end of file diff --git a/lib/BracketingNonlinearSolve/test/rootfind_tests.jl b/lib/BracketingNonlinearSolve/test/rootfind_tests.jl index 5b368bea2..633f7b2f8 100644 --- a/lib/BracketingNonlinearSolve/test/rootfind_tests.jl +++ b/lib/BracketingNonlinearSolve/test/rootfind_tests.jl @@ -1,4 +1,5 @@ @testsnippet RootfindingTestSnippet begin + linear_f(u, p) = u .- p quadratic_f(u, p) = u .* u .- p quadratic_f!(du, u, p) = (du .= u .* u .- p) quadratic_f2(u, p) = @. p[1] * u * u - p[2] @@ -74,15 +75,27 @@ end @test result_tol > ϵ end end +end + +@testitem "Zero-Tolerance Tests Interval Methods" setup=[RootfindingTestSnippet] tags=[:core] begin + prob=IntervalNonlinearProblem(quadratic_f, (1.0, 20.0), 2.0) + prob_lin=IntervalNonlinearProblem(linear_f, (-1.0, 1.0), 0.0) - # Some solvers support abstol=0.0 and converge to floating point precision - @testset for alg in (Bisection(), ITP(), Brent()) + @testset for alg in (Alefeld(), Bisection(), Brent(), ITP(), Ridder(), nothing) sol = solve(prob, alg; abstol = 0.0) # Test that solution is to floating point precision @test sol.retcode == ReturnCode.FloatingPointLimit @test quadratic_f(sol.left, 2.0) < 0 @test quadratic_f(sol.right, 2.0) > 0 @test nextfloat(sol.left) == sol.right + + # Solve problem with a root representable with floating point + sol = solve(prob_lin, alg; abstol = 0.0) + # Test that solution is exact + @test sol.retcode == ReturnCode.Success + @test sol.u == 0.0 + @test sol.left == 0.0 + @test sol.right == 0.0 end end @@ -93,14 +106,19 @@ end f2(u, p) = p - u * u for p in 1:4 - inp1 = IntervalNonlinearProblem(f1, (1.0, 2.0), p) - inp2 = IntervalNonlinearProblem(f2, (1.0, 2.0), 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) + sol1 = solve(IntervalNonlinearProblem(f1, (1.0, 2.0), p), alg) + sol2 = solve(IntervalNonlinearProblem(f2, (1.0, 2.0), p), alg) + sol3 = solve(IntervalNonlinearProblem(f1, (2.0, 1.0), p), alg) + sol4 = solve(IntervalNonlinearProblem(f2, (2.0, 1.0), p), alg) + @test abs.(sol1.u) ≈ sqrt.(p) + @test abs.(sol2.u) ≈ sqrt.(p) + @test abs.(sol3.u) ≈ sqrt.(p) + @test abs.(sol4.u) ≈ sqrt.(p) + # Test brackets consistency + @test sol1.left ≤ sol1.right + @test sol2.left ≤ sol2.right + @test sol3.left ≥ sol3.right + @test sol4.left ≥ sol4.right end end end