Reference/API

find_zero(fs, x0, M, [N::AbstractBracketing]; kwargs...)

Interface to one of several methods for find zeros of a univariate function.

Initial starting value

For most methods, x0 is a scalar value indicating the initial value in the iterative procedure. (Secant methods can have a tuple specify their initial values.) Values must be a subtype of Number and have methods for float, real, and oneunit defined.

For bracketing intervals, x0 is specified as a tuple or a vector. A bracketing interval, (a,b), is one where f(a) and f(b) have different signs.

Specifying a method

A method is specified to indicate which algorithm to employ:

• There are methods for bisection where a bracket is specified: Bisection, Roots.A42, FalsePosition

• There are several derivative-free methods: cf. Order0, Order1 (secant method), Order2 (Steffensen), Order5, Order8, and Order16, where the number indicates the order of the convergence. Methods Roots.Order1B and Roots.Order2B implement methods useful when the desired zero has a multiplicity.

• There are some classical methods where derivatives are required: Roots.Newton, Roots.Halley, Roots.Schroder. (The are not exported.)

For more detail, see the help page for each method (e.g., ?Order1).

If no method is specified, the default method depends on x0:

• If x0 is a scalar, the default is the slower, but more robust Order0 method.

• If x0 is a tuple or vector indicating a bracketing interval, the Bisection method is used. (The exact algorithm depends on the number type, the tolerances, and verbose.)

Specifying the function

The function(s) are passed as the first argument.

For the few methods that use a derivative (Newton, Halley, Shroder, and optionally Order5) a tuple of functions is used.

Optional arguments (tolerances, limit evaluations, tracing)

• xatol - absolute tolerance for x values. Passed to isapprox(x_n, x_{n-1})
• xrtol - relative tolerance for x values. Passed to isapprox(x_n, x_{n-1})
• atol - absolute tolerance for f(x) values.
• rtol - relative tolerance for f(x) values.
• maxevals - limit on maximum number of iterations
• maxfnevals - limit on maximum number of function evaluations
• strict - if false (the default), when the algorithm stops, possible zeros are checked with a relaxed tolerance
• verbose - if true a trace of the algorithm will be shown on successful completion.

See the help string for Roots.assess_convergence for details on convergence. See the help page for Roots.default_tolerances(method) for details on the default tolerances.

In general, with floating point numbers, convergence must be understood as not an absolute statement. Even if mathematically α is an answer and xstar the floating point realization, it may be that f(xstar) - f(α) ≈ xstar ⋅ f'(α) ⋅ eps(α), so tolerances must be appreciated, and at times specified.

For the Bisection methods, convergence is guaranteed, so the tolerances are set to be 0 by default.

If a bracketing method is passed in after the method specification, when a bracket is found, the bracketing method will used to identify the zero. This is what Order0 does by default, with a secant step initially and the A42 method when a bracket is encountered.

Note: The order of the method is hinted at in the naming scheme. A scheme is order r if, with eᵢ = xᵢ - α, eᵢ₊₁ = C⋅eᵢʳ. If the error eᵢ is small enough, then essentially the error will gain r times as many leading zeros each step. However, if the error is not small, this will not be the case. Without good initial guesses, a high order method may still converge slowly, if at all. The OrderN methods have some heuristics employed to ensure a wider range for convergence at the cost of not faithfully implementing the method, though those are available through unexported methods.

Examples:

Default methods.

julia> using Roots

julia> find_zero(sin, 3)  # use Order0()
3.141592653589793

julia> find_zero(sin, (3,4)) # use Bisection()
3.1415926535897936

Specifying a method,

julia> find_zero(sin, (3,4), Order1())            # can specify two starting points for secant method
3.141592653589793

julia> find_zero(sin, 3.0, Order2())              # Use Steffensen method
3.1415926535897936

julia> find_zero(sin, big(3.0), Order16())        # rapid convergence
3.141592653589793238462643383279502884197169399375105820974944592307816406286198

julia> find_zero(sin, (3, 4), Roots.A42()())      # fewer function calls than Bisection(), in this case
ERROR: MethodError: objects of type Roots.A42 are not callable
Stacktrace:
 [1] top-level scope at REPL[12]:1

julia> find_zero(sin, (3, 4), FalsePosition(8))   # 1 of 12 possible algorithms for false position
3.141592653589793

julia> find_zero((sin,cos), 3.0, Roots.Newton())  # use Newton's method
3.141592653589793

julia> find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley())  # use Halley's method
3.141592653589793

Changing tolerances.

julia> fn = x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1);

julia> x0, xstar = 3.0,  2.9947567209477;

julia> find_zero(fn, x0, Order2()) ≈ xstar       
true

julia> find_zero(fn, x0, Order2(), atol=0.0, rtol=0.0) # error: x_n ≉ x_{n-1}; just f(x_n) ≈ 0
ERROR: Roots.ConvergenceFailed("Stopped at: xn = 2.991488255523429")
Stacktrace:
 [1] find_zero(::Function, ::Float64, ::Order2, ::Nothing; tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol,Symbol},NamedTuple{(:atol, :rtol),Tuple{Float64,Float64}}}) at /Users/verzani/julia/Roots/src/find_zero.jl:670
 [2] top-level scope at REPL[12]:1

julia> fn = x -> (sin(x)*cos(x) - x^3 + 1)^9;

julia> x0, xstar = 1.0,  1.112243913023029;

julia> find_zero(fn, x0, Order2()) ≈ xstar           
true

julia> find_zero(fn, x0, Order2(), maxevals=3)    # Roots.ConvergenceFailed: 26 iterations needed, not 3
ERROR: Roots.ConvergenceFailed("Stopped at: xn = 1.0482748172022405")
Stacktrace:
 [1] find_zero(::Function, ::Float64, ::Order2, ::Nothing; tracks::Roots.NullTracks, verbose::Bool, kwargs::Base.Iterators.Pairs{Symbol,Int64,Tuple{Symbol},NamedTuple{(:maxevals,),Tuple{Int64}}}) at /Users/verzani/julia/Roots/src/find_zero.jl:670
 [2] top-level scope at REPL[12]:1

Tracing output.

julia> find_zero(x->sin(x), 3.0, Order2(), verbose=true)   # 3 iterations
Results of univariate zero finding:

* Converged to: 3.1415926535897936
* Algorithm: Order2()
* iterations: 2
* function evaluations: 5
* stopped as |f(x_n)| ≤ max(δ, max(1,|x|)⋅ϵ) using δ = atol, ϵ = rtol

Trace:
x_0 =  3.0000000000000000,	 fx_0 =  0.1411200080598672
x_1 =  3.1425464815525403,	 fx_1 = -0.0009538278181169
x_2 =  3.1415926535897936,	 fx_2 = -0.0000000000000003

3.1415926535897936

julia> find_zero(x->sin(x)^5, 3.0, Order2(), verbose=true) # 22 iterations
Results of univariate zero finding:

* Converged to: 3.140534939851113
* Algorithm: Order2()
* iterations: 22
* function evaluations: 46
* stopped as |f(x_n)| ≤ max(δ, max(1,|x|)⋅ϵ) using δ = atol, ϵ = rtol

Trace:
x_0 =  3.0000000000000000,	 fx_0 =  0.0000559684091879
x_1 =  3.0285315910604353,	 fx_1 =  0.0000182783542076
x_2 =  3.0512479368872216,	 fx_2 =  0.0000059780426727
x_3 =  3.0693685883136541,	 fx_3 =  0.0000019566875137
x_4 =  3.0838393517913989,	 fx_4 =  0.0000006407325974
x_5 =  3.0954031275790856,	 fx_5 =  0.0000002098675747
x_6 =  3.1046476918938040,	 fx_6 =  0.0000000687514870
x_7 =  3.1120400753639790,	 fx_7 =  0.0000000225247921
x_8 =  3.1179523212360416,	 fx_8 =  0.0000000073801574
x_9 =  3.1226812716693950,	 fx_9 =  0.0000000024181703
x_10 =  3.1264639996586729,	 fx_10 =  0.0000000007923528
x_11 =  3.1294899615147704,	 fx_11 =  0.0000000002596312
x_12 =  3.1319106162503414,	 fx_12 =  0.0000000000850746
x_13 =  3.1338470835729701,	 fx_13 =  0.0000000000278769
x_14 =  3.1353962361189192,	 fx_14 =  0.0000000000091346
x_15 =  3.1366355527270868,	 fx_15 =  0.0000000000029932
x_16 =  3.1376269806673180,	 fx_16 =  0.0000000000009808
x_17 =  3.1384199603056921,	 fx_17 =  0.0000000000003215
x_18 =  3.1390543962469541,	 fx_18 =  0.0000000000001054
x_19 =  3.1395625842920500,	 fx_19 =  0.0000000000000345
x_20 =  3.1399667213451634,	 fx_20 =  0.0000000000000114
x_21 =  3.1402867571209034,	 fx_21 =  0.0000000000000038
x_22 =  3.1405349398511131,	 fx_22 =  0.0000000000000013

3.140534939851113

julia> find_zero(x->sin(x)^5, 3.0, Roots.Order2B(), verbose=true) # 2 iterations
Results of univariate zero finding:

* Converged to: 3.1397074174874358
* Algorithm: Roots.Order2B()
* iterations: 2
* function evaluations: 7
* Note: Estimate for multiplicity had issues.
	Algorithm stopped early, but |f(xn)| < ϵ^(1/3), where ϵ depends on xn, rtol, and atol.

Trace:
x_0 =  3.0000000000000000,	 fx_0 =  0.0000559684091879
x_1 =  3.1397074174874358,	 fx_1 =  0.0000000000000238
x_2 =  3.1397074174874358,	 fx_2 =  0.0000000000000238

3.1397074174874358

find_zeros(f, a, b; [no_pts=12, k=8, naive=false, xatol, xrtol, atol, rtol])

Search for zeros of f in the interval [a,b].

Examples

julia> using Roots

julia> find_zeros(x -> exp(x) - x^4, -5, 20)        # a few well-spaced zeros
3-element Array{Float64,1}:
 -0.8155534188089606
  1.4296118247255556
  8.613169456441398

julia> find_zeros(x -> sin(x^2) + cos(x)^2, 0, 2pi)  # many zeros
12-element Array{Float64,1}:
 1.78518032659534
 2.391345462376604
 3.2852368649448853
 3.3625557095737544
 4.016412952618305
 4.325091924521049
 4.68952781386834
 5.00494459113514
 5.35145266881871
 5.552319796014526
 5.974560835055425
 6.039177477770888

julia> find_zeros(x -> cos(x) + cos(2x), 0, 4pi)    # mix of simple, non-simple zeros
6-element Array{Float64,1}:
  1.0471975511965976
  3.141592653589943
  5.235987755982988
  7.330382858376184
  9.424777960769228
 11.519173063162574

julia> f(x) = (x-0.5) * (x-0.5001) * (x-1)          # nearby zeros
f (generic function with 1 method)

julia> find_zeros(f, 0, 2)
3-element Array{Float64,1}:
 0.5
 0.5001
 1.0

julia> f(x) = (x-0.5) * (x-0.5001) * (x-4) * (x-4.001) * (x-4.2)
f (generic function with 1 method)

julia> find_zeros(f, 0, 10)
3-element Array{Float64,1}:
 0.5
 0.5001
 4.2

julia> f(x) = (x-0.5)^2 * (x-0.5001)^3 * (x-4) * (x-4.001) * (x-4.2)^2  # hard to identify
f (generic function with 1 method)

julia> find_zeros(f, 0, 10, no_pts=21)                # too hard for default
5-element Array{Float64,1}:
 0.49999999999999994
 0.5001
 4.0
 4.001
 4.200000000000001

The basic algorithm checks for zeros among the endpoints, and then divides the interval (a,b) into no_pts-1 subintervals and then proceeds to look for zeros through bisection or a derivative-free method. As checking for a bracketing interval is relatively cheap and bisection is guaranteed to converge, each interval has k pairs of intermediate points checked for a bracket.

If any zeros are found, the algorithm uses these to partition (a,b) into subintervals. Each subinterval is shrunk so that the endpoints are not zeros and the process is repeated on the subinterval. If the initial interval is too large, then the naive scanning for zeros may be fruitless and no zeros will be reported. If there are nearby zeros, the shrinking of the interval may jump over them, though as seen in the examples, nearby roots can be identified correctly, though for really nearby points, or very flat functions, it may help to increase no_pts.

The tolerances are used to shrink the intervals, but not to find zeros within a search. For searches, bisection is guaranteed to converge with no specified tolerance. For the derivative free search, a modification of the Order0 method is used, which at worst case compares |f(x)| <= 8*eps(x) to identify a zero. The algorithm might identify more than one value for a zero, due to floating point approximations. If a potential pair of zeros satisfy isapprox(a,b,atol=sqrt(xatol), rtol=sqrt(xrtol)) then they are consolidated.

The algorithm can make many function calls. When zeros are found in an interval, the naive search is carried out on each subinterval. To cut down on function calls, though with some increased chance of missing some zeros, the adaptive nature can be skipped with the argument naive=true or the number of points stepped down.

The algorithm is derived from one in a PR by @djsegal.

Roots.Newton()

Implements Newton's method: xᵢ₊₁ = xᵢ - f(xᵢ)/f'(xᵢ). This is a quadratically convergent method requiring one derivative. Two function calls per step.

Example

find_zero((sin,cos), 3.0, Roots.Newton())

If function evaluations are expensive one can pass in a function which returns (f, f/f') as follows

find_zero(x -> (sin(x), sin(x)/cos(x)), 3.0, Roots.Newton())

This can be advantageous if the derivative is easily computed from the value of f, but otherwise would be expensive to compute.

The error, eᵢ = xᵢ - α, can be expressed as eᵢ₊₁ = f[xᵢ,xᵢ,α]/(2f[xᵢ,xᵢ])eᵢ² (Sidi, Unified treatment of regula falsi, Newton-Raphson, secant, and Steffensen methods for nonlinear equations).

Roots.Halley()

Implements Halley's method, xᵢ₊₁ = xᵢ - (f/f')(xᵢ) * (1 - (f/f')(xᵢ) * (f''/f')(xᵢ) * 1/2)^(-1) This method is cubically converging, but requires more function calls per step (3) than other methods.

Example

find_zero((sin, cos, x->-sin(x)), 3.0, Roots.Halley())

If function evaluations are expensive one can pass in a function which returns (f, f/f',f'/f'') as follows

find_zero(x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x)), 3.0, Roots.Halley())

This can be advantageous if the derivatives are easily computed from the value of f, but otherwise would be expensive to compute.

The error, eᵢ = xᵢ - α, satisfies eᵢ₊₁ ≈ -(2f'⋅f''' -3⋅(f'')²)/(12⋅(f'')²) ⋅ eᵢ³ (all evaluated at α).

Roots.Schroder()

Schröder's method, like Halley's method, utilizes f, f', and f''. Unlike Halley it is quadratically converging, but this is independent of the multiplicity of the zero (cf. Schröder, E. "Über unendlich viele Algorithmen zur Auflösung der Gleichungen." Math. Ann. 2, 317-365, 1870; mathworld). Schröder's method applies Newton's method to f/f', a function with all simple zeros.

Example

m = 2
f(x) = (cos(x)-x)^m
fp(x) = (-x + cos(x))*(-2*sin(x) - 2)
fpp(x) = 2*((x - cos(x))*cos(x) + (sin(x) + 1)^2)
find_zero((f, fp, fpp), pi/4, Roots.Halley())    # 14 steps
find_zero((f, fp, fpp), 1.0, Roots.Schroder())    # 3 steps

(Whereas, when m=1, Halley is 2 steps to Schröder's 3.)

If function evaluations are expensive one can pass in a function which returns (f, f/f',f'/f'') as follows

find_zero(x -> (sin(x), sin(x)/cos(x), -cos(x)/sin(x)), 3.0, Roots.Schroder())

This can be advantageous if the derivatives are easily computed from the value of f, but otherwise would be expensive to compute.

The error, eᵢ = xᵢ - α, is the same as Newton with f replaced by f/f'.

Order1()
Roots.Secant()

The Order1() method is an alias for Secant. It specifies the secant method. This method keeps two values in its state, x_n and x_n1. The updated point is the intersection point of x axis with the secant line formed from the two points. The secant method uses 1 function evaluation per step and has order (1+sqrt(5))/2.

The error, eᵢ = xᵢ - α, satisfies eᵢ₊₂ = f[xᵢ₊₁,xᵢ,α] / f[xᵢ₊₁,xᵢ] * (xᵢ₊₁-α) * (xᵢ - α).

Order2()
Roots.Steffensen()

The quadratically converging Steffensen method is used for the derivative-free Order2() algorithm. Unlike the quadratically converging Newton's method, no derivative is necessary, though like Newton's method, two function calls per step are. Steffensen's algorithm is more sensitive than Newton's method to poor initial guesses when f(x) is large, due to how f'(x) is approximated. This algorithm, Order2, replaces a Steffensen step with a secant step when f(x) is large.

The error, eᵢ - α, satisfies eᵢ₊₁ = f[xᵢ, xᵢ+fᵢ, α] / f[xᵢ,xᵢ+fᵢ] ⋅ (1 - f[xᵢ,α] ⋅ eᵢ²

Order5()

Implements an order 5 algorithm from A New Fifth Order Derivative Free Newton-Type Method for Solving Nonlinear Equations by Manoj Kumar, Akhilesh Kumar Singh, and Akanksha, Appl. Math. Inf. Sci. 9, No. 3, 1507-1513 (2015), DOI: 10.12785/amis/090346. Four function calls per step are needed.

The error, eᵢ = xᵢ - α, satisfies eᵢ₊₁ = K₁ ⋅ K₅ ⋅ M ⋅ eᵢ⁵ + O(eᵢ⁶)

Order8()

Implements an eighth-order algorithm from New Eighth-Order Derivative-Free Methods for Solving Nonlinear Equations by Rajinder Thukral, International Journal of Mathematics and Mathematical Sciences Volume 2012 (2012), Article ID 493456, 12 pages DOI: 10.1155/2012/493456. Four function calls per step are required.

The error, eᵢ = xᵢ - α, is expressed as eᵢ₊₁ = K ⋅ eᵢ⁸ in (2.25) of the paper for an explicit K.

Order16()

Implements the order 16 algorithm from New Sixteenth-Order Derivative-Free Methods for Solving Nonlinear Equations by R. Thukral, American Journal of Computational and Applied Mathematics p-ISSN: 2165-8935; e-ISSN: 2165-8943; 2012; 2(3): 112-118 doi: 10.5923/j.ajcam.20120203.08.

Five function calls per step are required. Though rapidly converging, this method generally isn't faster (fewer function calls/steps) over other methods when using Float64 values, but may be useful for solving over BigFloat.

The error, eᵢ = xᵢ - α, is expressed as eᵢ₊₁ = K⋅eᵢ¹⁶ for an explicit K in equation (50) of the paper.

Roots.Order1B()
Roots.King()

A superlinear (order 1.6...) modification of the secant method for multiple roots. Presented in A SECANT METHOD FOR MULTIPLE ROOTS, by RICHARD F. KING, BIT 17 (1977), 321-328

The basic idea is similar to Schroder's method: apply the secant method to f/f'. However, this uses f' ~ fp = (fx - f(x-fx))/fx (a Steffensen step). In this implementation, Order1B, when fx is too big, a single secant step of f is used.

The asymptotic error, eᵢ = xᵢ - α, is given by eᵢ₊₂ = 1/2⋅G''/G'⋅ eᵢ⋅eᵢ₊₁ + (1/6⋅G'''/G' - (1/2⋅G''/G'))^2⋅eᵢ⋅eᵢ₊₁⋅(eᵢ+eᵢ₊₁).

Roots.Order2B()
Roots.Esser()

Esser's method. This is a quadratically convergent method that, like Schroder's method, does not depend on the multiplicity of the zero. Schroder's method has update step x - r2/(r2-r1) * r1, where ri = f^(i-1)/f^(i). Esser approximates f' ~ f[x-h, x+h], f'' ~ f[x-h,x,x+h], where h = fx, as with Steffensen's method, Requiring 3 function calls per step. The implementation Order2B uses a secant step when |fx| is considered too large.

Esser, H. Computing (1975) 14: 367. https://doi.org/10.1007/BF02253547 Eine stets quadratisch konvergente Modifikation des Steffensen-Verfahrens

Example

f(x) = cos(x) - x
g(x) = f(x)^2
x0 = pi/4
find_zero(f, x0, Order2(), verbose=true)        #  3 steps / 7 function calls
find_zero(f, x0, Roots.Order2B(), verbose=true) #  4 / 9
find_zero(g, x0, Order2(), verbose=true)        #  22 / 45
find_zero(g, x0, Roots.Order2B(), verbose=true) #  4 / 10

Bisection()
Roots.BisectionExact()

If possible, will use the bisection method over Float64 values. The bisection method starts with a bracketing interval [a,b] and splits it into two intervals [a,c] and [c,b], If c is not a zero, then one of these two will be a bracketing interval and the process continues. The computation of c is done by _middle, which reinterprets floating point values as unsigned integers and splits there. It was contributed by Jason Merrill. This method avoids floating point issues and when the tolerances are set to zero (the default) guarantees a "best" solution (one where a zero is found or the bracketing interval is of the type [a, nextfloat(a)]).

When tolerances are given, this algorithm terminates when the midpoint is approximately equal to an endpoint using absolute tolerance xatol and relative tolerance xrtol.

When a zero tolerance is given and the values are not Float64 values, this will call the A42 method.

Roots.A42()

Bracketing method which finds the root of a continuous function within a provided interval [a, b], without requiring derivatives. It is based on algorithm 4.2 described in: G. E. Alefeld, F. A. Potra, and Y. Shi, "Algorithm 748: enclosing zeros of continuous functions," ACM Trans. Math. Softw. 21, 327–344 (1995), DOI: 10.1145/210089.210111. Originally by John Travers.

Roots.AlefeldPotraShi()

Follows algorithm in "ON ENCLOSING SIMPLE ROOTS OF NONLINEAR EQUATIONS", by Alefeld, Potra, Shi; DOI: 10.1090/S0025-5718-1993-1192965-2. Efficiency is 1.618. Less efficient, but can be faster than the A42 method. Originally by John Travers.

Roots.Brent()

An implementation of Brent's (or Brent-Dekker) method. This method uses a choice of inverse quadratic interpolation or a secant step, falling back on bisection if necessary.

FalsePosition()

Use the false position method to find a zero for the function f within the bracketing interval [a,b].

The false position method is a modified bisection method, where the midpoint between [a_k, b_k] is chosen to be the intersection point of the secant line with the x axis, and not the average between the two values.

To speed up convergence for concave functions, this algorithm implements the 12 reduction factors of Galdino (A family of regula falsi root-finding methods). These are specified by number, as in FalsePosition(2) or by one of three names FalsePosition(:pegasus), FalsePosition(:illinois), or FalsePosition(:anderson_bjork) (the default). The default choice has generally better performance than the others, though there are exceptions.

For some problems, the number of function calls can be greater than for the bisection64 method, but generally this algorithm will make fewer function calls.

Examples

find_zero(x -> x^5 - x - 1, (-2, 2), FalsePosition())

Order0()

The Order0 method is engineered to be a more robust, though possibly slower, alternative to the other derivative-free root-finding methods. The implementation roughly follows the algorithm described in Personal Calculator Has Key to Solve Any Equation f(x) = 0, the SOLVE button from the HP-34C. The basic idea is to use a secant step. If along the way a bracket is found, switch to bisection, using AlefeldPotraShi. If the secant step fails to decrease the function value, a quadratic step is used up to 4 times.

This is not really 0-order: the secant method has order 1.6...[Wikipedia](https://en.wikipedia.org/wiki/Secantmethod#Comparisonwithotherroot-findingmethods and the the bracketing method has order 1.6180...Wikipedia so for reasonable starting points, this algorithm should be superlinear, and relatively robust to non-reasonable starting points.

Roots.assess_convergence(method, state, options)

Assess if algorithm has converged.

If alogrithm hasn't converged returns false.

If algorithm has stopped or converged, return true and sets one of state.stopped, state.x_converged, state.f_converged, or state.convergence_failed; as well, a message may be set.

• state.x_converged = true if abs(xn1 - xn0) < max(xatol, max(abs(xn1), abs(xn0)) * xrtol)

• state.f_converged = true if |f(xn1)| < max(atol, |xn1|*rtol)

• state.convergence_failed = true if xn1 or fxn1 is NaN or an infinity

• state.stopped = true if the number of steps exceed maxevals or the number of function calls exceeds maxfnevals.

In find_zero, stopped values (and x_converged) are checked for convergence with a relaxed tolerance.

default_tolerances(M::AbstractUnivariateZeroMethod, [T], [S])

The default tolerances for most methods are xatol=eps(T), xrtol=eps(T), atol=4eps(S), and rtol=4eps(S), with the proper units (absolute tolerances have the units of x and f(x); relative tolerances are unitless). For Complex{T} values, T is used.

The number of iterations is limited by maxevals=40, the number of function evaluations is not capped.

default_tolerances(M::Bisection, [T], [S])

For Bisection (or BisectionExact), when the x values are of type Float64, Float32, or Float16, the default tolerances are zero and there is no limit on the number of iterations or function evalutions. In this case, the algorithm is guaranteed to converge to an exact zero, or a point where the function changes sign at one of the answer's adjacent floating point values.

For other types, the A42 method (with its tolerances) is used.

default_tolerances(::AbstractAlefeldPotraShi, T, S)

The default tolerances for Alefeld, Potra, and Shi methods are xatol=zero(T), xrtol=eps(T)/2, atol= zero(S), and rtol=zero(S), with appropriate units; maxevals=45, maxfnevals = Inf; and strict=true.

secant_method(f, xs; [atol=0.0, rtol=8eps(), maxevals=1000])

Perform secant method to solve f(x) = 0.

The secant method is an iterative method with update step given by b - fb/m where m is the slope of the secant line between (a,fa) and (b,fb).

The inital values can be specified as a pair of 2, as in (a,b) or [a,b], or as a single value, in which case a value of b is chosen.

The algorithm returns m when abs(fm) <= max(atol, abs(m) * rtol). If this doesn't occur before maxevals steps or the algorithm encounters an issue, a value of NaN is returned. If too many steps are taken, the current value is checked to see if there is a sign change for neighboring floating point values.

The Order1 method for find_zero also implements the secant method. This one will be faster, as there are fewer setup costs.

Examples:

Roots.secant_method(sin, (3,4))
Roots.secant_method(x -> x^5 -x - 1, 1.1)

Note:

This function will specialize on the function f, so that the inital call can take more time than a call to the Order1() method, though subsequent calls will be much faster. Using FunctionWrappers.jl can ensure that the initial call is also equally as fast as subsequent ones.

bisection(f, a, b; [xatol, xrtol])

Performs bisection method to find a zero of a continuous function.

It is assumed that (a,b) is a bracket, that is, the function has different signs at a and b. The interval (a,b) is converted to floating point and shrunk when a or b is infinite. The function f may be infinite for the typical case. If f is not continuous, the algorithm may find jumping points over the x axis, not just zeros.

If non-trivial tolerances are specified, the process will terminate when the bracket (a,b) satisfies isapprox(a, b, atol=xatol, rtol=xrtol). For zero tolerances, the default, for Float64, Float32, or Float16 values, the process will terminate at a value x with f(x)=0 or f(x)*f(prevfloat(x)) < 0 or f(x) * f(nextfloat(x)) < 0. For other number types, the A42 method is used.

muller(f, xᵢ; xatol=nothing, xrtol=nothing, maxevals=100)
muller(f, xᵢ₋₂, xᵢ₋₁, xᵢ; xatol=nothing, xrtol=nothing, maxevals=100)

Muller’s method generalizes the secant method, but uses quadratic interpolation among three points instead of linear interpolation between two. Solving for the zeros of the quadratic allows the method to find complex pairs of roots. Given three previous guesses for the root xᵢ₋₂, xᵢ₋₁, xᵢ, and the values of the polynomial f at those points, the next approximation xᵢ₊₁ is produced.

Excerpt and the algorithm taken from

W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery Numerical Recipes in C, Cambridge University Press (2002), p. 371

Convergence here is decided by xᵢ₊₁ ≈ xᵢ using the tolerances specified, which both default to eps(one(typeof(abs(xᵢ))))^4/5 in the appropriate units. Each iteration performs three evaluations of f. The first method picks two remaining points at random in relative proximity of xᵢ.

Note that the method may return complex result even for real intial values as this depends on the function.

Examples:

muller(x->x^3-1, 0.5, 0.5im, -0.5) # → -0.500 + 0.866…im
muller(x->x^2+2, 0.0, 0.5, 1.0) # → ≈ 0.00 - 1.41…im
muller(x->(x-5)*x*(x+5), rand(3)...) # → ≈ 0.00
muller(x->x^3-1, 1.5, 1.0, 2.0) # → 2.0, Not converged

Roots.newton(f, fp, x0; kwargs...)

Implementation of Newton's method: x_n1 = x_n - f(x_n)/ f'(x_n)

Arguments:

• f::Function – function to find zero of

• fp::Function – the derivative of f.

• x0::Number – initial guess. For Newton's method this may be complex.

With the FowardDiff package derivatives may be computed automatically. For example, defining D(f) = x -> ForwardDiff.derivative(f, float(x)) allows D(f) to be used for the first derivative.

Keyword arguments are passed to find_zero using the Roots.Newton() method.

See also Roots.newton((f,fp), x0) andRoots.newton(fΔf, x0)` for simpler implementations.

dfree(f, xs)

A more robust secant method implementation

Solve for f(x) = 0 using an alogorithm from Personal Calculator Has Key to Solve Any Equation f(x) = 0, the SOLVE button from the HP-34C.

This is also implemented as the Order0 method for find_zero.

The inital values can be specified as a pair of two values, as in (a,b) or [a,b], or as a single value, in which case a value of b is computed, possibly from fb. The basic idea is to follow the secant method to convergence unless:

• a bracket is found, in which case bisection is used;

• the secant method is not converging, in which case a few steps of a quadratic method are used to see if that improves matters.

Convergence occurs when f(m) == 0, there is a sign change between m and an adjacent floating point value, or f(m) <= 2^3*eps(m).

A value of NaN is returned if the algorithm takes too many steps before identifying a zero.

Examples

Roots.dfree(x -> x^5 - x - 1, 1.0)

fzero(f, x0; order=0; kwargs...)
fzero(f, x0, M; kwargs...)
fzero(f, x0, M, N; kwargs...)
fzero(f, x0; kwargs...)
fzero(f, a::Number, b::Numbers; kwargs...)
fzero(f, a::Number, b::Numbers; order=?, kwargs...)
fzero(f, fp, a::Number; kwargs...)

Find zero of a function using one of several iterative algorithms.

• f: a scalar function or callable object

• x0: an initial guess, a scalar value or tuple of two values

• order: An integer, symbol, or string indicating the algorithm to use for find_zero. The Order0 default may be specified directly by order=0, order=:0, or order="0"; Order1() by order=1, order=:1, order="1", or order=:secant; Order1B() by order="1B", etc.

• M: a specific method, as would be passed to find_zero, bypassing the use of the order keyword

• N: a specific bracketing method. When given, if a bracket is identified, method N will be used to finish instead of method M.

• a, b: When two values are passed along, if no order value is specified, Bisection will be used over the bracketing interval (a,b). If an order value is specified, the value of x0 will be set to (a,b) and the specified method will be used.

• fp: when fp is specified (assumed to compute the derivative of f), Newton's method will be used

• kwargs...: See find_zero for the specification of tolerances and other keyword arguments

Examples:

fzero(sin, 3)                  # use Order0() method, the default
fzero(sin, 3, order=:secant)   # use secant method (also just `order=1`)
fzero(sin, 3, Roots.Order1B()) # use secant method variant for multiple roots.
fzero(sin, 3, 4)               # use bisection method over (3,4)
fzero(sin, 3, 4, xatol=1e-6)   # use bisection method until |x_n - x_{n-1}| <= 1e-6
fzero(sin, 3, 3.1, order=1)    # use secant method with x_0=3.0, x_1 = 3.1
fzero(sin, (3, 3.1), order=2)  # use Steffensen's method with x_0=3.0, x_1 = 3.1
fzero(sin, cos, 3)             # use Newton's method

Note: unlike find_zero, fzero does not specialize on the type of the function argument. This has the advantage of making the first use of the function f faster, but subsequent uses slower.