inital commit, started to implement Broyden

Author: Deltadahl

src/SimpleNonlinearSolve.jl

@@ -4,6 +4,7 @@ using Reexport
 using FiniteDiff, ForwardDiff
 using ForwardDiff: Dual
 using StaticArraysCore
+using LinearAlgebra # TODO check if it is ok to add this
 import ArrayInterfaceCore
 
 @reexport using SciMLBase
@@ -18,12 +19,13 @@ include("bisection.jl")
 include("falsi.jl")
 include("raphson.jl")
 include("ad.jl")
+include("broyden.jl")
 
 import SnoopPrecompile
 
 SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
     prob_no_brack = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2))
-    for alg in (SimpleNewtonRaphson,)
+    for alg in (SimpleNewtonRaphson, Broyden)
         solve(prob_no_brack, alg(), tol = T(1e-2))
     end
 
@@ -44,6 +46,6 @@ SnoopPrecompile.@precompile_all_calls begin for T in (Float32, Float64)
 end end
 
 # DiffEq styled algorithms
-export Bisection, Falsi, SimpleNewtonRaphson
+export Bisection, Broyden, Falsi, SimpleNewtonRaphson
 
 end # module

src/broyden.jl

@@ -0,0 +1,43 @@
+# TODO add docstrings
+
+# TODO check what the supertype should be
+# TODO check if this should be defined as in raphson.jl
+struct Broyden <: AbstractSimpleNonlinearSolveAlgorithm end
+
+function SciMLBase.solve(prob::NonlinearProblem,
+                         alg::Broyden, args...; abstol = nothing,
+                         reltol = nothing,
+                         maxiters = 1000, kwargs...)
+    f = Base.Fix2(prob.f, prob.p)
+    xₙ = float(prob.u0)
+    T = typeof(xₙ)
+    J⁻¹ = Matrix{T}(I, length(xₙ), length(xₙ))
+
+    if SciMLBase.isinplace(prob)
+        error("Broyden 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)
+
+    for _ in 1:maxiters
+        # TODO check if nameing with heavy use of subscrips is ok
+        fₙ = f(xₙ)
+        xₙ₊₁ = xₙ + J⁻¹ * fₙ
+        Δxₙ = xₙ₊₁ - xₙ
+        Δfₙ = f(xₙ₊₁) - fₙ
+        J⁻¹ .+= ((Δxₙ .- J⁻¹ * Δfₙ) ./ (Δxₙ' * J⁻¹ * Δfₙ)) * (Δxₙ' * J⁻¹)
+
+        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ₙ₊₁
+    end
+
+    @show xₙ, fₙ
+    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
+end

test/basictests.jl

@@ -77,7 +77,7 @@ for alg in [Bisection(), Falsi()]
     global g, p
     g = function (p)
         probN = IntervalNonlinearProblem{false}(f, tspan, p)
-        sol = solve(probN, Bisection())
+        sol = solve(probN, Bisection())  # TODO check if "alg" should replace "Bisection()"
         return [sol.left]
     end
 
@@ -102,13 +102,15 @@ probN = NonlinearProblem(f, u0)
 @test solve(probN, SimpleNewtonRaphson(); immutable = false).u[end] ≈ sqrt(2.0)
 @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0)
 @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0)
+# TODO check why the 2 lines above are identical
 
 for u0 in [1.0, [1, 1.0]]
     local f, probN, sol
     f = (u, p) -> u .* u .- 2.0
     probN = NonlinearProblem(f, u0)
     sol = sqrt(2) * u0
 
+    # TODO check why the two lines below are identical
     @test solve(probN, SimpleNewtonRaphson()).u ≈ sol
     @test solve(probN, SimpleNewtonRaphson()).u ≈ sol
     @test solve(probN, SimpleNewtonRaphson(; autodiff = false)).u ≈ sol