Fix Issue #34: Add new method Alefeld, Potra, Shi (1995)

Project.toml

@@ -43,4 +43,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]
+test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "StaticArrays", "NNlib"]

src/SimpleNonlinearSolve.jl

@@ -38,6 +38,7 @@ include("brent.jl")
 include("dfsane.jl")
 include("ad.jl")
 include("halley.jl")
+include("alefeld.jl")
 
 import SnoopPrecompile
 
@@ -59,13 +60,13 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
     =#
 
     prob_brack = IntervalNonlinearProblem{false}((u, p) -> u * u - p, T.((0.0, 2.0)), T(2))
-    for alg in (Bisection, Falsi, Ridder, Brent)
+    for alg in (Bisection, Falsi, Ridder, Brent, Alefeld)
         solve(prob_brack, alg(), abstol = T(1e-2))
     end
 end end
 
 # DiffEq styled algorithms
 export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, Halley, Klement,
-       Ridder, SimpleNewtonRaphson, SimpleTrustRegion
+       Ridder, SimpleNewtonRaphson, SimpleTrustRegion, Alefeld
 
 end # module

src/alefeld.jl

@@ -0,0 +1,165 @@
+"""
+`Alefeld()` 
+
+An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145/210089.210111).
+
+The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than 
+algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
+"""
+struct Alefeld <: AbstractBracketingAlgorithm end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem,
+                            alg::Alefeld, args...; abstol = nothing,
+                            reltol = nothing,
+                            maxiters = 1000, kwargs...)
+
+    f = Base.Fix2(prob.f, prob.p)
+    a, b = prob.tspan
+    c = a - (b - a) / (f(b) - f(a)) * f(a)
+    @show c
+
+    fc = f(c)
+    if iszero(fc)
+        return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+    end
+    a, b, d = _bracket(f, a, b, c)
+    e = 0   # Set e as 0 before iteration to avoid a non-value f(e)
+
+    # Begin of algorithm iteration
+    for i in 2:maxiters
+        # The first bracketing block
+        f₁, f₂, f₃, f₄ = f(a), f(b), f(d), f(e)
+        if i == 2 || (f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄)
+            c = _newton_quadratic(f, a, b, d, 2)
+        else 
+            c = _ipzero(f, a, b, d, e)
+            if (c - a) * (c - b) ≥ 0
+                c = _newton_quadratic(f, a, b, d, 2)
+            end
+        end 
+        ē, fc = d, f(c)   
+        iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+        ā, b̄, d̄ = _bracket(f, a, b, c) 
+
+        # The second bracketing block
+        f₁, f₂, f₃, f₄ = f(ā), f(b̄), f(d̄), f(ē)
+        if f₁ == f₂ || f₁ == f₃ || f₁ == f₄ || f₂ == f₃ || f₂ == f₄ || f₃ == f₄
+            c = _newton_quadratic(f, ā, b̄, d̄, 3)
+        else 
+            c = _ipzero(f, ā, b̄, d̄, ē)
+            if (c - ā) * (c - b̄) ≥ 0
+                c = _newton_quadratic(f, ā, b̄, d̄, 3)
+            end
+        end
+        fc = f(c)
+        iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+        ā, b̄, d̄ = _bracket(f, ā, b̄, c) 
+
+        # The third bracketing block
+        if abs(f(ā)) < abs(f(b̄))
+            u = ā
+        else
+            u = b̄
+        end
+        c = u - 2 * (b̄ - ā) / (f(b̄) - f(ā)) * f(u)
+        if (abs(c - u)) > 0.5 * (b̄ - ā)
+            c = 0.5 * (ā + b̄)
+        end
+        fc = f(c)
+        iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+        ā, b̄, d = _bracket(f, ā, b̄, c) 
+
+        # The last bracketing block
+        if b̄ - ā < 0.5 * (b - a)
+            a, b, e = ā, b̄, d̄
+        else
+            e = d
+            c = 0.5 * (ā + b̄)
+            fc = f(c)
+            iszero(fc) &&
+            return SciMLBase.build_solution(prob, alg, c, fc;
+                                        retcode = ReturnCode.Success, 
+                                        left = a,
+                                        right = b)
+            a, b, d = _bracket(f, ā, b̄, c)
+        end
+    end
+
+    # Reassign the value a, b, and c
+    if b == c
+        b = d
+    elseif a == c
+        a = d
+    end
+    fc = f(c)
+
+    # Reuturn solution when run out of max interation
+    return SciMLBase.build_solution(prob, alg, c, fc; retcode = ReturnCode.MaxIters,
+                                    left = a, right = b)
+end
+
+# Define subrotine function bracket, check fc before bracket to return solution
+function _bracket(f::F, a, b, c) where F
+    if iszero(f(c))
+        ā, b̄, d = a, b, c
+    else
+        if f(a) * f(c) < 0 
+            ā, b̄, d = a, c, b
+        elseif f(b) * f(c) < 0
+            ā, b̄, d = c, b, a
+        end
+    end
+
+    return ā, b̄, d 
+end
+
+# Define subrotine function newton quadratic, return the approximation of zero
+function _newton_quadratic(f::F, a, b, d, k) where F
+    A = ((f(d) - f(b)) / (d - b) - (f(b) - f(a)) / (b - a)) / (d - a) 
+    B = (f(b) - f(a)) / (b - a)
+
+    if iszero(A)
+        return a - (1 / B) * f(a)
+    elseif A * f(a) > 0
+        rᵢ₋₁ = a 
+    else 
+        rᵢ₋₁ = b
+    end 
+
+    for i in 1:k
+        rᵢ = rᵢ₋₁ - B * rᵢ₋₁ / (B + A * (2 * rᵢ₋₁ - a - b))
+        rᵢ₋₁ = rᵢ
+    end
+
+    return rᵢ₋₁
+end
+
+# Define subrotine function ipzero, also return the approximation of zero
+function _ipzero(f::F, a, b, c, d) where F
+    Q₁₁ = (c - d) * f(c) / (f(d) - f(c))
+    Q₂₁ = (b - c) * f(b) / (f(c) - f(b))
+    Q₃₁ = (a - b) * f(a) / (f(b) - f(a))
+    D₂₁ = (b - c) * f(c) / (f(c) - f(b))
+    D₃₁ = (a - b) * f(b) / (f(b) - f(a))
+    Q₂₂ = (D₂₁ - Q₁₁) * f(b) / (f(d) - f(b))
+    Q₃₂ = (D₃₁ - Q₂₁) * f(a) / (f(c) - f(a))
+    D₃₂ = (D₃₁ - Q₂₁) * f(c) / (f(c) - f(a))
+    Q₃₃ = (D₃₂ - Q₂₂) * f(a) / (f(d) - f(a))
+
+    return a + Q₃₁ + Q₃₂ + Q₃₃
+end

test/basictests.jl

@@ -210,6 +210,18 @@ for p in 1.1:0.1:100.0
     @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
 end
 
+# Alefeld
+g = function (p)
+    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
+    sol = solve(probN, Alefeld())
+    return sol.u
+end
+
+for p in 1.1:0.1:100.0
+    @test g(p) ≈ sqrt(p)
+    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+end
+
 f, tspan = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
 t = (p) -> [sqrt(p[2] / p[1])]
 p = [0.9, 50.0]