How to show there is no roots on unbounded domain? #149
Comments
|
You can use infinity as follows: julia> f(x) = x
f (generic function with 2 methods)
julia> roots(f, 1..∞) # type as \infty<TAB>. Or use Inf
0-element Array{Root{Interval{Float64}},1}The fact that the list of roots is empty means that there is no root. The second question is much harder. |
|
Your function is particularly nasty due to those hyperbolic functions, e.g. julia> @time roots(f, -5..6)
0.691032 seconds (4.82 M allocations: 148.984 MiB, 2.00% gc time)
1-element Array{Root{Interval{Float64}},1}:
Root([0.597845, 0.597846], :unique)
julia> @time roots(f, -5..7)
5.044978 seconds (34.90 M allocations: 1.053 GiB, 1.84% gc time)
1-element Array{Root{Interval{Float64}},1}:
Root([0.597845, 0.597846], :unique)As you increase the right limit of the interval the time to prove that there is a unique root increases dramatically fast. This is basically because of the cancellation in your function. Is there a way to rewrite the function somehow? In general, showing that a function is positive on a range using interval arithmetic should be "easy" by just evaluating the function over the interval, e.g. julia> g(x) = x^2 - 2
f (generic function with 2 methods)
julia> g(3..4)
[7, 14]This proves that the image of the function However, again the nature of your function makes this very difficult, e.g. for your function julia> f(3..4)
[-2.66365e+06, 27.2897]even though the true range of the function is much smaller than that! For example, it may be possible to approximate your function using a Taylor series in some range and bound that more easily; see |
|
Another option is to calculate the derivative and show that the derivative is negative over some region, implying that the function is decreasing there. |
|
So I transfer the problem to showing that the following is non-negative between Can we show something like it's positive on something like 1/2*t^2*(log(2)*log(t^2 - 1)/log(t) - 2*log(-2/(t^2 + 2) + 2))I tried |
|
When I evaluate that function at But it still seems to have the dependency problem. Basically when a variable is repeated more than once, interval arithmetic starts to give intervals that are too wide. |
|
You can always try to force the use of smaller intervals to try to avoid the dependency problem. In your original example: julia> f(x, y) = sinh(x) - 1/2 * (cosh(x)/cosh(x-y))^(x/y)
f (generic function with 1 method)
julia> f(x) = f(x, 0.5*log(2))
f (generic function with 2 methods)
julia> f(1..2)
[-329.129, 3.6074]But if instead of trying julia> bounds = 1:0.0001:2
1.0:0.0001:2.0
julia> Xs = Interval.(bounds[1:end-1], bounds[2:end]) ;
julia> res = f.(Xs) ;
julia> minimum(inf.(res))
0.19284681924470437It may take a while and may be more efficient to bisect the interval until the required precision is matched, rather than setting a constant step. However, you won't be able to do it up to infinity since the function is not well define there julia> f(Inf)
NaN |
|
Basically this bisection is what |
|
I managed to transform the problem into this showing that -(log(2*s + 1) - log(s + 1))*log(s)/log(2) + log(-s + 1)is non-negative on -(log(2*s + 1) - log(s + 1))*log(1/2/exp(2))/log(2) + log(-s + 1)which goes to 2*(s^2*log(2) + s*(2*log(2) + 1) - 1)/((2*s + 1)*(s + 1)*(s - 1)*log(2))This we can use So it seems that to prove such inequalities mathematically using interval arithmetic, one still has to
|
|
Actually f(s) = -(log(2*s + 1) - log(s + 1))*log(1/2/exp(2))/log(2) + log(-s + 1)is well-defined on julia> roots(f, -exp(-2)..exp(2))
2-element Array{Root{Interval{Float64}},1}:
Root([0.749847, 0.749848], :unique)
Root([-1.44908e-09, 1.38232e-09], :unique)We know analytically that f(0)=0. julia> (@interval f(exp(-2))) > 0
trueTherefore, we can conclude that f(x) >= 0 on 0 to e^-2. |
|
Yes you will often need to "massage" your problem first. It's unfortunate. There is some work going on on symbolic manipulation which might help, but it will still need to be driven by a human. Note that you can use julia> f(s) = -(log(2*s + 1) - log(s + 1))*log(1/2/exp(2))/log(2) + log(-s + 1)
f (generic function with 1 method)
julia> ForwardDiff.derivative(f, 0.1..0.2)
[0.768384, 2.12672] |
|
I have written this up -- https://newptcai.github.io/proving-an-inequality-with-julia-and-interval-arithmetic.html Hope it would be helpful for others. :) |
|
BTW, this is probably a bug in IntervalArithemtic julia> a = @biginterval log(2)/2
[0.346573, 0.346574]₂₅₆
julia> f(x) = sinh(x) - 1/2 * (cosh(x)/cosh(x-a))^(x/a)
f (generic function with 1 method)
julia> f(50)
[3.32447e+07, 3.32448e+07]₂₅₆I have tried Mathematica, Maple, and SageMath. All suggests that f(50) = 1.381*10^-20 |
julia> f(@biginterval(50), @biginterval(a))
[1.38165e-20, 1.38166e-20]₂₅₆
julia> f(big(50), @biginterval(a))
[1.38165e-20, 1.38166e-20]₂₅₆
In general IntervalArithmetic results are only guaranteed if all floating point number involved have been replaced by intervals. |
|
In my code above I used integer x=50 and a was a big interval. So couldn’t 50 be turned into big interval automatically? That was what I was expecting. |
|
It will be turn into a big interval when it first interact with an interval. But |
|
That's a great blog post! The way Julia works, it just does what you tell it. You start off with an integer x and the first thing it does is compute sinh(x) for that integer x, by converting x to a Float64. Only later do you ask it to subtract a from x, at which point it will convert x to an Interval. |
I see. Thanks. |
We must be using different versions. I got this julia> using IntervalArithmetic
julia> a = @biginterval log(2)/2
[0.346573, 0.346574]₂₅₆
julia> f(x, y) = sinh(x) - 1/2 * (cosh(x)/cosh(x-y))^(x/y)
f (generic function with 1 method)
julia> f(@biginterval(50), @biginterval(a))
[7.95003e-14, 7.95004e-14]₂₅₆
julia> f(@biginterval(50), a)
[7.95003e-14, 7.95004e-14]₂₅₆
julia> f(big(50), a)
[7.95003e-14, 7.95004e-14]₂₅₆ |
|
That's worrying. Which version and machine are you on? |
|
Status shell> lsb_release --all |
|
What happens if you update to the latest version of IntervalArithmetic? I don't see why it should make a difference, but... |
|
By the way, in Julia you can do versioninfo()to get info about your machine |
|
Also you can test this without needing to update all packages by making a new directory and doing ]activate .there, then adding |
|
OK weird, I see the same results as @newptcai when I make a new environment and load the released version of IntervalArithmetic |
Yeah, I think I have the latest version. |
|
What happens if you use the branch
|
|
Maybe there's a problem with the definition of |
Yes. This branch is correct now. |
|
OK thanks, will look into this. |
Let's take the simplest example. If I take
f(x) = xand I want to show that there is root of this function on 1 to infinity, how can I do this?A more complicated example is
Taking y equals
1/2*log(2), this function remains positive but goes to 0 asxgoes to infinity. Is there anyway that we can show this with interval arithmetics?The text was updated successfully, but these errors were encountered: