From 69ba15ce18359968dd1b9f1770238b3c3d2bdbd6 Mon Sep 17 00:00:00 2001 From: Fabian Gittins Date: Sat, 5 Apr 2025 12:37:28 +0200 Subject: [PATCH 1/5] Add Muller's method --- .../src/SimpleNonlinearSolve.jl | 2 + lib/SimpleNonlinearSolve/src/muller.jl | 52 +++++++++++++++++ .../test/core/muller_tests.jl | 56 +++++++++++++++++++ 3 files changed, 110 insertions(+) create mode 100644 lib/SimpleNonlinearSolve/src/muller.jl create mode 100644 lib/SimpleNonlinearSolve/test/core/muller_tests.jl diff --git a/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl b/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl index 495cadef2..1e284502f 100644 --- a/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl +++ b/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl @@ -54,6 +54,7 @@ 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( @@ -166,6 +167,7 @@ export SimpleBroyden, SimpleKlement, SimpleLimitedMemoryBroyden export SimpleDFSane export SimpleGaussNewton, SimpleNewtonRaphson, SimpleTrustRegion export SimpleHalley +export SimpleMuller export solve diff --git a/lib/SimpleNonlinearSolve/src/muller.jl b/lib/SimpleNonlinearSolve/src/muller.jl new file mode 100644 index 000000000..ed8874e8e --- /dev/null +++ b/lib/SimpleNonlinearSolve/src/muller.jl @@ -0,0 +1,52 @@ +""" + 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 = 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 + 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ᵢ₋₂ + + 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) + 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) +end diff --git a/lib/SimpleNonlinearSolve/test/core/muller_tests.jl b/lib/SimpleNonlinearSolve/test/core/muller_tests.jl new file mode 100644 index 000000000..fab2fbd0e --- /dev/null +++ b/lib/SimpleNonlinearSolve/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 From fb9d5a4a3639e532ce9ab65b74071353157da270 Mon Sep 17 00:00:00 2001 From: Fabian Gittins Date: Sat, 5 Apr 2025 13:14:06 +0200 Subject: [PATCH 2/5] Modify Muller to be an IntervalNonlinearProblem --- .../src/BracketingNonlinearSolve.jl | 3 +- .../src/muller.jl | 17 +++--- .../test/muller_tests.jl | 56 +++++++++++++++++++ .../src/SimpleNonlinearSolve.jl | 2 - .../test/core/muller_tests.jl | 56 ------------------- 5 files changed, 68 insertions(+), 66 deletions(-) rename lib/{SimpleNonlinearSolve => BracketingNonlinearSolve}/src/muller.jl (76%) create mode 100644 lib/BracketingNonlinearSolve/test/muller_tests.jl delete mode 100644 lib/SimpleNonlinearSolve/test/core/muller_tests.jl diff --git a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl index ae95eecba..3367c71c2 100644 --- 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, Muller, ITP, Ridder end diff --git a/lib/SimpleNonlinearSolve/src/muller.jl b/lib/BracketingNonlinearSolve/src/muller.jl similarity index 76% rename from lib/SimpleNonlinearSolve/src/muller.jl rename to lib/BracketingNonlinearSolve/src/muller.jl index ed8874e8e..ce4acafb5 100644 --- a/lib/SimpleNonlinearSolve/src/muller.jl +++ b/lib/BracketingNonlinearSolve/src/muller.jl @@ -1,18 +1,18 @@ """ - 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) @@ -20,6 +20,9 @@ function SciMLBase.solve(prob::NonlinearProblem, alg::SimpleMuller, args...; 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ᵢ₋₂ diff --git a/lib/BracketingNonlinearSolve/test/muller_tests.jl b/lib/BracketingNonlinearSolve/test/muller_tests.jl new file mode 100644 index 000000000..69144b841 --- /dev/null +++ b/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 diff --git a/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl b/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl index 1e284502f..495cadef2 100644 --- a/lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl +++ b/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 diff --git a/lib/SimpleNonlinearSolve/test/core/muller_tests.jl b/lib/SimpleNonlinearSolve/test/core/muller_tests.jl deleted file mode 100644 index fab2fbd0e..000000000 --- a/lib/SimpleNonlinearSolve/test/core/muller_tests.jl +++ /dev/null @@ -1,56 +0,0 @@ -@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 From ccae74424a0345c57c5ab7bb779acc78f8811917 Mon Sep 17 00:00:00 2001 From: Fabian Gittins Date: Sat, 5 Apr 2025 14:26:26 +0200 Subject: [PATCH 3/5] Fix bug with complex roots --- lib/BracketingNonlinearSolve/src/muller.jl | 4 ++-- .../test/muller_tests.jl | 16 ++++++++++++++++ 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/lib/BracketingNonlinearSolve/src/muller.jl b/lib/BracketingNonlinearSolve/src/muller.jl index ce4acafb5..6773ba4cf 100644 --- a/lib/BracketingNonlinearSolve/src/muller.jl +++ b/lib/BracketingNonlinearSolve/src/muller.jl @@ -20,8 +20,8 @@ function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Muller, args...; xᵢ₊₁, fxᵢ₊₁ = xᵢ₋₂, fxᵢ₋₂ - abstol = NonlinearSolveBase.get_tolerance( - xᵢ₋₂, abstol, promote_type(eltype(xᵢ₋₂), eltype(xᵢ))) + abstol = abs(NonlinearSolveBase.get_tolerance( + xᵢ₋₂, abstol, promote_type(eltype(xᵢ₋₂), eltype(xᵢ)))) for _ ∈ 1:maxiters q = (xᵢ - xᵢ₋₁)/(xᵢ₋₁ - xᵢ₋₂) diff --git a/lib/BracketingNonlinearSolve/test/muller_tests.jl b/lib/BracketingNonlinearSolve/test/muller_tests.jl index 69144b841..e1cf519d6 100644 --- a/lib/BracketingNonlinearSolve/test/muller_tests.jl +++ b/lib/BracketingNonlinearSolve/test/muller_tests.jl @@ -53,4 +53,20 @@ @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) + sol = solve(prob, Muller()) + + @test sol.u ≈ (-1 + √3*im)/2 + + tspan = (-1.0, -1.0*im) + prob = IntervalNonlinearProblem{false}(f, tspan) + sol = solve(prob, Muller()) + + @test sol.u ≈ (-1 - √3*im)/2 + end end From 297db290172c07fc464419342d66d9e8a11d1bc3 Mon Sep 17 00:00:00 2001 From: Fabian Gittins Date: Sat, 5 Apr 2025 17:10:59 +0200 Subject: [PATCH 4/5] Incorporate Muller into test suite Add `left` and `right` fields to solution --- .../src/BracketingNonlinearSolve.jl | 2 +- lib/BracketingNonlinearSolve/src/muller.jl | 6 ++++-- lib/BracketingNonlinearSolve/test/rootfind_tests.jl | 8 ++++---- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl index 3367c71c2..64db84621 100644 --- a/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl +++ b/lib/BracketingNonlinearSolve/src/BracketingNonlinearSolve.jl @@ -45,6 +45,6 @@ end @reexport using SciMLBase, NonlinearSolveBase -export Alefeld, Bisection, Brent, Falsi, Muller, 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 index 6773ba4cf..5f1291a35 100644 --- a/lib/BracketingNonlinearSolve/src/muller.jl +++ b/lib/BracketingNonlinearSolve/src/muller.jl @@ -43,7 +43,8 @@ function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Muller, args...; # Termination Check if abstol ≥ abs(fxᵢ₊₁) return SciMLBase.build_solution(prob, alg, xᵢ₊₁, fxᵢ₊₁; - retcode = ReturnCode.Success) + retcode = ReturnCode.Success, + left = xᵢ₊₁, right = xᵢ₊₁) end xᵢ₋₂, xᵢ₋₁, xᵢ = xᵢ₋₁, xᵢ, xᵢ₊₁ @@ -51,5 +52,6 @@ function CommonSolve.solve(prob::IntervalNonlinearProblem, alg::Muller, args...; end return SciMLBase.build_solution(prob, alg, xᵢ₊₁, fxᵢ₊₁; - retcode = ReturnCode.MaxIters) + retcode = ReturnCode.MaxIters, + left = xᵢ₊₁, right = xᵢ₊₁) end diff --git a/lib/BracketingNonlinearSolve/test/rootfind_tests.jl b/lib/BracketingNonlinearSolve/test/rootfind_tests.jl index 9161bdfe6..e8666a8c1 100644 --- 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 From adfde65a2dc8dce6216dc59476e7b7726028d2d1 Mon Sep 17 00:00:00 2001 From: Fabian Gittins Date: Sun, 6 Apr 2025 14:33:03 +0200 Subject: [PATCH 5/5] Add functionality to specify `middle` guess --- lib/BracketingNonlinearSolve/src/muller.jl | 17 ++++-- .../test/muller_tests.jl | 54 +++++++++++++------ 2 files changed, 52 insertions(+), 19 deletions(-) diff --git a/lib/BracketingNonlinearSolve/src/muller.jl b/lib/BracketingNonlinearSolve/src/muller.jl index 5f1291a35..edde280de 100644 --- a/lib/BracketingNonlinearSolve/src/muller.jl +++ b/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) diff --git a/lib/BracketingNonlinearSolve/test/muller_tests.jl b/lib/BracketingNonlinearSolve/test/muller_tests.jl index e1cf519d6..f4800b9d4 100644 --- a/lib/BracketingNonlinearSolve/test/muller_tests.jl +++ b/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