Skip to content

Commit

Permalink
Merge pull request #74 from jverzani/nopolys2
Browse files Browse the repository at this point in the history
Nopolys2
  • Loading branch information
jverzani authored Jun 22, 2017
2 parents 10a101f + ed57417 commit 33b37b3
Show file tree
Hide file tree
Showing 14 changed files with 93 additions and 1,320 deletions.
1 change: 0 additions & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@ os:
- osx
- linux
julia:
- 0.4
- 0.5
- 0.6
- nightly
Expand Down
89 changes: 11 additions & 78 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
[![Roots](http://pkg.julialang.org/badges/Roots_0.4.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.4)
[![Roots](http://pkg.julialang.org/badges/Roots_0.5.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.5)
[![Roots](http://pkg.julialang.org/badges/Roots_0.5.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.5)
[![Roots](http://pkg.julialang.org/badges/Roots_0.6.svg)](http://pkg.julialang.org/?pkg=Roots&ver=0.6)
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/ChrisRackauckas/roots-jl)
Windows: [![Build status](https://ci.appveyor.com/api/projects/status/goteuptn5kypafyl?svg=true)](https://ci.appveyor.com/project/jverzani/roots-jl)

# Root finding functions for Julia

Expand Down Expand Up @@ -34,36 +34,10 @@ algorithm based on its argument(s):
* `fzeros(f, a::Real, b::Real; no_pts::Int=200)` will split
the interval `[a,b]` into many subintervals and search for zeros in
each using a bracketing method if possible. This naive algorithm
will miss double zeros that lie within the same subinterval.
may miss double zeros that lie within the same subinterval and zeros
where there is no crossing of the x-axis.


For polynomials either of class `Poly` (from the `Polynomials`
package) or from functions which are of polynomial type there are
specializations:

* The `roots` function will dispatch to the `roots` function of the
`Polynomials` package to return all roots (including
complex ones) of the polynomial.


* `fzeros(f::Function)` calls `real_roots` to find the real roots of
the polynomial. For polynomials with integer coefficients, the
rational roots are found first, and then numeric approximations to
the remaining real roots are returned.

* For polynomial functions over the integers or rational numbers, the
`factor` function will return a dictionary of factors (as
`Polynomials`) and their multiplicities.

* For polynomials over the real numbers, the `multroot` function will
return the roots and their multiplicities through a dictionary. The
`roots` function from the `Polynomials` package will find all the
roots of a polynomial. Its performance degrades when the polynomial
has high multiplicities. The `multroot` function is provided to
handle this case a bit better. The function follows algorithms due
to Zeng,
["Computing multiple roots of inexact polynomials", Math. Comp. 74 (2005), 869-903](http://www.ams.org/journals/mcom/2005-74-250/S0025-5718-04-01692-8/home.html).




Expand Down Expand Up @@ -108,53 +82,6 @@ fzero(x^5 - x - 1, 1.0)



All real roots of a polynomial can be found at once:

```
f(x) = x^5 - x - 1
fzeros(f)
```

Or using an explicit polynomial:

```
using Polynomials
x = poly([0])
fzeros(x^5 -x - 1)
fzeros(x*(x-1)*(x-2)*(x^2 + x + 1))
```

The `factor` command will factor a polynomial or polynomial function with integer or rational
coefficients over the integers:

```
factor(x -> (2x-1)^2 * (4x-3)^5) # Dict(Poly(1) => 1, Poly(-1 + 2x) =>2, Poly(-3 + 4x) => 5)
```


Polynomial root finding using `multroot` is a bit better when multiple roots are present.

```
x = poly([0.0])
p = (x-1)^2 * (x-3)
U = multroot(p)
collect(keys(U)) # compare to roots(p)
```

Again, a polynomial function may be passed in

```
f(x) = (x-1)*(x-2)^2*(x-3)^3
multroot(f)
```

The `factor` function will factor polynomials with rational and integer coefficients over the integers:

```
factor(f)
```


The well-known methods can be used with or without supplied
derivatives. If not specified, the `ForwardDiff` package is used for
automatic differentiation.
Expand Down Expand Up @@ -189,3 +116,9 @@ fzero(D(m), 0, 1) - median(as) # 0.0
```

Some additional documentation can be read [here](http://nbviewer.ipython.org/url/github.com/JuliaLang/Roots.jl/blob/master/doc/roots.ipynb?create=1).


## Polynomials

Special methods for finding roots of polynomials have been moved to
the `PolynomialZeros` package and its `polyroots(f, domain)` function.
4 changes: 1 addition & 3 deletions REQUIRE
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
julia 0.4
julia 0.5
ForwardDiff 0.2.0
Polynomials 0.0.5
PolynomialFactors 0.0.3
Compat 0.17.0
2 changes: 0 additions & 2 deletions appveyor.yml
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
environment:
matrix:
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.4/julia-0.4-latest-win32.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.4/julia-0.4-latest-win64.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5-latest-win32.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5-latest-win64.exe"
- JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.6/julia-0.6-latest-win32.exe"
Expand Down
61 changes: 21 additions & 40 deletions doc/roots.ipynb

Large diffs are not rendered by default.

98 changes: 16 additions & 82 deletions doc/roots.md
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ To illustrate, consider the function $f(x) = \cos(x) - x$. From the
graph we see readily that $[0,1]$ is a bracket (which we emphasize
with an overlay):

```
```figure
f(x) = cos(x) - x
plot(f, -2, 2)
plot!([0,1], [0,0], linewidth=2)
Expand Down Expand Up @@ -78,6 +78,12 @@ The basic algorithm used for bracketing when the values are simple
floating point values is the bisection method. For big float values, an
algorithm due to Alefeld, Potra, and Shi is used.

The endpoints can even be infinite:

```
fzero(x -> Inf*sign(x), -Inf, Inf) # Float64 only
```

## Using an initial guess

Bracketing methods have guaranteed convergence, but in general require
Expand All @@ -87,9 +93,10 @@ interface to some different iterative algorithms that are more
efficient. Unlike bracketing methods, these algorithms may not
converge to the desired root if the initial guess is not well chosen.

The basic algorithm is modeled after an algorithm used for
The default algorithm is modeled after an algorithm used for
[HP-34 calculators](http://www.hpl.hp.com/hpjournal/pdfs/IssuePDFs/1979-12.pdf). This
algorithm is more forgiving of the quality of the initial guess. In
algorithm is designed to be more forgiving of the quality of the
initial guess at the cost of possible performing many more steps. In
many cases it satisfies the criteria for a bracketing solution, as it
will use bracketing if within the algorithm a bracket is identified.

Expand All @@ -114,15 +121,15 @@ For even more precision, `BigFloat` numbers can be used

```
x = fzero(sin, big(3))
x, f(x), x - pi
x, sin(x), x - pi
```

### Higher order methods

The default call to `fzero` uses a first order method and
then possibly bracketing, which involves potentially many more function
calls. Though specifying a initial value is more convenient than a
bracket, there may be times where a more efficient algorithm is sough.
bracket, there may be times where a more efficient algorithm is sought.
For such, a higher-order method might be better
suited. There are algorithms of order 1 (secant method), 2
([Steffensen](http://en.wikipedia.org/wiki/Steffensen's_method)), 5,
Expand Down Expand Up @@ -211,7 +218,7 @@ defined below comes from a [test
suite](http://people.sc.fsu.edu/~jburkardt/cpp_src/test_zero/test_zero.html)
of difficult functions. The default method finds the zero starting at 0:

```
```figure
f(x) = cos(100*x)-4*erf(30*x-10)
plot(f, -2, 2)
```
Expand Down Expand Up @@ -252,7 +259,7 @@ fzero(f, x0, order=2)
A graph shows the issue. We have overlayed a 15 steps of Newton's
method, the other algorithms being somewhat similar:

```nocode
```figure,nocode
xs = [x0]
n = 15
for i in 1:(n-1) push!(xs, xs[end] - f(xs[end])/D(f)(xs[end])) end
Expand Down Expand Up @@ -376,7 +383,7 @@ end
To see the trajectory if shot at 45 degrees, we have:


```
```figure
theta = 45
plot(x -> flight(x, theta), 0, howfar(theta))
```
Expand All @@ -396,81 +403,8 @@ tstar = fzero(howfarp, 45)

This shows the differences in the trajectories:

```
```figure
plot(x -> flight(x, tstar), 0, howfar(tstar))
plot!(x -> flight(x, 45), 0, howfar(45))
```

## Polynomials

The `Polynomials` package provides a type for working with polynomial
expressions. In that package, the `roots` function is used to find the
roots of a polynomial expression.


For example,

```
using Polynomials
x = poly([0.0]) # (x - 0.0)
roots((x-1)*(x-2)*(x-3))
```

The `Roots` package provides a few interfaces for finding roots of
polynomial functions.

As a convenience, this package adds a function interface to `roots`:

```
f(x) = (x-1)*(x-2)*(x-3)
roots(f)
```


The `fzeros` function will find the real roots of a univariate polynomial:

```
fzeros(x -> (x-1)*(x-2)*(x^2 + x + 1))
```

(For polynomial functions, no search interval is needed, as it is
for other functions.)

The algorithm can have numeric issues when the polynomial degree gets
too large, or the roots are too close together.


The `polyfactor` function provides a related task for
polynomials with integer or rational coefficients:

```
polyfactor(x -> (x-1)^2*(x-2)^3*(x-3)^4)
```

The linear factors correspond to the *rational* roots. (For such
polynomial functions, `fzeros` will first look for rational roots and
return them as rational numbers before searching for approximate
values for the non-rational roots.)

The `multroot` function will also find the roots. The algorithm can do a
better job when there are multiple roots, as it first identifies the multiplicity structure of the
roots, and then tries to improve these values.


```
multroot((x-1)*(x-2)*(x-3)) # roots, multiplicity
```
The `roots` function degrades as there are multiplicities:

```
p = (x-1)^2*(x-2)^3*(x-3)^4
roots(p)
```

Whereas, `multroot` gets it right.

```
multroot(p)
```

The difference grows as the multiplicities get large.
Loading

0 comments on commit 33b37b3

Please sign in to comment.