Implementation of Ridder

Author: Deltadahl

src/SimpleNonlinearSolve.jl

@@ -22,6 +22,7 @@ include("raphson.jl")
 include("broyden.jl")
 include("klement.jl")
 include("trustRegion.jl")
+include("ridder.jl")
 include("ad.jl")
 
 import SnoopPrecompile
@@ -43,12 +44,12 @@ 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)
+    for alg in (Bisection, Falsi, Ridder)
         solve(prob_brack, alg(), abstol = T(1e-2))
     end
 end end
 
 # DiffEq styled algorithms
-export Bisection, Broyden, Falsi, Klement, SimpleNewtonRaphson, SimpleTrustRegion
+export Bisection, Broyden, Falsi, Klement, Ridder, SimpleNewtonRaphson, SimpleTrustRegion
 
 end # module

src/ridder.jl

@@ -0,0 +1,84 @@
+"""
+`Ridder()`
+
+A non-allocating ridder method
+
+"""
+struct Ridder <: AbstractBracketingAlgorithm end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...;
+                         maxiters = 1000,
+                         kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+
+    if iszero(fl)
+        return SciMLBase.build_solution(prob, alg, left, fl;
+                                        retcode = ReturnCode.ExactSolutionLeft, left = left,
+                                        right = right)
+    end
+
+    xo = oftype(left, Inf)
+    i = 1
+    if !iszero(fr)
+        while i < maxiters
+            mid = (left + right) / 2
+            (mid == left || mid == right) &&
+                return SciMLBase.build_solution(prob, alg, left, fl;
+                                                retcode = ReturnCode.FloatingPointLimit,
+                                                left = left, right = right)
+            fm = f(mid)
+            s = sqrt(fm^2 - fl * fr)
+            iszero(s) &&
+                return SciMLBase.build_solution(prob, alg, left, fl;
+                                                retcode = ReturnCode.Failure,
+                                                left = left, right = right)
+            x = mid + (mid - left) * sign(fl - fr) * fm / s
+            fx = f(x)
+            xo = x
+            if iszero(fx)
+                right = x
+                fr = fx
+                break
+            end
+            if sign(fx) != sign(fm)
+                left = mid
+                fl = fm
+                right = x
+                fr = fx
+            elseif sign(fx) != sign(fl)
+                right = x
+                fr = fx
+            else
+                @assert sign(fx) != sign(fr)
+                left = x
+                fl = fx
+            end
+            i += 1
+        end
+    end
+
+    while i < maxiters
+        mid = (left + right) / 2
+        (mid == left || mid == right) &&
+            return SciMLBase.build_solution(prob, alg, left, fl;
+                                            retcode = ReturnCode.FloatingPointLimit,
+                                            left = left, right = right)
+        fm = f(mid)
+        if iszero(fm)
+            right = mid
+            fr = fm
+        elseif sign(fm) == sign(fl)
+            left = mid
+            fl = fm
+        else
+            right = mid
+            fr = fm
+        end
+        i += 1
+    end
+
+    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
+                                    left = left, right = right)
+end

test/basictests.jl

@@ -109,10 +109,22 @@ for p in 1.1:0.1:100.0
     @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
 end
 
+# Ridder
+g = function (p)
+    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
+    sol = solve(probN, Ridder())
+    return sol.left
+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]
-for alg in [Bisection(), Falsi()]
+for alg in [Bisection(), Falsi(), Ridder()]
     global g, p
     g = function (p)
         probN = IntervalNonlinearProblem{false}(f, tspan, p)
@@ -173,6 +185,10 @@ probB = IntervalNonlinearProblem(f, tspan)
 sol = solve(probB, Falsi())
 @test sol.left ≈ sqrt(2.0)
 
+# Ridder
+sol = solve(probB, Ridder())
+@test sol.left ≈ sqrt(2.0)
+
 sol = solve(probB, Bisection())
 @test sol.left ≈ sqrt(2.0)