Docs (#184)

* use documenter files for documentation

* edit docs; tweak travis.yml

* edits test; chane doc link

* remove problematic test

.travis.yml

@@ -2,20 +2,39 @@ language: julia
 os:
   - osx
   - linux
+  - windows
+
 julia:
   - 1.0
   - 1.1
   - 1.2
   - 1.3
   - 1.4
   - nightly
-matrix:
-  allow_failures:
-    - julia: nightly
+
 notifications:
   email: false
-#script:
-#  - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
-#  - julia -e 'Pkg.clone(pwd()); Pkg.build("Roots"); Pkg.test("Roots"; coverage=true)';
+
+codecov: true
+
+jobs:
+  allow_failures:
+    - julia: nightly
+  include:
+    - stage: "Documentation"
+      julia: 1.4
+      os: linux
+      env:
+        - GKSwstype="100"
+      script:
+        - julia --project=docs/ -e 'using Pkg; Pkg.instantiate(); Pkg.develop(PackageSpec(path=pwd()))'
+        - julia --project=docs/ docs/make.jl
+      after_success: skip
+
 after_success:
-  - julia -e 'cd(Pkg.dir("Roots")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())';
+    - |
+      julia -e '
+        using Pkg
+        Pkg.add("Coverage")
+        using Coverage
+        Codecov.submit(process_folder())'

Project.toml

@@ -1,6 +1,6 @@
 name = "Roots"
 uuid = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-version = "1.0.1"
+version = "1.0.2"
 
 [deps]
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

README.md

@@ -1,6 +1,10 @@
+[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://juliahub.com/docs/Roots/)
 Linux: [![Build Status](https://travis-ci.org/JuliaMath/Roots.jl.svg?branch=master)](https://travis-ci.org/JuliaMath/Roots.jl)
 Windows: [![Build status](https://ci.appveyor.com/api/projects/status/goteuptn5kypafyl?svg=true)](https://ci.appveyor.com/project/jverzani/roots-jl)
 
+https://juliahub.com/docs/Roots/
+Documentation
+
 # Root finding functions for Julia
 
 

docs/.gitignore

@@ -0,0 +1,2 @@
+build/
+site/

docs/Project.toml

@@ -0,0 +1,4 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

docs/make.jl

@@ -0,0 +1,15 @@
+using Documenter
+using Roots
+
+makedocs(
+    sitename = "Roots",
+    format = Documenter.HTML(),
+    modules = [Roots]
+)
+
+# Documenter can also automatically deploy documentation to gh-pages.
+# See "Hosting Documentation" and deploydocs() in the Documenter manual
+# for more information.
+#=deploydocs(
+    repo = "<repository url>"
+)=#

docs/src/index.md

@@ -0,0 +1,83 @@
+# Roots.jl
+
+Documentation for [Roots.jl](https://github.com/JuliaMath/Roots.jl)
+
+
+## About
+
+`Roots` is  a `Julia` package  for finding zeros of continuous
+scalar functions of a single real variable. That  is solving $f(x)=0$ for $x$.
+The `find_zero`function provides the
+primary interface. It supports various algorithms through the
+specification of a method. These include:
+
+* Bisection-like algorithms. For functions where a bracketing interval
+  is known (one where $f(a)$ and $f(b)$ have alternate signs), the
+  `Bisection` method can be specified. For most floating point number
+  types, bisection occurs in a manner exploiting floating point
+  storage conventions. For others, an algorithm of Alefeld, Potra, and
+  Shi is used. These methods are guaranteed to converge.
+
+
+* Several derivative-free methods are implemented. These are specified
+  through the methods `Order0`, `Order1` (the secant method), `Order2`
+  (the Steffensen method), `Order5`, `Order8`, and `Order16`. The
+  number indicates, roughly, the order of convergence. The `Order0`
+  method is the default, and the most robust, but may take many more
+  function calls to converge. The higher order methods promise higher
+  order (faster) convergence, though don't always yield results with
+  fewer function calls than `Order1` or `Order2`. The methods
+  `Roots.Order1B` and `Roots.Order2B` are superlinear and quadratically converging
+  methods independent of the multiplicity of the zero.
+
+
+* There are historic methods that require a derivative or two:
+  `Roots.Newton` and `Roots.Halley`.  `Roots.Schroder` provides a
+  quadratic method, like Newton's method, which is independent of the
+  multiplicity of the zero.
+
+
+
+## Basic usage
+
+Consider  the polynomial   function  $f(x) = x^5 - x + 1/2$. As a polynomial,  its roots, or  zeros, could  be identified with the  `roots` function of  the `Polynomials` package. However, even  that function uses a numeric method to identify   the values, as no  solution with radicals is available. That is, even for polynomials, non-linear root finders are needed to solve $f(x)=0$. 
+
+The `Roots` package provides a variety algorithms for this  task. In this overview, only the  default ones  are illustrated.
+
+For  the function $f(x) = x^5 - x + 1/2$ a simple plot will show a zero  somewhere between $-1.2$ and $-1.0$ and two zeros near $0.6$. 
+
+For the zero between two values at which the function changes sign, a
+bracketing method is useful, as bracketing methods are guaranteed to
+converge for continuous functions by the intermediate value
+theorem. A bracketing algorithm will be used when the initial data is
+passed as a tuple:
+
+```julia
+julia> using Roots
+
+julia> f(x) =  x^5 - x + 1/2
+f (generic function with 1 method)
+
+julia> find_zero(f, (-1.2,  -1))
+-1.0983313019186334
+```
+
+The default algorithm is guaranteed to have an  answer nearly as accurate as is  possible  given the limitations of floating point  computations. 
+
+For the zeros "near" a point,  a non-bracketing method is often used, as generally  the algorithms are more efficient and can be  used in cases where a zero does  not. Passing just  the initial point will dispatch to  such a method:
+
+```julia
+julia> find_zero(f,  0.6)
+0.550606579334135
+```
+
+
+This finds  the answer  to the left of the starting point. To get the other nearby zero, a starting point closer to the answer cana be used.  However,  an initial graph might convince one  that any of the upto 5 reaal  roots  will   occur between `-5`  and `5`.  The `find_zeros` function uses  heuristics and a few of the  algorithms to  identify  all zeros between the specified range. Here  we see  there  are 3:
+
+```julia
+julia> find_zeros(f, -5,  5)
+3-element Array{Float64,1}:
+ -1.0983313019186334
+  0.550606579334135
+  0.7690997031778959
+```

docs/src/reference.md

@@ -0,0 +1,220 @@
+# Reference/API
+
+The `Roots`  package provides several  different  algorithms  to solve `f(x)=0`.
+
+
+```@index
+Pages = ["reference.md"]
+```
+
+```@setup reference
+using Roots
+```
+
+```@meta
+DocTestSetup = quote
+  using Roots
+end
+```
+
+```@meta
+CurrentModule = Roots
+```
+
+##  The `find_zero`  and `find_zeros` functions
+
+There are  two main  functions:  `find_zero`   to  identify  a  zero  of  $f$  given  some initial starting  value or  bracketing i nterval and  `find_zeros` to heuristically identify  all  zeros in  a specified interval.
+
+
+```@docs
+find_zero
+find_zeros
+```
+
+
+##  Classical  methods  based on derivatives
+
+We begin  by  describing  the classical methods even though they are not necessarily  recommended  because they require more work of the  user,  as they give insight into  why there  are a variety  of methods available.
+
+The classical  methods of [Newton](https://en.wikipedia.org/wiki/Newton%27s_method) and  [Halley](https://en.wikipedia.org/wiki/Halley%27s_method) utilize information about  the function  and  its derivative(s) in  an  interative manner  to converge to  a zero of  `f` given an initial starting value. 
+
+Newton's method is   easily described:
+
+From  an initial point,  the  next  point  in  the iterative algorithm is found by identifying the  intersection of  the $x$    axis  with  the tangent line of `f` at the initial  point. This is repeated until convergence  or the realization that   convergence won't happen for the  initial point. Mathematically,  
+
+$x_{n+1}  =  x_{n}  - f(x_n)/f'(x_n).$ 
+
+Some facts  are helpful  to  understand the different   methods  available in `Roots`:
+
+* For Newton's method  there  is a formula  for the error:  Set  $\epsilon_n  = \alpha - x_n$, where  $\alpha$ is the zero, then 
+
+$\epsilon_{n+1} = -f''(\xi_n)/(2f'(\xi_n) \cdot \epsilon_n^2,$ 
+
+here $\xi_n$ is some value  between  $\alpha$ and  $x_n$.
+
+* The error  term, when of the form $|\epsilon_{n+1}| \leq C\cdot|\epsilon_n|^2$, can be used to identify an  interval around $\alpha$ for  which convergence  is  guaranteed. Such convergence is  termed *quadratic*  (order  2).  For  floating point solutions, quadratic convergence and a well chosen initial point can lead to convergence in 4 or 5 iterations. In  general, convergence is termed order  $q$  when $|\epsilon_{n+1}| \approx C\cdot|\epsilon_n|^q$
+
+*  The term $-f''(\xi_n)/(2f'(\xi_n)$ indicates possible issues  when $f''$  is  too big  near $\alpha$  or  $f'$ is too small  near $\alpha$. In particular if $f'(\alpha)  =  0$, there need  not be quadratic  convergence, and convergence  can   take many  iterations. A  zero   for which $f(x) = (x-\alpha)^{1+\beta}\cdot g(x)$, with $g(\alpha) \neq 0$  is called *simple* when $\beta=0$ and  non-simple when  $\beta >  0$. Newton's method is quadratic near *simple  zeros* and need not be quadratic  near  *non-simple* zeros.
+As well,  if  $f''$ is too  big near $\alpha$, or  $f'$ too small near  $\alpha$, or $x_n$  too  far  from  $\alpha$ (that is,  $|\epsilon_n|>1$) the  error  might actually increase and convergence is not guaranteed.
+
+* The explicit form of  the error function can  be used to guarantee convergence for functions with a certain shape (monotonic, convex functions where the sign of $f''$ and $f'$ don't change). Quadratic convergence may only occur once the algorithm is near the zero.
+
+* The number of function evaluations  per step for Newton's method is 2. 
+
+----
+
+```@docs
+Roots.Newton
+Roots.Halley
+```
+
+For non-simple zeros, Schroder showed an additional  derivative can  be used to yield quadratic convergence.
+
+```@docs
+Roots.Schroder
+```
+
+## Derivative free methods
+
+The [secant](https://en.wikipedia.org/wiki/Secant_method) method replaces the  derivative term in Newton's method with the slope of the secant line.  
+
+$x_{n+1} = x_n - (\frac{f(x_n)-f(x_{n-1})}{x_n - x_{n-1}})^{-1}\cdot  f(x_n).$
+
+Though the secant  method   has  convergence  rate of  order $\approx 1.618$ and is not quadratic,  it
+only requires one new  function call per  step  so  can be very effective. Often  function evaluations are the  slowest part of  the computation and, as  well, no derivative is  needed. Because  it  can be  very efficient, the secant  method  is used in  the default method  of `find_zero` when  called with a single initial starting point.
+
+[Steffensen's](https://en.wikipedia.org/wiki/Steffensen%27s_method) method is a quadratially  converging derivative free method  which uses a secant  line  based on $x_n$ and $x_n + f(x_n)$.  Though of  higher  order, it requires  additional function calls per step and depends on a  good initial starting value. Other  derivative free methods are available, trading off  increased function calls for higher-order convergence. They may be  of interest when arbitrary  precision is needed. A  measure of efficiency is $q^{1/r}$ where $q$ is the order of convergence and $r$ the number of function calls per step.   With this measure, the secant method  would be $\approx (1.618)^{1/1}$ and Steffensen's  would be less ($2^{1/2}$).
+
+----
+
+```@docs
+Order1
+Order2
+Order5
+Order8
+Order16
+```
+
+As with  Newton's method, convergence rates are for simple  zeros. For  non simple zeros there are related methods with varying convergence rates:
+
+```@docs
+Roots.King
+Roots.Esser
+```
+
+
+
+## Bracketing methods
+
+The [bisection](https://en.wikipedia.org/wiki/Bisection_method) identifies a zero of a *continuous* function bewteen $a$ and $b$  when  $f(a)$ and  $f(b)$ have different  signs. (The interval $[a,b]$ is called a bracketing interval when $f(a)\cdot  f(b)  <0$.)  The basaic  alorithm is particularly simple, an interval  $[a_i,b_i]$ is  split  at  $c =  (a_i+b_i)/2$. Either  $f(c)=0$,  or one  of  $[a_i,c]$  or  $[c,b_i]$ is a bracketing  interval,  which is  called  $[a_{i+1},b_{i+1}]$. From this  description,  we  see  thaat  $[a_i,b_i]$ has length  $2^{-i}$ times the length of $[a_0,b_0]$, so  the intervals will eventually terminate by finding  a zero, $c$,  or converge  to a zero. This convergence is slow (the efficiency  is only $1$, but guaranteed. For  64-bit  floating point  values, a  reinterpretation  of  how the midpoint  ($c$) is found  leads  to convergence  is  no more  than   $64$ iterations. 
+
+In floating point,  by  guaranteed  convergence we have either an exact zero or a bracketing interval  consisting   of  two  adjacent floating point values. When applied to *non*-continuous  functions,  this algorithm  will identify   an exact  zero or  a zero crossing   of the function. (E.g., applied  to  $f(x)=1/x$ it  will  find  $0$.)
+
+The default selection of  midpoint described above includes no information  about the function $f$ beyonds its  sign. Algorithms exploiting  the shape of the function  can be significantly more efficient. The default bracketing method is due to [Alefeld, Potra, and Shi](https://dl.acm.org/doi/10.1145/210089.210111), and has  efficiency $\approx 1.6686$. It  is  also   used in the  default method for `find_zero` when a  single initial starting point is given if a bracketing interval is identified.
+
+----
+
+```@docs
+Bisection
+Roots.A42
+Roots.AlefeldPotraShi
+Roots.Brent
+FalsePosition
+```
+
+## Hybrid  methods
+
+A useful  strategy  is   to  begin with a non-bracketing  method and switch to a bracketing method should a bracket be encoutered. This  allows   for the identification of zeros which are not surrounded by a bracket, and have guaranteed convergence  should a bracket be  encountered.  Is  is  used  by default by `find_zero(f,a)`.
+
+```@docs
+Roots.Order0
+```
+
+## Convergence
+
+Identifying when an algorithm converges or diverges requires specifications of tolerances  and convergence criteria.
+
+In the case of exact bisection, convergence is mathematically
+guaranteed. For floating point numbers, either an *exact* zero is
+found, or the bracketing interval can be subdivided into $[a_n,b_n]$
+with $a_n$ and $b_n$ being adjacent floating point values. That is
+$b_n-a_n$ is as small as possible in floating point numbers. This can
+be considered a stopping criteria in $\Delta x$. For early termination
+(less precision but fewer function calls) a tolerance can be given so
+that if $\Delta_n=b_n-a_n$ is small enough the algorithm stops
+successfully.  In floating point assessing if $b_n \approx a_n$
+requires two tolerances: a *relative* tolerance, as the minimial
+differences in floating point values depend on the size of $b_n$ and
+$a_n$, and an absolute tolerancef or values near $0$. The valuese
+`xrtol` anad `xatol` are passed to the `Base.isapprox` function to
+determine this.
+
+Relying on the closeness of two $x$ values will not be adequate for
+all problems, sa there are examples where the difference
+$\Delta_n=|x_n-x_{n-1}|$ can be quite small, $0$ even, yet $f(x_n)$ is
+not near a $0$. As such, for non-bracketing methods, a check on the
+size of $f(x_n)$ is also used. As we find floating point
+approximations to $\alpha$, the zero, we must consider values small
+when $f(\alpha(1+\epsilon))$ is small. By Taylor's aapproximation, we
+can expect this to be around 
+$\alpha\cdot \epsilon \cdot f'(\alpha)$. 
+That is, small depends on the size of $\alpha$ and the
+derivative at $\alpha$.  The former is handled by both relative and absolute
+tolerances (`rtol` and `atol`).  The size of $f'(\alpha)$ is problem
+dependent, and can be accommodated by larger relative or absolute
+tolerances.
+
+When an algorithm returns  a  `NaN` value,  it terminates. This  can  happen near convergence or  may indicate some issues.  Early termination is checked for convergence  in the  size  of $f(x_{n-1})$ with a relaxed tolerance when `strict=false` is specified (the default).
+
+!!! note "Relative tolerances  and assessing  $f(x) \approx  0$"
+    The use of  relative tolerances  to  check  if   $f(x)  \approx  0$ can lead  to spurious  answers  where  $x$ is very large   (and  hence the relative  tolerance  is large). The return of  very  large solutions  should  be checked against expectations  of the  answer.
+
+
+Deciding if an algorithm won't  terminate   is  done  through  counting the number or  iterations performed or  the number of  function   evaluations.   These are  adjusted through `maxevals` and `maxfnevals`. As  most  algorithms are superlinear, convergence happens rapidly near  the answer, but  all the algorithms  can take  a while  to  get near  an  answer, even   when progress  is made. As  such, the valuese can indicate non-convergence by being too slow.  On the other hand, leaving this too  large  can slow the  determination.
+
+
+Convergence criteria are method dependent and are determined by  the  `Roots.assess_convergence`  methods.
+
+
+```@docs
+Roots.assess_convergence
+```
+
+Default tolerances  are specified through the `Roots.default_tolerances` methods.
+
+```@docs
+Roots.default_tolerances
+```
+
+
+
+
+## Performance considerations
+
+
+The abstractions and many checks for  convergence employed by `find_zero` have a performance cost. When that is a critical concern, there are  several "simple" methods provided which can offer significantly improved performaance.
+
+```@docs
+Roots.secant_method
+Roots.bisection
+Roots.muller
+Roots.newton
+Roots.dfree
+```
+
+## Matlab interface
+
+The initial naming scheme used `fzero` instead  of `fzeros`, following the name of the  MATLAB function [fzero](https://www.mathworks.com/help/matlab/ref/fzero.html). This interface  is not recommended, but, for now, still maintained. 
+
+```@docs
+Roots.fzero
+```
+
+!!! note
+    At one point there  were technical reasons (specialization on functions) to have both `fzero` and `find_zero` interfaces, but that is no longer the case.
+
+
+
+
+
+