Cannot find root in two dimensional case

[newptcai]

You can try the following code

using IntervalArithmetic, IntervalRootFinding, IntervalOptimisation, StaticArrays
a = 1/2*log(@biginterval(2))
f(x, y)=sinh(x) - 1/2 * (cosh(x)/cosh(x-y))^(x/y)

Then

rts1 = roots(x->f(x, a.lo), @biginterval(0, 1))

finds a root

1-element Array{Root{Interval{BigFloat}},1}:
 Root([0.597845, 0.597846]₂₅₆, :unique)

But the two dimensional version

g((x, y)) = SVector(f(x,y))
roots(g, @biginterval(0,1) × a)

fails find anything.

[dpsanders]

As far as I can see, g is a function from R^2 to R, so it does not make sense to talk about its (isolated) roots -- in general the set where g == 0 will be a curve.

If you are looking to find that curve, you need https://github.com/JuliaIntervals/IntervalConstraintProgramming.jl

It may be that a zero_set function which chooses one of these options would be useful in the future.

[dpsanders]

Ah I see that a is a tiny interval. But still.

[dpsanders]

It is perhaps concerning that it just returns an empty list. But I think the routines are designed to only work when the input dimension = output dimension.

This should probably be checked explicitly. cc @Kolaru

[dpsanders]

There is also a parametric interval Newton method that in principle could prove that there is a unique solution for every y in a, but this is not yet implemented.

I believe @mewilhel has worked on that.

[newptcai]

I don't need it to work for every pair of x and a. I only want to show that the function f(x, a) has one unique root between 0 and 1. But then a is not a rational number, if I approximate it by any float, then I cannot really say that this is true.

Is there any way to actually prove this with numeric methods?

[dpsanders]

Ah, I see.

You can just put the interval right into the roots command!

julia> rts1 = roots(x->f(x, a), @biginterval(0, 1))
1-element Array{Root{Interval{BigFloat}},1}:
 Root([0.597845, 0.597846]₂₅₆, :unique)

[dpsanders]

By the way, you don't need BigFloats for this:

julia> a = (1/2) * log(2..2)
[0.346573, 0.346574]

julia> roots(x -> f(x, a), 0..1)
1-element Array{Root{Interval{Float64}},1}:
 Root([0.597845, 0.597846], :unique)

[dpsanders]

You can also use @interval outside the expression, which will "magically" wrap each constant inside an interval.

julia> @format full
Display parameters:
- format: full
- decorations: false
- significant figures: 6

julia> a = (1/2) * log(interval(2))
Interval(0.34657359027997264, 0.3465735902799727)

julia> a = @interval 1/2 * log(2)
Interval(0.34657359027997264, 0.3465735902799727)

[newptcai]

I see. Thanks. That solves my problem.

But I tried something similar, but I got an error:

f(x::Interval{BigFloat}) = f(x, a)
roots(f, @biginterval(0, 1))

MethodError: no method matching f(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f),Interval{BigFloat}},Interval{BigFloat},1})

[dpsanders]

As it says, there's no way it can apply f to a Dual number, which is what ForwardDiff uses to calculate the derivatives needed by the Newton method.

That's because when you defined your new f, you restricted the type of argument it accepts unnecessarily (to only accept an interval). Just leave off the restriction:

julia> f(x) = f(x, a)
f (generic function with 2 methods)

julia> roots(f, @biginterval(0, 1))
1-element Array{Root{Interval{BigFloat}},1}:
 Root(Interval(0.5978454838339249935256456020018315465663452836738495130173055992072950702816678, 0.597845484773853995337508310718679163254084878360925861130856532225309042455311), :unique)

[dpsanders]

You' re welcome. That' s a nice example function, by the way!

[Kolaru]

Technically Newton and Krawczyck methods require the Jacobian of the function to be invertible close to the solution. This can never be the case for function from R^2 to R, the Jacobian is not even a square matrix.

I wonder why it doesn't error, and I too believe we should explicitly fail here.