setup with IntervalNonlinearProblem

Author: ChrisRackauckas

Project.toml

@@ -1,6 +1,6 @@
 name = "SimpleNonlinearSolve"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-authors = ["Kanav Gupta <kanav0610@gmail.com>"]
+authors = ["SciML"]
 version = "0.1.0"
 
 [deps]
@@ -21,7 +21,7 @@ FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 RecursiveArrayTools = "2"
 Reexport = "0.2, 1"
-SciMLBase = "1.32"
+SciMLBase = "1.73"
 Setfield = "0.7, 0.8, 1"
 StaticArrays = "0.12,1.0"
 UnPack = "1.0"

src/ad.jl

@@ -1,9 +1,15 @@
 function scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     f = prob.f
     p = value(prob.p)
-    u0 = value(prob.u0)
 
-    newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
+    if prob isa IntervalNonlinearProblem
+        tspan = value(prob.tspan)
+        newprob = IntervalNonlinearProblem(f, tspan, p; prob.kwargs...)
+    else
+        u0 = value(prob.u0)
+        newprob = NonlinearProblem(f, u0, p; prob.kwargs...)
+    end
+
     sol = solve(newprob, alg, args...; kwargs...)
 
     uu = sol.u
@@ -39,7 +45,8 @@ end
 
 # avoid ambiguities
 for Alg in [Bisection]
-    @eval function SciMLBase.solve(prob::NonlinearProblem{uType, iip, <:Dual{T, V, P}},
+    @eval function SciMLBase.solve(prob::IntervalNonlinearProblem{uType, iip,
+                                                                  <:Dual{T, V, P}},
                                    alg::$Alg, args...;
                                    kwargs...) where {uType, iip, T, V, P}
         sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
@@ -49,8 +56,12 @@ for Alg in [Bisection]
                                         right = Dual{T, V, P}(sol.right, partials))
         #return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode, sol.resid)
     end
-    @eval function SciMLBase.solve(prob::NonlinearProblem{uType, iip,
-                                                          <:AbstractArray{<:Dual{T, V, P}}},
+    @eval function SciMLBase.solve(prob::IntervalNonlinearProblem{uType, iip,
+                                                                  <:AbstractArray{
+                                                                                  <:Dual{T,
+                                                                                         V,
+                                                                                         P}
+                                                                                  }},
                                    alg::$Alg, args...;
                                    kwargs...) where {uType, iip, T, V, P}
         sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)

src/bisection.jl

@@ -7,10 +7,11 @@ function Bisection(; exact_left = false, exact_right = false)
     Bisection(exact_left, exact_right)
 end
 
-function SciMLBase.solve(prob::NonlinearProblem, alg::Bisection, args...; maxiters = 1000,
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args...;
+                         maxiters = 1000,
                          kwargs...)
     f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.u0
+    left, right = prob.tspan
     fl, fr = f(left), f(right)
 
     if iszero(fl)

src/falsi.jl

@@ -1,9 +1,10 @@
 struct Falsi <: AbstractBracketingAlgorithm end
 
-function SciMLBase.solve(prob::NonlinearProblem, alg::Falsi, args...; maxiters = 1000,
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...;
+                         maxiters = 1000,
                          kwargs...)
     f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.u0
+    left, right = prob.tspan
     fl, fr = f(left), f(right)
 
     if iszero(fl)
@@ -15,14 +16,14 @@ function SciMLBase.solve(prob::NonlinearProblem, alg::Falsi, args...; maxiters =
     i = 1
     if !iszero(fr)
         while i < maxiters
-            if nextfloat_tdir(left, prob.u0...) == right
+            if nextfloat_tdir(left, prob.tspan...) == right
                 return SciMLBase.build_solution(prob, alg, left, fl;
                                                 retcode = ReturnCode.FloatingPointLimit,
                                                 left = left, right = right)
             end
             mid = (fr * left - fl * right) / (fr - fl)
             for i in 1:10
-                mid = max_tdir(left, prevfloat_tdir(mid, prob.u0...), prob.u0...)
+                mid = max_tdir(left, prevfloat_tdir(mid, prob.tspan...), prob.tspan...)
             end
             if mid == right || mid == left
                 break

src/raphson.jl

@@ -1,13 +1,14 @@
 struct SimpleNewtonRaphson{CS, AD, FDT} <: AbstractNewtonAlgorithm{CS, AD, FDT}
     function SimpleNewtonRaphson(; chunk_size = Val{0}(), autodiff = Val{true}(),
-                                   diff_type = Val{:forward})
+                                 diff_type = Val{:forward})
         new{SciMLBase._unwrap_val(chunk_size), SciMLBase._unwrap_val(autodiff),
-        SciMLBase._unwrap_val(diff_type)}()
+            SciMLBase._unwrap_val(diff_type)}()
     end
 end
 
 function SciMLBase.solve(prob::NonlinearProblem,
-                         alg::SimpleNewtonRaphson, args...; xatol = nothing, xrtol = nothing,
+                         alg::SimpleNewtonRaphson, args...; xatol = nothing,
+                         xrtol = nothing,
                          maxiters = 1000, kwargs...)
     f = Base.Fix2(prob.f, prob.p)
     x = float(prob.u0)

test/basictests.jl

@@ -57,10 +57,10 @@ for p in 1.1:0.1:100.0
     @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
 end
 
-u0 = (1.0, 20.0)
+tspan = (1.0, 20.0)
 # Falsi
 g = function (p)
-    probN = NonlinearProblem{false}(f, typeof(p).(u0), p)
+    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
     sol = solve(probN, Falsi())
     return sol.left
 end
@@ -70,13 +70,13 @@ for p in 1.1:0.1:100.0
     @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
 end
 
-f, u0 = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
+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()]
     global g, p
     g = function (p)
-        probN = NonlinearProblem{false}(f, u0, p)
+        probN = IntervalNonlinearProblem{false}(f, tspan, p)
         sol = solve(probN, Bisection())
         return [sol.left]
     end
@@ -115,8 +115,8 @@ for u0 in [1.0, [1, 1.0]]
 end
 
 # Bisection Tests
-f, u0 = (u, p) -> u .* u .- 2.0, (1.0, 2.0)
-probB = NonlinearProblem(f, u0)
+f, tspan = (u, p) -> u .* u .- 2.0, (1.0, 2.0)
+probB = IntervalNonlinearProblem(f, tspan)
 
 # Falsi
 sol = solve(probB, Falsi())
@@ -135,7 +135,7 @@ f = function (u, p)
         return 0.0
     end
 end
-probB = NonlinearProblem(f, (0.0, 4.0))
+probB = IntervalNonlinearProblem(f, (0.0, 4.0))
 
 sol = solve(probB, Bisection(; exact_left = true))
 @test f(sol.left, nothing) < 0.0