Add Muller's method

src/SimpleNonlinearSolve.jl

@@ -37,6 +37,7 @@ include("nlsolve/klement.jl")
 include("nlsolve/trustRegion.jl")
 include("nlsolve/halley.jl")
 include("nlsolve/dfsane.jl")
+include("nlsolve/muller.jl")
 
 ## Interval Nonlinear Solvers
 include("bracketing/bisection.jl")
@@ -115,7 +116,7 @@ end
 end
 
 export SimpleBroyden, SimpleDFSane, SimpleGaussNewton, SimpleHalley, SimpleKlement,
-       SimpleLimitedMemoryBroyden, SimpleNewtonRaphson, SimpleTrustRegion
+       SimpleLimitedMemoryBroyden, SimpleNewtonRaphson, SimpleTrustRegion, SimpleMuller
 export Alefeld, Bisection, Brent, Falsi, ITP, Ridder
 
 end # module

src/nlsolve/muller.jl

@@ -0,0 +1,55 @@
+"""
+    SimpleMuller()
+
+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 `(xᵢ₋₂, xᵢ₋₁, xᵢ)` for the root.
+"""
+struct SimpleMuller <: AbstractSimpleNonlinearSolveAlgorithm end
+
+function SciMLBase.solve(prob::NonlinearProblem, alg::SimpleMuller, args...;
+    abstol = nothing, maxiters = 1000, kwargs...)
+    @assert !isinplace(prob) "`SimpleMuller` only supports OOP problems."
+    @assert length(prob.u0) == 3 "`SimpleMuller` requires three initial guesses."
+    xᵢ₋₂, xᵢ₋₁, xᵢ = prob.u0
+    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ᵢ)
+
+    abstol = __get_tolerance(nothing, abstol,
+            promote_type(eltype(fxᵢ₋₂), eltype(xᵢ₋₂)))
+
+    xᵢ₊₁, fxᵢ₊₁ = xᵢ₋₂, fxᵢ₋₂
+
+    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 build_solution(prob, alg, xᵢ₊₁, fxᵢ₊₁;
+                                  retcode = ReturnCode.Success)
+        end
+
+        xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, xᵢ₊₁
+        fxᵢ₋₂, fxᵢ₋₁, fxᵢ = fxᵢ₋₁, fxᵢ, fxᵢ₊₁
+    end
+
+    return build_solution(prob, alg, xᵢ₊₁, fxᵢ₊₁;
+                          retcode = ReturnCode.MaxIters)
+end

test/core/muller_tests.jl

@@ -0,0 +1,56 @@
+@testitem "SimpleMuller" begin
+    @testset "Quadratic function" begin
+        f(u, p) = u^2 - p
+
+        u0 = (10.0, 20.0, 30.0)
+        p = 612.0
+        prob = NonlinearProblem{false}(f, u0, p)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ √612
+
+        u0 = (-10.0, -20.0, -30.0)
+        prob = NonlinearProblem{false}(f, u0, p)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ -√612
+    end
+
+    @testset "Sine function" begin
+        f(u, p) = sin(u)
+
+        u0 = (1.0, 2.0, 3.0)
+        prob = NonlinearProblem{false}(f, u0)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ π
+
+        u0 = (2.0, 4.0, 6.0)
+        prob = NonlinearProblem{false}(f, u0)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ 2*π
+    end
+
+    @testset "Exponential-sine function" begin
+        f(u, p) = exp(-u)*sin(u)
+
+        u0 = (-2.0, -3.0, -4.0)
+        prob = NonlinearProblem{false}(f, u0)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ -π
+
+        u0 = (-1.0, 0.0, 1/2)
+        prob = NonlinearProblem{false}(f, u0)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ 0
+
+        u0 = (-1.0, 0.0, 1.0)
+        prob = NonlinearProblem{false}(f, u0)
+        sol = solve(prob, SimpleMuller())
+
+        @test sol.u ≈ π
+    end
+end