Add functionality to specify middle guess

Author: fgittins

lib/BracketingNonlinearSolve/src/muller.jl

@@ -1,18 +1,27 @@
 """
-    Muller()
+    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 `(xᵢ₋₂, xᵢ₋₁, xᵢ)` for the root.
+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 <: AbstractBracketingAlgorithm end
+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ᵢ₋₁ = (xᵢ₋₂ + xᵢ) / 2  # Use midpoint for middle guess
+    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)

lib/BracketingNonlinearSolve/test/muller_tests.jl

@@ -1,7 +1,10 @@
 @testitem "Muller" begin
-    @testset "Quadratic function" begin
-        f(u, p) = u^2 - p
+    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)
@@ -17,56 +20,77 @@
     end
 
     @testset "Sine function" begin
-        f(u, p) = sin(u)
-
         tspan = (1.0, 3.0)
-        prob = IntervalNonlinearProblem{false}(f, tspan)
+        prob = IntervalNonlinearProblem{false}(g, tspan)
         sol = solve(prob, Muller())
 
         @test sol.u ≈ π
 
         tspan = (2.0, 6.0)
-        prob = IntervalNonlinearProblem{false}(f, tspan)
+        prob = IntervalNonlinearProblem{false}(g, tspan)
         sol = solve(prob, Muller())
 
         @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}(f, tspan)
+        prob = IntervalNonlinearProblem{false}(h, tspan)
         sol = solve(prob, Muller())
 
         @test sol.u ≈ -π
 
         tspan = (-3.0, 1.0)
-        prob = IntervalNonlinearProblem{false}(f, tspan)
+        prob = IntervalNonlinearProblem{false}(h, tspan)
         sol = solve(prob, Muller())
 
         @test sol.u ≈ 0 atol = 1e-15
 
         tspan = (-1.0, 1.0)
-        prob = IntervalNonlinearProblem{false}(f, tspan)
+        prob = IntervalNonlinearProblem{false}(h, tspan)
         sol = solve(prob, Muller())
 
         @test sol.u ≈ π
     end
 
     @testset "Complex roots" begin
-        f(u, p) = u^3 - 1
-
         tspan = (-1.0, 1.0*im)
-        prob = IntervalNonlinearProblem{false}(f, tspan)
+        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}(f, tspan)
+        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