Rigorous computation of function range
Download the notebook for this tutorial.
Introduction
The range of a function is defined as the set , where is the domain of the function. This tutorial will show how to to exploit interval arithmetic to give a rigorous estimate of the function range.
First, let's import all the packages we need
using IntervalArithmetic, Plots
A first example
To begin with, let us compute the range of the function over the interval .
f(x) = x^2+2x
X = -5..5
@show range_1_interval = f(X)
plot(f, -5, 5, leg=false)
plot!(IntervalBox(X, f(X)))
range_1_interval = f(X) = [-10, 35]
As can be noticed, in this case we overestimate the range of the function. This is a common issue in interval arithmetic and it is known as dependency problem. If a variable is present more than once in the expression, then evaluating the function over an interval will produce larger interval. For this simple function, the problem can be solved with simple algebraic manipulation observing that
The variable is now present only once and no overestimation occurs, as the following code confirms
f1(x) = (x+1)^2-1
@show range_f = f1(X)
plot(f1, -5, 5, leg=false)
plot!(IntervalBox(X, f1(X)))
range_f = f1(X) = [-1, 35]
Despite for this simple example we could analytically solve the dependency problem, this is not possible in general. The strategy to reduce overestimation is divide and conquer. The domain is divided into smaller intervals and the function is evaluated at each interval. The process of partitioning the interval into smaller piecing is called "mincing" and IntervalArithmetic
offers the function mince(X, N)
to partition the interval into intervals. This function returns an array of intervals. For example using 10 intervals
Xs = mince(X, 10)
Y = f.(Xs)
plot(f, -5, 5, leg=false)
plot!(IntervalBox.(Xs, f.(Xs)))
Assuming the function is continuous, we can now compute the function range by taking the union of the computed intervals. This can be done using the function reduce
range_10_intervals = reduce(∪, Y)
overestimate_1_interval = (diam(range_1_interval)-diam(range_f))/diam(range_f)
overestimate_10_intervals = (diam(range_10_intervals)-diam(range_f))/diam(range_f)
@show overestimate_1_interval, overestimate_10_intervals
(overestimate_1_interval, overestimate_10_intervals) = (0.25, 0.05555555555555555)
Observe that increasing the number of intervals from 1 to 10 reduced the relative error from 25% to 5.5%. The following animation shows how the range approximation improves doubling the number of intervals at each iteration.
anim = @animate for i in 0:10
Xs = mince(X, 2^i)
plot(f, -5, 5, leg=false, ylims=(-5, 40), xlims=(X.lo, X.hi), lw=2)
plot!(IntervalBox.(Xs, f.(Xs)))
end
We are now ready to write our function range(f, X, tol)
which estimates the range of a function over an interval . The function will take a third parameter an error tolerance tol
and keep increasing the number of intervals until the relative change will become smaller than tol
.
function range(f, X, N, tol=0.01)
Xs = mince(X, N)
Ynext = f.(Xs)
Ynext = reduce(∪, Ynext)
err = tol+1
while err > tol
Yprev = Ynext
N *= 2
Xs = mince(X, N)
Ynext = f.(Xs)
Ynext = reduce(∪, Ynext)
err = (diam(Yprev) - diam(Ynext))/diam(Yprev)
end
return (Ynext, N, err)
end
range (generic function with 2 methods)
Now we are ready to compute the range of our function
Y, N, err = range(f, X, 3)
@show (Y, N, err)
plot(f, -5, 5, legend=false)
plot!(IntervalBox(X, Y))
(Y, N, err) = ([-1.2066, 35], 96, 0.005721096543504041)
A more challenging example (1D)
Now that we have developed our tool to estimate a function range, we are ready to test it with more challenging functions. Let's estimate the range of the function over the interval .
X = -10..10
g(x) = -sum([k*x*sin(k*(x-3)/3) for k in 1:5])
Y, N, err = range(g, X, 3)
@show (Y, N, err)
plot(g, -10, 10, legend=false)
plot!(IntervalBox(X, Y))
(Y, N, err) = ([-55.8697, 33.9264], 1536, 0.009006235498525415)
anim = @animate for i in 2 .^(0:10)
Xs = mince(X, i)
plot(g, -10, 10, leg=false, xlims=(X.lo, X.hi), ylims=(-60, 50), lw=2)
plot!(IntervalBox.(Xs, g.(Xs)))
end