|
4 | 4 |
|
5 | 5 | [](http://codecov.io/github/JuliaIntervals/IntervalRootFinding.jl?branch=master)
|
6 | 6 |
|
7 |
| -Guaranteed root-finding methods using interval arithmetic from the [`IntervalArithmetic.jl`](https://github.com/JuliaIntervals/IntervalArithmetic.jl) package. |
| 7 | +This package provides guaranteed methods for finding **roots** of functions $f: \mathbb{R}^n \to \mathbb{R}^n$ with $n \ge 1$, i.e. vectors (or scalars, for $n=1$) $\mathbb{x}$ for which $f(\mathbb{x}) = \mathbb{0}$. It guarantees to find *all* roots inside a given box in $\mathbb{R}^n$, or report subboxes for which it is unable to provide guarantees. |
8 | 8 |
|
9 |
| -Basic documentation is available [here](https://juliaintervals.github.io/IntervalRootFinding.jl/latest/). |
| 9 | +To do so, it uses methods from interval analysis, using interval arithmetic from the [`IntervalArithmetic.jl`](https://github.com/JuliaIntervals/IntervalArithmetic.jl) package by the same authors. |
| 10 | + |
| 11 | + |
| 12 | +## Basic usage examples |
| 13 | + |
| 14 | +The basic function is `roots`. A standard Julia function and a box (interval in 1D) are supplied. Optional search methods (currently `Bisection` or `Newton`) and tolerances may be provided. |
| 15 | + |
| 16 | +### 1D |
| 17 | + |
| 18 | +```jl |
| 19 | +julia> using IntervalArithmetic, IntervalRootFinding |
| 20 | + |
| 21 | +julia> rts = roots(x->x^2 - 2, -10..10, Bisection) |
| 22 | +2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}: |
| 23 | + Root([1.41418, 1.4148], :unknown) |
| 24 | + Root([-1.4148, -1.41418], :unknown) |
| 25 | +``` |
| 26 | +An array is returned that consists of `Root` objects, containing an interval and the status of that interval. Here, `:unknown` indicates that there *may* be a root in the interval. Any region not contained in one of the intervals is guaranteed *not* to contain a root. |
| 27 | + |
| 28 | +```jl |
| 29 | +julia> rts = roots(x->x^2 - 2, -10..10) # default is Newton |
| 30 | +2-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}: |
| 31 | + Root([1.41421, 1.41422], :unique) |
| 32 | + Root([-1.41422, -1.41421], :unique) |
| 33 | +``` |
| 34 | + |
| 35 | +The interval Newton method used here can *guarantee* that there exists a unique root in each of these intervals. Again, other regions have been excluded. |
| 36 | + |
| 37 | +Interval methods are not able to control multiple roots: |
| 38 | + |
| 39 | +```jl |
| 40 | +julia> g(x) = (x^2-2)^2 * (x^2 - 3) |
| 41 | +g (generic function with 1 method) |
| 42 | + |
| 43 | +julia> roots(g, -10..10) |
| 44 | +4-element Array{IntervalRootFinding.Root{IntervalArithmetic.Interval{Float64}},1}: |
| 45 | + Root([1.73205, 1.73206], :unique) |
| 46 | + Root([1.41418, 1.4148], :unknown) |
| 47 | + Root([-1.4148, -1.41418], :unknown) |
| 48 | + Root([-1.73206, -1.73205], :unique) |
| 49 | + ``` |
| 50 | + |
| 51 | + The two double roots are reported as being possible roots, but no guarantees are given. The single roots are guaranteed to exist and be unique within the corresponding intervals. |
| 52 | + |
| 53 | + |
| 54 | +### nD |
| 55 | + |
| 56 | +For dimensions $n>1$, functions must return a Julia vector or an `SVector` from the `StaticArrays.jl` package. |
| 57 | + |
| 58 | +The Rosenbrock function is well known to be difficult to optimize, but is not a problem for this package: |
| 59 | + |
| 60 | +```jl |
| 61 | +julia> using StaticArrays; |
| 62 | + |
| 63 | +julia> rosenbrock(xx) = ( (x, y) = xx; SVector( 1 - x, 100 * (y - x^2) ) ); |
| 64 | + |
| 65 | +julia> X = IntervalBox(-1e5..1e5, 2) # 2D IntervalBox; |
| 66 | + |
| 67 | +julia> rts = roots(rosenbrock, X) |
| 68 | +1-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}: |
| 69 | + Root([1, 1] × [1, 1], :unique) |
| 70 | + ``` |
| 71 | + |
| 72 | + Again, a unique root has been found. |
| 73 | + |
| 74 | + |
| 75 | +A 3D example: |
| 76 | +```jl |
| 77 | +julia> function g(x) |
| 78 | + (x1, x2, x3) = x |
| 79 | + SVector( x1^2 + x2^2 + x3^2 - 1, |
| 80 | + x1^2 + x3^2 - 0.25, |
| 81 | + x1^2 + x2^2 - 4x3 |
| 82 | + ) |
| 83 | + end |
| 84 | +g (generic function with 1 method) |
| 85 | + |
| 86 | +julia> X = (-5..5) |
| 87 | +[-5, 5] |
| 88 | + |
| 89 | +julia> @time rts = roots(g, X × X × X) |
| 90 | + 0.843263 seconds (766.37 k allocations: 36.338 MiB, 1.12% gc time) |
| 91 | +4-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{3,Float64}},1}: |
| 92 | + Root([0.440762, 0.440763] × [0.866025, 0.866026] × [0.236067, 0.236068], :unique) |
| 93 | + Root([0.440762, 0.440763] × [-0.866026, -0.866025] × [0.236067, 0.236068], :unique) |
| 94 | + Root([-0.440763, -0.440762] × [0.866025, 0.866026] × [0.236067, 0.236068], :unique) |
| 95 | + Root([-0.440763, -0.440762] × [-0.866026, -0.866025] × [0.236067, 0.236068], :unique) |
| 96 | + ``` |
| 97 | + |
| 98 | + There are guaranteed to be four unique roots. |
| 99 | + |
| 100 | +### Stationary points |
| 101 | + |
| 102 | +Stationary points of a function $f:\mathbb{R}^n \to \mathbb{R}$ may be found as zeros of the gradient. |
| 103 | +The package exports the `∇` operator to calculate gradients using `ForwardDiff.jl`: |
| 104 | + |
| 105 | +```jl |
| 106 | +julia> f(xx) = ( (x, y) = xx; sin(x) * sin(y) ) |
| 107 | +f (generic function with 1 method) |
| 108 | + |
| 109 | +julia> ∇f = ∇(f) # gradient operator from the package |
| 110 | +(::#53) (generic function with 1 method) |
| 111 | + |
| 112 | +julia> rts = roots(∇f, IntervalBox(-5..6, 2), Newton, 1e-5) |
| 113 | +25-element Array{IntervalRootFinding.Root{IntervalArithmetic.IntervalBox{2,Float64}},1}: |
| 114 | + Root([4.71238, 4.71239] × [4.71238, 4.71239], :unique) |
| 115 | + Root([4.71238, 4.71239] × [1.57079, 1.5708], :unique) |
| 116 | + ⋮ |
| 117 | + [output snipped for brevity] |
| 118 | +``` |
| 119 | +
|
| 120 | +
|
| 121 | +
|
| 122 | +Some basic documentation (mostly now out of date) is available at [here](https://juliaintervals.github.io/IntervalRootFinding.jl/latest/). |
10 | 123 |
|
11 | 124 | ## Authors
|
12 | 125 | - [Luis Benet](http://www.cicc.unam.mx/~benet/), Instituto de Ciencias Físicas,
|
|
0 commit comments