Merge branch 'master' into qqy/source

ChrisRackauckas · web-flow · commit 0df023fc8a74 · 2025-04-06T11:27:49.000-04:00
diff --git a/.github/workflows/CI_BracketingNonlinearSolve.yml b/.github/workflows/CI_BracketingNonlinearSolve.yml
@@ -31,7 +31,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_NonlinearSolve.yml b/.github/workflows/CI_NonlinearSolve.yml
@@ -44,7 +44,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_NonlinearSolveBase.yml b/.github/workflows/CI_NonlinearSolveBase.yml
@@ -30,7 +30,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_NonlinearSolveFirstOrder.yml b/.github/workflows/CI_NonlinearSolveFirstOrder.yml
@@ -31,7 +31,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_NonlinearSolveHomotopyContinuation.yml b/.github/workflows/CI_NonlinearSolveHomotopyContinuation.yml
@@ -30,7 +30,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_NonlinearSolveQuasiNewton.yml b/.github/workflows/CI_NonlinearSolveQuasiNewton.yml
@@ -31,7 +31,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_NonlinearSolveSpectralMethods.yml b/.github/workflows/CI_NonlinearSolveSpectralMethods.yml
@@ -31,7 +31,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_SCCNonlinearSolve.yml b/.github/workflows/CI_SCCNonlinearSolve.yml
@@ -31,7 +31,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_SciMLJacobianOperators.yml b/.github/workflows/CI_SciMLJacobianOperators.yml
@@ -29,7 +29,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/.github/workflows/CI_SimpleNonlinearSolve.yml b/.github/workflows/CI_SimpleNonlinearSolve.yml
@@ -32,7 +32,6 @@ jobs:
         os:
           - ubuntu-latest
           - macos-latest
-          - windows-latest
         group:
           - core
           - adjoint
diff --git a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl b/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, ITP, Muller, Ridder
 
 end
diff --git a/lib/BracketingNonlinearSolve/src/muller.jl b/lib/BracketingNonlinearSolve/src/muller.jl
@@ -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
diff --git a/lib/BracketingNonlinearSolve/test/muller_tests.jl b/lib/BracketingNonlinearSolve/test/muller_tests.jl
@@ -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
diff --git a/lib/BracketingNonlinearSolve/test/rootfind_tests.jl b/lib/BracketingNonlinearSolve/test/rootfind_tests.jl
@@ -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)
@@ -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))
@@ -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]
@@ -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
 
diff --git a/lib/SCCNonlinearSolve/Project.toml b/lib/SCCNonlinearSolve/Project.toml
@@ -18,6 +18,7 @@ ExplicitImports = "1.5"
 Hwloc = "3"
 InteractiveUtils = "<0.0.1, 1"
 NonlinearProblemLibrary = "0.1.2"
+NonlinearSolveBase = "1.5.1"
 NonlinearSolveFirstOrder = "1"
 Pkg = "1.10"
 PrecompileTools = "1.2"
@@ -36,13 +37,14 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
+NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "BenchmarkTools", "ExplicitImports", "Hwloc", "InteractiveUtils", "NonlinearSolveFirstOrder", "NonlinearProblemLibrary", "Pkg", "ReTestItems", "StableRNGs", "StaticArrays", "Test"]
+test = ["Aqua", "BenchmarkTools", "ExplicitImports", "Hwloc", "InteractiveUtils", "NonlinearProblemLibrary", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "Pkg", "ReTestItems", "StableRNGs", "StaticArrays", "Test"]
