Skip to content

Commit

Permalink
Merge e195e51 into 125ecfe
Browse files Browse the repository at this point in the history
  • Loading branch information
jverzani committed Aug 6, 2018
2 parents 125ecfe + e195e51 commit 2a8c51c
Show file tree
Hide file tree
Showing 21 changed files with 1,346 additions and 990 deletions.
96 changes: 54 additions & 42 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -9,33 +9,32 @@ Windows: [![Build status](https://ci.appveyor.com/api/projects/status/goteuptn5k
This package contains simple routines for finding roots of continuous
scalar functions of a single real variable. The `find_zero`function provides the
primary interface. It supports various algorithms through the
specification of an method. These include:
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 with a guaranteed
convergence. For most floating point number types, bisection occurs
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.
others, an algorithm of Alefeld, Potra, and Shi is used. These
methods are guaranteed to converge.

For typically faster convergence -- though not guaranteed -- the
`FalsePosition` method can be specified. This method has one of 12
implementations for a modified secant method to
implementations for a modified secant method to potentially
accelerate convergence.

* 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 generally takes many more
function calls. The higher order methods promise higer order
convergence, though don't always yield results with fewer function
calls than `Order1` or `Order2`.
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`.

* There are two historic methods that require a derivative:
`Roots.Newton` and `Roots.Halley`. (Neither is currently exported.)
If a derivative is not given, an automatic derivative is found using
the `ForwardDiff` package.

* There are two historic methods that require a derivative or two:
`Roots.Newton` and `Roots.Halley`. (Neither is exported.)

Each method's documentation has additional detail.

Expand All @@ -56,7 +55,12 @@ julia> find_zero(f, (-10, 0)) # Bisection if x is a tuple and no method

julia> find_zero(f, (-10, 0), FalsePosition()) # just 11 function evaluations
-0.8155534188089607
```

For non-bracketing methods, the initial position is passed in as a
scalar:

```julia
## find_zero(f, x0::Number) will use Order0()
julia> find_zero(f, 3) # default is Order0()
1.4296118247255556
Expand All @@ -82,7 +86,7 @@ Or,
```julia
using Polynomials
x = variable(Int)
fzero(x^5 - x - 1, 1.0) # 1.1673039782614185
find_zero(x^5 - x - 1, 1.0) # 1.1673039782614185
```

The function should respect the units of the `Unitful` package:
Expand All @@ -97,12 +101,22 @@ y(t) = -g*t^2 + v0*t + y0
find_zero(y, 1s) # 1.886053370668014 s
```

Newton's method can be used without taking derivatives:
Newton's method can be used without taking derivatives, if the
`ForwardDiff` package is available:



```julia
using ForwardDiff
D(f) = x -> ForwardDiff.derivative(f,float(x))
```

Now we have:

```
f(x) = x^3 - 2x - 5
x0 = 2
find_zero(f, x0, Roots.Newton()) # 2.0945514815423265
find_zero((f,D(f)) x0, Roots.Newton()) # 2.0945514815423265
```

Automatic derivatives allow for easy solutions to finding critical
Expand All @@ -127,10 +141,10 @@ fzero(D(m), 0, 1) - median(as) # 0.0
### Multiple zeros

The `find_zeros` function can be used to search for all zeros in a
specified interval. The basic algorithm splits the interval into many
subintervals. For each, if there is a bracket a bracketing algorithm
specified interval. The basic algorithm essentially splits the interval into many
subintervals. For each, if there is a bracket, a bracketing algorithm
is used to identify a zero, otherwise a derivative free method is used
to check. This algorithm can miss zeros for various reasons, so the
to search for zeros. This algorithm can miss zeros for various reasons, so the
results should be confirmed by other means.

```julia
Expand All @@ -141,32 +155,32 @@ find_zeros(f, -10, 10)

### Convergence

For most algorithms (besides the `Bisection` ones) convergence is decided when
For most algorithms, convergence is decided when

* The value f(x_n) ≈ 0 with tolerances `atol` and `rtol` *or*
* The value |f(x_n)| < tol with `tol = max(atol, abs(x_n)*rtol)`, or

* the values x_n ≈ x_{n-1} with tolerances `xatol` and `xrtol` *and*
f(x_n) ≈ 0 with a *relaxed* tolerance based on `atol` and `rtol`.
f(x_n) ≈ 0 with a *relaxed* tolerance based on `atol` and `rtol`.

The algorithm stops if

* it encounters an `NaN` or an `Inf`, or

* the number of iterations exceed `maxevals`, or

* an algorithm encounters an `NaN` or `Inf` and yet f(x_n) ≈ 0 with a *relaxed* tolerance based on `atol` and `rtol`.
* the number of function calls exceeds `maxfnevals`.

There is no convergence if the number of iterations exceed `maxevals`,
or the number of function calls exceeds `maxfnevals`.
If the algorithm stops and the relaxed convergence criteria is met,
the suspected zero is returned. Otherwise an error is thrown
indicating no convergence. To adjust the tolerances, `find_zero`
accepts keyword arguments `atol`, `rtol`, `xatol`, and `xrtol`.

The tolerances may need to be adjusted. To determine if convergence
occurs due to f(x_n) ≈ 0, it is necessary to consider that even if
`xstar` is the correct answer mathematically, due to floating point
roundoff it is expected that f(xstar) ≈ f'(xstar) ⋅ xstar ⋅ ϵ. The
relative error used accounts for the value of `x`, but the default
tolerance may need adjustment if the derivative is large near the
zero, as the default is a bit aggressive. On the other hand, the
absolute tolerance might seem too relaxed.

To determine if convergence is determined as x_n ≈ x_{n-1} the check
on f(x_n) ≈ 0 is done as algorithms can be fooled by asymptotes, or
other areas where the tangent lines have large slopes.
The `Bisection` and `Roots.A42` methods are guaranteed to converge
even if the tolerances are set to zero, so these are the
defaults. Non-zero values for `xatol` and `xrtol` can be specified to
reduce the number of function calls when lower precision is required.

The `Bisection` and `Roots.A42` methods will converge, so the tolerances are ignored.

## An alternate interface

Expand All @@ -182,10 +196,8 @@ function. `Roots` also provides this alternative interface:
derivative-free method. with the order specified matching one of
`Order0`, `Order1`, etc.

* `fzeros(f, a::Real, b::Real; no_pts::Int=200)` will call `find_zeros`.
* `fzeros(f, a::Real, b::Real)` will call `find_zeros`.

* The function `secant_method`, `newton`, and `halley` provide direct
access to those methods.


## Usage examples
Expand Down
1 change: 0 additions & 1 deletion REQUIRE
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
julia 0.6.0
ForwardDiff 0.6.0
Compat 0.59.0
Missings 0.2.0
35 changes: 17 additions & 18 deletions appveyor.yml
Original file line number Diff line number Diff line change
@@ -1,13 +1,19 @@
environment:
matrix:
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.6/julia-0.6-latest-win64.exe"
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
- julia_version: 0.6
- julia_version: 0.7
- julia_version: latest

platform:
- x86 # 32-bit
- x64 # 64-bit

## uncomment the following lines to allow failures on nightly julia
## (tests will run but not make your overall status red)
matrix:
allow_failures:
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x86/julia-latest-win32.exe"
- JULIA_URL: "https://julialangnightlies-s3.julialang.org/bin/winnt/x64/julia-latest-win64.exe"
- julia_version: latest


branches:
only:
Expand All @@ -21,19 +27,12 @@ notifications:
on_build_status_changed: false

install:
- ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12"
# Download most recent Julia Windows binary
- ps: (new-object net.webclient).DownloadFile(
$env:JULIA_URL,
"C:\projects\julia-binary.exe")
# Run installer silently, output to C:\projects\julia
- C:\projects\julia-binary.exe /S /D=C:\projects\julia
- ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1"))

build_script:
# Need to convert from shallow to complete for Pkg.clone to work
- IF EXIST .git\shallow (git fetch --unshallow)
- C:\projects\julia\bin\julia -e "versioninfo();
Pkg.clone(pwd(), \"Roots\"); Pkg.build(\"Roots\")"
- echo "%JL_BUILD_SCRIPT%"
- C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%"

test_script:
- C:\projects\julia\bin\julia -e "Pkg.test(\"Roots\")"
- echo "%JL_TEST_SCRIPT%"
- C:\julia\bin\julia -e "%JL_TEST_SCRIPT%"
198 changes: 116 additions & 82 deletions doc/roots.ipynb

Large diffs are not rendered by default.

0 comments on commit 2a8c51c

Please sign in to comment.