Interval root finding tutorial

Download the notebook for this tutorial.

Setup

The IntervalRootFinding.jl package can be installed with

julia> using Pkg; Pkg.add("IntervalRootFinding");

Once the package is installed, it can be imported. Note that you will need also the IntervalArithmetic.jl package.

using IntervalArithmetic, IntervalRootFinding

Introduction

The IntervalRootFinding.jl can be used to rigorously compute the zeros of a function $f: \mathbb{R}^m\rightarrow\mathbb{R}^n$ over a given interval (or interval box) $X\subset\mathbb{R}^m$.

Basic usage

To get started, let's compute the roots of the simple function $f(x)=x^2-2x$ over the interval $[-4, -4]$. To fullfil the task, we can use the function roots from the IntervalRootFinding package which has syntax

rts = roots(f, X, method),
where

• f is the function whose roots we want to find

• X is the interval over which we look for the roots

• method is a third optional parameter that specifies what method is used to compute the roots. The following methods are supported: Newton (default), Krawczyk and Bisection.

f(x) = x^2-2x
X = -4..4

rts = roots(f, X) # using newton method

2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([1.99999, 2.00001], :unique)
 Root([-1.72982e-09, 1.38862e-09], :unique)

As can be noticed, the result of the function is an array of elements of type Root. Each Root-element has two fields: interval, containing the interval contain a zero of the function, and status, a "label" telling whether the interval is guaranteed to contain a single root (:unique) or whether the interval may possibly contain multiple or zero roots (:unknown). It is good to notice that only Newton and Krawczyk methods can make some considerations on the uniqueness of the roots. If Bisection method is used, the status of all roots will automatically be unknown, as the following example shows.

roots(f, X, Bisection)

4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([-2.3921e-08, 3.24651e-08], :unknown)
 Root([1.99999, 2.00001], :unknown)
 Root([1.99999, 2], :unknown)
 Root([2, 2.00001], :unknown)

Note also how the Bisection method may return some extra intervals not containing the roots.

Explicit derivatives

Both Newton and Krawczyk method need to compute the derivative of the function. This is done under the hood using the ForwardDiff.jl package. You can help the solver by explicitly giving the derivative of the function as a second parameter. This is particularly useful if ForwardDiff is not able to compute the derivative.

roots(log, x->1/x, -2..2)

1-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([0.999999, 1.00001], :unique)

Note that this will (generally) speed up the computations a little, but it will not affect the accuracy of the result

using BenchmarkTools
@btime roots(log, x->1/x, -2..2)

@btime roots(log, -2..2)

  4.214 μs (54 allocations: 4.17 KiB)
  4.983 μs (70 allocations: 4.67 KiB)
1-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([0.999999, 1.00001], :unique)

Multidimensional example

The function roots works also for multidimensional functions $\mathbb{R}^m\rightarrow\mathbb{R}^n$, with a few adjustments:

• the function should return a SVector, an efficient data structure for arrays whose size is known at compile time. SVector requires the package StaticArrays.

• the function can take as input an array or a tuple of scalars

• The roots are now searched over an IntervalBox

using StaticArrays
g( (x, y) ) = SVector(sin(x), cos(y))
X = IntervalBox(-3..3, 2)

rts = roots(g, X)

2-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}:
 Root([-1.19115e-21, 1.0588e-22] × [-1.5708, -1.57079], :unique)
 Root([-1.19115e-21, 1.0588e-22] × [1.57079, 1.5708], :unique)

Similarly to 1D, if Newton or Krawczyk method are used, the solver can be helped by explicitly providing the Jacobian, which should be a SMatrix

function dg( (x, y) )
    return SMatrix{2, 2}(cos(x), 0, 0, -sin(y))
end

@btime roots(g, dg, X)
@btime roots(g, X)

  51.299 μs (461 allocations: 25.36 KiB)
  49.899 μs (462 allocations: 25.41 KiB)

Tolerance

By default, the roots are found with a tolerance of 1e-15. This can be changed by giving the solver an extra parameter. A Higher tolerance can significantly speed up the computations. However, as a drawback, more intervals might have an :unknown status. Observe the following example

h(x) = cos(x) * sin(1 / x)

@btime roots(h, 0.05..1)

  36.499 μs (140 allocations: 12.98 KiB)
6-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([0.106103, 0.106104], :unique)
 Root([0.318309, 0.31831], :unique)
 Root([0.0795774, 0.0795775], :unique)
 Root([0.0636619, 0.063662], :unique)
 Root([0.0530516, 0.0530517], :unique)
 Root([0.159154, 0.159155], :unique)

@btime roots(h, 0.05..1, Newton, 1e-2)

  22.299 μs (101 allocations: 10.20 KiB)
6-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([0.317179, 0.319299], :unique)
 Root([0.157209, 0.165989], :unknown)
 Root([0.104739, 0.107542], :unknown)
 Root([0.0515382, 0.0531049], :unique)
 Root([0.0785458, 0.0797755], :unknown)
 Root([0.0570253, 0.0641614], :unknown)

@btime roots(h, 0.05..1, Newton, 1e-1)

  8.266 μs (46 allocations: 3.92 KiB)
3-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([0.107541, 0.165989], :unknown)
 Root([0.283803, 0.3555], :unknown)
 Root([0.0499999, 0.107542], :unknown)