Modify Muller to be an IntervalNonlinearProblem

Author: fgittins

lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl

@@ -45,6 +45,6 @@ end
 
 @reexport using SciMLBase, NonlinearSolveBase
 
-export Alefeld, Bisection, Brent, Falsi, ITP, Muller, Ridder
+export Alefeld, Bisection, Brent, Falsi, Muller, ITP, Ridder
 
 end

lib/BracketingNonlinearSolve/src/muller.jl

@@ -1,66 +1,55 @@
 """
-    Muller(; middle = nothing)
+    Muller()
 
 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.
+initial guesses `(xᵢ₋₂, xᵢ₋₁, xᵢ)` for the root.
 """
-struct Muller{T} <: AbstractBracketingAlgorithm
-    middle::T
-end
-
-Muller() = Muller(nothing)
+struct Muller <: AbstractBracketingAlgorithm end
 
 function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Muller, args...;
-        abstol = nothing, maxiters = 1000, kwargs...)
+    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ᵢ) / 2  # Use midpoint for middle guess
     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ᵢ))))
+    abstol = NonlinearSolveBase.get_tolerance(
+        xᵢ₋₂, abstol, promote_type(eltype(xᵢ₋₂), eltype(xᵢ)))
 
-    for _ in 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ᵢ
+    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)
+        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₊
+            xᵢ₊₁ = xᵢ - (xᵢ - xᵢ₋₁)*2*C/denom₊
         else
-            xᵢ₊₁ = xᵢ - (xᵢ - xᵢ₋₁) * 2 * C / denom₋
+            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ᵢ₊₁)
+                                            retcode = ReturnCode.Success)
         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ᵢ₊₁)
+                                    retcode = ReturnCode.MaxIters)
 end

lib/BracketingNonlinearSolve/test/muller_tests.jl

@@ -1,10 +1,7 @@
 @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
+        f(u, p) = u^2 - p
+
         tspan = (10.0, 30.0)
         p = 612.0
         prob = IntervalNonlinearProblem{false}(f, tspan, p)
@@ -20,77 +17,40 @@
     end
 
     @testset "Sine function" begin
+        f(u, p) = sin(u)
+
         tspan = (1.0, 3.0)
-        prob = IntervalNonlinearProblem{false}(g, tspan)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
         sol = solve(prob, Muller())
 
         @test sol.u ≈ π
 
         tspan = (2.0, 6.0)
-        prob = IntervalNonlinearProblem{false}(g, tspan)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
         sol = solve(prob, Muller())
 
-        @test sol.u ≈ 2 * π
+        @test sol.u ≈ 2*π
     end
 
     @testset "Exponential-sine function" begin
+        f(u, p) = exp(-u)*sin(u)
+
         tspan = (-2.0, -4.0)
-        prob = IntervalNonlinearProblem{false}(h, tspan)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
         sol = solve(prob, Muller())
 
         @test sol.u ≈ -π
 
         tspan = (-3.0, 1.0)
-        prob = IntervalNonlinearProblem{false}(h, tspan)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
         sol = solve(prob, Muller())
 
-        @test sol.u≈0 atol=1e-15
+        @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)
+        prob = IntervalNonlinearProblem{false}(f, 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

lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl

@@ -54,7 +54,6 @@ include("klement.jl")
 include("lbroyden.jl")
 include("raphson.jl")
 include("trust_region.jl")
-include("muller.jl")
 
 # By Pass the highlevel checks for NonlinearProblem for Simple Algorithms
 function CommonSolve.solve(
@@ -167,7 +166,6 @@ export SimpleBroyden, SimpleKlement, SimpleLimitedMemoryBroyden
 export SimpleDFSane
 export SimpleGaussNewton, SimpleNewtonRaphson, SimpleTrustRegion
 export SimpleHalley
-export SimpleMuller
 
 export solve