Merge pull request #18 from CCsimon123/main

Implemented a Klement solver

lib/SimpleNonlinearSolve/src/SimpleNonlinearSolve.jl

@@ -20,6 +20,7 @@ include("falsi.jl")
 include("raphson.jl")
 include("ad.jl")
 include("broyden.jl")
+include("klement.jl")
 
 import SnoopPrecompile
 
@@ -46,6 +47,6 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
 end end
 
 # DiffEq styled algorithms
-export Bisection, Broyden, Falsi, SimpleNewtonRaphson
+export Bisection, Broyden, Falsi, Klement, SimpleNewtonRaphson
 
 end # module

lib/SimpleNonlinearSolve/src/klement.jl

@@ -0,0 +1,62 @@
+"""
+```julia
+Klement()
+```
+
+A low-overhead implementation of [Klement](https://jatm.com.br/jatm/article/view/373).
+This method is non-allocating on scalar and static array problems.
+"""
+struct Klement <: AbstractSimpleNonlinearSolveAlgorithm end
+
+function SciMLBase.solve(prob::NonlinearProblem,
+                         alg::Klement, args...; abstol = nothing,
+                         reltol = nothing,
+                         maxiters = 1000, kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    x = float(prob.u0)
+    fₙ = f(x)
+    T = eltype(x)
+    J = ArrayInterfaceCore.zeromatrix(x) + I
+
+    if SciMLBase.isinplace(prob)
+        error("Klement currently only supports out-of-place nonlinear problems")
+    end
+
+    atol = abstol !== nothing ? abstol :
+           real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5)
+    rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5)
+
+    xₙ = x
+    xₙ₋₁ = x
+    fₙ₋₁ = fₙ
+    for _ in 1:maxiters
+        xₙ = xₙ₋₁ - inv(J) * fₙ₋₁
+        fₙ = f(xₙ)
+        Δxₙ = xₙ - xₙ₋₁
+        Δfₙ = fₙ - fₙ₋₁
+
+        # Prevent division by 0
+        denominator = max.(J' .^ 2 * Δxₙ .^ 2, 1e-9)
+
+        k = (Δfₙ - J * Δxₙ) ./ denominator
+        J += (k * Δxₙ' .* J) * J
+
+        # Prevent inverting singular matrix
+        if det(J) ≈ 0
+            J = ArrayInterfaceCore.zeromatrix(x) + I
+        end
+
+        iszero(fₙ) &&
+            return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
+                                            retcode = ReturnCode.Success)
+
+        if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol)
+            return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
+                                            retcode = ReturnCode.Success)
+        end
+        xₙ₋₁ = xₙ
+        fₙ₋₁ = fₙ
+    end
+
+    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
+end

lib/SimpleNonlinearSolve/test/basictests.jl

@@ -35,13 +35,24 @@ sol = benchmark_scalar(sf, csu0)
 @test sol.u * sol.u - 2 < 1e-9
 @test (@ballocated benchmark_scalar(sf, csu0)) == 0
 
+# Klement
+function benchmark_scalar(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    sol = (solve(probN, Klement()))
+end
+
+sol = benchmark_scalar(sf, csu0)
+@test sol.retcode === ReturnCode.Success
+@test sol.u * sol.u - 2 < 1e-9
+@test (@ballocated benchmark_scalar(sf, csu0)) == 0
+
 # AD Tests
 using ForwardDiff
 
 # Immutable
 f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0]
 
-for alg in [SimpleNewtonRaphson(), Broyden()]
+for alg in [SimpleNewtonRaphson(), Broyden(), Klement()]
     g = function (p)
         probN = NonlinearProblem{false}(f, csu0, p)
         sol = solve(probN, alg, tol = 1e-9)
@@ -56,7 +67,7 @@ end
 
 # Scalar
 f, u0 = (u, p) -> u * u - p, 1.0
-for alg in [SimpleNewtonRaphson(), Broyden()]
+for alg in [SimpleNewtonRaphson(), Broyden(), Klement()]
     g = function (p)
         probN = NonlinearProblem{false}(f, oftype(p, u0), p)
         sol = solve(probN, alg)
@@ -97,7 +108,7 @@ for alg in [Bisection(), Falsi()]
     @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
 end
 
-for alg in [SimpleNewtonRaphson(), Broyden()]
+for alg in [SimpleNewtonRaphson(), Broyden(), Klement()]
     global g, p
     g = function (p)
         probN = NonlinearProblem{false}(f, 0.5, p)
@@ -120,6 +131,9 @@ probN = NonlinearProblem(f, u0)
 @test solve(probN, Broyden()).u[end] ≈ sqrt(2.0)
 @test solve(probN, Broyden(); immutable = false).u[end] ≈ sqrt(2.0)
 
+@test solve(probN, Klement()).u[end] ≈ sqrt(2.0)
+@test solve(probN, Klement(); immutable = false).u[end] ≈ sqrt(2.0)
+
 for u0 in [1.0, [1, 1.0]]
     local f, probN, sol
     f = (u, p) -> u .* u .- 2.0
@@ -131,6 +145,8 @@ for u0 in [1.0, [1, 1.0]]
     @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u ≈ sol
 
     @test solve(probN, Broyden()).u ≈ sol
+
+    @test solve(probN, Klement()).u ≈ sol
 end
 
 # Bisection Tests