Skip to content

Add Muller's method #568

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Apr 6, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ include("bisection.jl")
include("brent.jl")
include("falsi.jl")
include("itp.jl")
include("muller.jl")
include("ridder.jl")

# Default Algorithm
Expand Down Expand Up @@ -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
66 changes: 66 additions & 0 deletions lib/BracketingNonlinearSolve/src/muller.jl
Original file line number Diff line number Diff line change
@@ -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
96 changes: 96 additions & 0 deletions lib/BracketingNonlinearSolve/test/muller_tests.jl
Original file line number Diff line number Diff line change
@@ -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
8 changes: 4 additions & 4 deletions lib/BracketingNonlinearSolve/test/rootfind_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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))
Expand All @@ -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]
Expand All @@ -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

Expand Down
Loading