Modified Muller to be an IntervalNonlinearProblem

Author: fgittins

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, Muller, ITP, Ridder
 
 end

lib/SimpleNonlinearSolve/src/muller.jl lib/BracketingNonlinearSolve/src/muller.jl

@@ -1,25 +1,28 @@
 """
-    SimpleMuller()
+    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 `(xᵢ₋₂, xᵢ₋₁, xᵢ)` for the root.
 """
-struct SimpleMuller <: AbstractSimpleNonlinearSolveAlgorithm end
+struct Muller <: AbstractBracketingAlgorithm end
 
-function SciMLBase.solve(prob::NonlinearProblem, alg::SimpleMuller, args...;
-    abstol = 1e-3, 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
+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ᵢ₋₂, 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 = 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ᵢ₋₂

lib/BracketingNonlinearSolve/test/muller_tests.jl

@@ -0,0 +1,56 @@
+@testitem "Muller" begin
+    @testset "Quadratic function" begin
+        f(u, p) = u^2 - p
+
+        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
+        f(u, p) = sin(u)
+
+        tspan = (1.0, 3.0)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
+        sol = solve(prob, Muller())
+
+        @test sol.u ≈ π
+
+        tspan = (2.0, 6.0)
+        prob = IntervalNonlinearProblem{false}(f, 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)
+        sol = solve(prob, Muller())
+
+        @test sol.u ≈ -π
+
+        tspan = (-3.0, 1.0)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
+        sol = solve(prob, Muller())
+
+        @test sol.u ≈ 0 atol = 1e-15
+
+        tspan = (-1.0, 1.0)
+        prob = IntervalNonlinearProblem{false}(f, tspan)
+        sol = solve(prob, Muller())
+
+        @test sol.u ≈ π
+    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