Download the notebook for this tutorial.
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
To begin with, let us compute the range of the function over the interval .
f(x) = x^2+2x X = -5..5 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 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
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)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 = 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
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) (Y, N, err) plot(f, -5, 5, legend=false) plot!(IntervalBox(X, Y))
(Y, N, err) = ([-1.2066, 35], 96, 0.005721096543504041)
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) (Y, N, err) plot(g, -10, 10, legend=false) plot!(IntervalBox(X, Y))
(Y, N, err) = ([-55.8697, 33.9264], 1536, 0.009006235498525415)
anim = 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