diff --git a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl index ae95eecba..64db84621 100644 --- a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl +++ b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl @@ -17,6 +17,7 @@ include("bisection.jl") include("brent.jl") include("falsi.jl") include("itp.jl") +include("muller.jl") include("ridder.jl") # Default Algorithm @@ -44,6 +45,6 @@ end @reexport using SciMLBase, NonlinearSolveBase -export Alefeld, Bisection, Brent, Falsi, ITP, Ridder +export Alefeld, Bisection, Brent, Falsi, ITP, Muller, Ridder end diff --git a/lib/BracketingNonlinearSolve/src/muller.jl b/lib/BracketingNonlinearSolve/src/muller.jl new file mode 100644 index 000000000..edde280de --- /dev/null +++ b/lib/BracketingNonlinearSolve/src/muller.jl @@ -0,0 +1,66 @@ +""" + Muller(; middle = nothing) + +Muller's method for determining a root of a univariate, scalar function. The +algorithm, described in Sec. 9.5.2 of +[Press et al. (2007)](https://numerical.recipes/book.html), requires three +initial guesses `(left, middle, right)` for the root. + +### Keyword Arguments + +- `middle`: the initial guess for the middle point. If not provided, the + midpoint of the interval `(left, right)` is used. +""" +struct Muller{T} <: AbstractBracketingAlgorithm + middle::T +end + +Muller() = Muller(nothing) + +function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Muller, args...; + abstol = nothing, maxiters = 1000, kwargs...) + @assert !SciMLBase.isinplace(prob) "`Muller` only supports out-of-place problems." + xᵢ₋₂, xᵢ = prob.tspan + xᵢ₋₁ = isnothing(alg.middle) ? (xᵢ₋₂ + xᵢ) / 2 : alg.middle + xᵢ₋₂, xᵢ₋₁, xᵢ = promote(xᵢ₋₂, xᵢ₋₁, xᵢ) + @assert xᵢ₋₂ ≠ xᵢ₋₁ ≠ xᵢ ≠ xᵢ₋₂ + f = Base.Fix2(prob.f, prob.p) + fxᵢ₋₂, fxᵢ₋₁, fxᵢ = f(xᵢ₋₂), f(xᵢ₋₁), f(xᵢ) + + xᵢ₊₁, fxᵢ₊₁ = xᵢ₋₂, fxᵢ₋₂ + + abstol = abs(NonlinearSolveBase.get_tolerance( + xᵢ₋₂, abstol, promote_type(eltype(xᵢ₋₂), eltype(xᵢ)))) + + for _ ∈ 1:maxiters + q = (xᵢ - xᵢ₋₁)/(xᵢ₋₁ - xᵢ₋₂) + A = q*fxᵢ - q*(1 + q)*fxᵢ₋₁ + q^2*fxᵢ₋₂ + B = (2*q + 1)*fxᵢ - (1 + q)^2*fxᵢ₋₁ + q^2*fxᵢ₋₂ + C = (1 + q)*fxᵢ + + denom₊ = B + √(B^2 - 4*A*C) + denom₋ = B - √(B^2 - 4*A*C) + + if abs(denom₊) ≥ abs(denom₋) + xᵢ₊₁ = xᵢ - (xᵢ - xᵢ₋₁)*2*C/denom₊ + else + xᵢ₊₁ = xᵢ - (xᵢ - xᵢ₋₁)*2*C/denom₋ + end + + fxᵢ₊₁ = f(xᵢ₊₁) + + # Termination Check + if abstol ≥ abs(fxᵢ₊₁) + return SciMLBase.build_solution(prob, alg, xᵢ₊₁, fxᵢ₊₁; + retcode = ReturnCode.Success, + left = xᵢ₊₁, right = xᵢ₊₁) + end + + xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, xᵢ₊₁ + fxᵢ₋₂, fxᵢ₋₁, fxᵢ = fxᵢ₋₁, fxᵢ, fxᵢ₊₁ + end + + return SciMLBase.build_solution(prob, alg, xᵢ₊₁, fxᵢ₊₁; + retcode = ReturnCode.MaxIters, + left = xᵢ₊₁, right = xᵢ₊₁) +end diff --git a/lib/BracketingNonlinearSolve/test/muller_tests.jl b/lib/BracketingNonlinearSolve/test/muller_tests.jl new file mode 100644 index 000000000..f4800b9d4 --- /dev/null +++ b/lib/BracketingNonlinearSolve/test/muller_tests.jl @@ -0,0 +1,96 @@ +@testitem "Muller" begin + f(u, p) = u^2 - p + g(u, p) = sin(u) + h(u, p) = exp(-u)*sin(u) + i(u, p) = u^3 - 1 + + @testset "Quadratic function" begin + tspan = (10.0, 30.0) + p = 612.0 + prob = IntervalNonlinearProblem{false}(f, tspan, p) + sol = solve(prob, Muller()) + + @test sol.u ≈ √612 + + tspan = (-10.0, -30.0) + prob = IntervalNonlinearProblem{false}(f, tspan, p) + sol = solve(prob, Muller()) + + @test sol.u ≈ -√612 + end + + @testset "Sine function" begin + tspan = (1.0, 3.0) + prob = IntervalNonlinearProblem{false}(g, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ π + + tspan = (2.0, 6.0) + prob = IntervalNonlinearProblem{false}(g, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ 2*π + end + + @testset "Exponential-sine function" begin + tspan = (-2.0, -4.0) + prob = IntervalNonlinearProblem{false}(h, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ -π + + tspan = (-3.0, 1.0) + prob = IntervalNonlinearProblem{false}(h, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ 0 atol = 1e-15 + + tspan = (-1.0, 1.0) + prob = IntervalNonlinearProblem{false}(h, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ π + end + + @testset "Complex roots" begin + tspan = (-1.0, 1.0*im) + prob = IntervalNonlinearProblem{false}(i, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ (-1 + √3*im)/2 + + tspan = (-1.0, -1.0*im) + prob = IntervalNonlinearProblem{false}(i, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ (-1 - √3*im)/2 + end + + @testset "Middle" begin + tspan = (10.0, 30.0) + p = 612.0 + prob = IntervalNonlinearProblem{false}(f, tspan, p) + sol = solve(prob, Muller(20.0)) + + @test sol.u ≈ √612 + + tspan = (1.0, 3.0) + prob = IntervalNonlinearProblem{false}(g, tspan) + sol = solve(prob, Muller(2.0)) + + @test sol.u ≈ π + + tspan = (-2.0, -4.0) + prob = IntervalNonlinearProblem{false}(h, tspan) + sol = solve(prob, Muller(-3.0)) + + @test sol.u ≈ -π + + tspan = (-1.0, 1.0*im) + prob = IntervalNonlinearProblem{false}(i, tspan) + sol = solve(prob, Muller(0.0)) + + @test sol.u ≈ (-1 + √3*im)/2 + end +end diff --git a/lib/BracketingNonlinearSolve/test/rootfind_tests.jl b/lib/BracketingNonlinearSolve/test/rootfind_tests.jl index 9161bdfe6..e8666a8c1 100644 --- a/lib/BracketingNonlinearSolve/test/rootfind_tests.jl +++ b/lib/BracketingNonlinearSolve/test/rootfind_tests.jl @@ -7,7 +7,7 @@ end @testitem "Interval Nonlinear Problems" setup=[RootfindingTestSnippet] tags=[:core] begin using ForwardDiff - @testset for alg in (Bisection(), Falsi(), Ridder(), Brent(), ITP(), Alefeld(), nothing) + @testset for alg in (Alefeld(), Bisection(), Brent(), Falsi(), ITP(), Muller(), Ridder(), nothing) tspan = (1.0, 20.0) function g(p) @@ -52,7 +52,7 @@ end prob = IntervalNonlinearProblem(quadratic_f, (1.0, 20.0), 2.0) ϵ = eps(Float64) # least possible tol for all methods - @testset for alg in (Bisection(), Falsi(), ITP(), nothing) + @testset for alg in (Bisection(), Falsi(), ITP(), Muller(), nothing) @testset for abstol in [0.1, 0.01, 0.001, 0.0001, 1e-5, 1e-6] sol = solve(prob, alg; abstol) result_tol = abs(sol.u - sqrt(2)) @@ -62,7 +62,7 @@ end end end - @testset for alg in (Ridder(), Brent()) + @testset for alg in (Brent(), Ridder()) # Ridder and Brent converge rapidly so as we lower tolerance below 0.01, it # converges with max precision to the solution @testset for abstol in [0.1] @@ -76,7 +76,7 @@ end end @testitem "Flipped Signs and Reversed Tspan" setup=[RootfindingTestSnippet] tags=[:core] begin - @testset for alg in (Alefeld(), Bisection(), Falsi(), Brent(), ITP(), Ridder(), nothing) + @testset for alg in (Alefeld(), Bisection(), Brent(), Falsi(), ITP(), Muller(), Ridder(), nothing) f1(u, p) = u * u - p f2(u, p) = p - u * u