Rigorous computation of function range

Download the notebook for this tutorial.

Introduction

The range of a function ff is defined as the set f(X)={f(x)xX}f(X) = \{f(x) | x\in X\}, where XX 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 f(x)=x2+2xf(x)=x^2+2x over the interval X=[5,5]X=[-5, 5].

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

x2+2x=x2+2x+11=(x+1)21 x^2+2x=x^2+2x+1-1=(x+1)^2-1

The variable xx 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 XX into NN 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 ff 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 ff over an interval XX. 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 g(x)=k=15kxsin(k(x3)3)g(x) = -\sum_{k=1}^5kx\sin\left(\frac{k(x-3)}{3}\right) over the interval X=[10,10]X=[-10,10].

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