Add new root status when the presence (but not uniqueness) of a solution is known

[Kolaru]

For Krawczyk and Newton method, the fact that a contracted interval is subset of the original interval is a sufficient condition for the presence of a solution.

The proposition here is to add a new root status :exist that correspoond to this information.

This can be usefull when solving for equations with a continuous set of solution, for which this situation is common, it allows to prove the presence of a solution in some region.

See a dummy example below.

julia> import ForwardDiff: jacobian

julia> using StaticArrays

julia> f(xx) = SVector(xx[1] + xx[2], xx[1] + xx[2])
f (generic function with 1 method)

julia> rts = roots(f, IntervalBox(-1..1, 2), Krawczyk)
4559-element Array{Root{IntervalBox{2,Float64}},1}:
 Root([0.00721069, 0.00813498] × [-0.0078125, -0.00690255], :unknown)
 Root([-0.0078125, -0.00690255] × [0.00721069, 0.00813498], :unknown)
 Root([-0.000359588, 0.000564693] × [-0.000359588, 0.000564693], :unknown)
 Root([0.00339598, 0.00433482] × [-0.00411516, -0.00319087], :unknown)
 Root([0.00528843, 0.00624206] × [-0.00597828, -0.00505399], :unknown)
 Root([0.00624205, 0.0072107] × [-0.00690256, -0.00597827], :unknown)
 Root([0.00624205, 0.0072107] × [-0.0078125, -0.00690255], :unknown)
 Root([0.00528843, 0.00624206] × [-0.00690256, -0.00597827], :unknown)
 Root([0.00433481, 0.00528844] × [-0.005054, -0.00411515], :unknown)
 Root([0.00433481, 0.00528844] × [-0.00597828, -0.00505399], :unknown)
 ⋮
 Root([-0.324746, -0.323745] × [0.32355, 0.324535], :unknown)
 Root([-0.323243, -0.32273] × [0.322566, 0.323551], :unknown)
 Root([-0.323746, -0.323242] × [0.322566, 0.323551], :unknown)
 Root([-0.325745, -0.324745] × [0.324534, 0.325535], :unknown)
 Root([-0.332648, -0.331663] × [0.332452, 0.332957], :unknown)
 Root([-0.331664, -0.330664] × [0.331453, 0.332453], :unknown)
 Root([-0.332648, -0.331663] × [0.331453, 0.332453], :unknown)
 Root([-0.333632, -0.332647] × [0.332956, 0.333468], :unknown)
 Root([-0.333632, -0.332647] × [0.332452, 0.332957], :unknown)
 Root([-0.331664, -0.330664] × [0.330453, 0.331454], :unknown)

julia> X = interval(rts[1])
[0.00433481, 0.00528844] × [-0.00597828, -0.00505399]

julia> C = Krawczyk(f, x -> jacobian(f, x))
Krawczyk{typeof(f),getfield(Main, Symbol("##9#10"))}(f, getfield(Main, Symbol("##9#10"))())

julia> status, CX = C(X, 1e-15)

julia> issubset(CX, X)
true

PS: I'm actually volunteer to implement it when I found a bit of time, but I prefer to first raise it as an issue as I suspect I may have overlooked something preventing this.

[dpsanders]

As far as I understand, if N(X) is a subset of X, this proves that there exists a solution inside X and that that solution is unique inside the interval X.

You're right that there are parametric versions of Newton that prove that there is a unique solution of F(x, p) for each p in the interval, but I do not know how those work. Is that what you are aiming for here?

[Kolaru]

As far as I understand, if N(X) is a subset of X, this proves that there exists a solution inside X and that that solution is unique inside the interval X.

My understanding is that there is a subtle:

• If N(X) is a subset of X it proves that a solution exists (by the fixpoint theorem).
• If N(X) is in the interior of X then it proves that this solution is unique (by consideration on the derivative of the function).

So what I would like is to be able to have a status representing the case where only the first condition is true.

[dpsanders]

That could well be. Do you have a reference for that?

[Kolaru]

Theorem 8.2. in Moore Introduction to Interval Analysis (2009) states it for Krawczyk method. He does not go into details for the Newton method though, but the fixpoint argument should hold as well. I should be able to find a reference in that case soon.

[Kolaru]

So I work a bit about the concept, and realised that in most case this will not be useful: if the set of solutions is continuous, the Jacobian is singular in the interval box (the singularity of the Jacobian is implied by the fact that the function does not vary along the set of solutions). As a consequence Newton and Krawczyck contractors always return the entire space and can never conclude to the existence of solution.

The only application is if there is a solution on the border of the interval, e.g. in the following

julia> roots(sin, 0..1)
1-element Array{Root{Interval{Float64}},1}:
 Root([0, 4.06359e-16], :unknown)

julia> C = Newton(sin, cos)
Newton{typeof(sin),typeof(cos)}(sin, cos)

julia> X = 0..0.01
[0, 0.0100001]

julia> CX = C(X, 1e-15)
(:unknown, [0, 2.03489e-08])

The fact that CX is a subset of X proves that CX contains a solution. Since it is rather unlikely, it may not be worth checking it explicitely.

[dpsanders]

Here is a reference that confirms your first statement in this issue:

https://www.researchgate.net/publication/285598481_Krawczyk_operator_revised