# References and implementation for safe/hybrid univariate Newton's method
**URL:** https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638
**Category:** Numerics
**Tags:** question
**Created:** [September 22, 2022, 12:21pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638 "2022-09-22T12:21:23Z")
**Posts on this page:** 15
**Page:** 1
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [September 22, 2022, 12:21pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/1 "2022-09-22T12:21:23Z")
There exist variants of Newton’s method with safeguards against non-convergence (escape bracket, not converging “fast enough”, etc, in which case they switch to a bracketing method like Brent’s or bisection). Various implementations I found differ in details, and I could not find a paper with convergence analysis.
1. Is there a reference that discusses these methods in detail?
2. Is there a Julia implementation?
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [September 22, 2022, 12:36pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/2 "2022-09-22T12:36:20Z")
For a while, anytime I needed to solve a new univariate problem I tried different univariate solvers and always found that Brent’s method (from Optim.jl) was performing the best. So now I just always reach for that.
I do remember seeing recently a post here about someone implementing [the ITP method](https://en.wikipedia.org/wiki/ITP_method) in a new package but I have not given it a try. This method would be more similar to what you described.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [September 22, 2022, 12:52pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/3 "2022-09-22T12:52:12Z")
No, the ITP method is not really similar, as it does not use derivatives, and it is [already implemented](https://github.com/JuliaMath/Roots.jl/pull/274) in Roots.jl (thanks @TheLateKronos!)
Just to clarify: I am not asking for recommendations on rootfinding methods.
The question is about a very specific family of methods, and is mainly motivated by curiosity. Two general recommendations I encounter are
1. get close to a root, then use Newton-Raphson for “polishing”,
2. try Newton-Raphson, fall back to some bracketing method if it “fails”
But in practice I find it hard to implement these in a robust and practical way, because the criteria for switching are quite fuzzy. This is a shame, since for some functions derivatives are readily available, so a bracketing method that uses derivatives would be useful.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [September 22, 2022, 1:18pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/4 "2022-09-22T13:18:27Z")
The Numerical Recipes book, chapter 9, covers this topic in detail.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [September 22, 2022, 1:24pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/5 "2022-09-22T13:24:12Z")
If you are referring to `rtsafe` in Section 9.4 of the 3rd edition, then it doesn’t. It is _one_ of the _ad-hoc_ implementations floating around. All you get is the description
> The hybrid algorithm takes a bisection step whenever Newton-Raphson would take the solution out of bounds, or whenever Newton-Raphson is not reducing the size of the brackets rapidly enough.
but a convergence analysis or a justification for the criteria is missing.
Again, it is very easy (and tempting) to make ad hoc modifications to existing algorithms. Unfortunately, these can easily backfire.
---
### Author:  [@dlakelan](https://discourse.julialang.org/u/dlakelan)
#### Post date: [September 22, 2022, 1:32pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/6 "2022-09-22T13:32:20Z")
> [@Tamas\_Papp](#):
>
> convergence analysis
If we are talking about a continuous function and a bracket which is always at least cut in half, then we have a guarantee that the method converges by the intermediate value theorem. Taking the newton step when it makes the interval smaller than half it’s current size and the bisection step otherwise can only make convergence faster as the interval will always be smaller than or equal to half it’s initial size at every step.
Not sure what else you are looking for.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [September 22, 2022, 2:05pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/7 "2022-09-22T14:05:12Z")
> [@dlakelan](#):
>
> Taking the newton step when it makes the interval smaller than half it’s current size
I am not sure how that follows. In my copy of the book the step is taken if it does not go out of the bracket _and_
|2f(x\_m)| \< |x\_2 - x\_1||f'(x\_m)|
where (x\_1, x\_2) is the original bracket, and x\_m = (x\_1 + x\_2)/2.
The new value is x\_m - f(x\_m)/f'(x\_m), which could be close to x\_1 or x\_2.
To make it a bit more concrete, consider x\_1 = 0, x\_2 = 1, f(0) = -1, f(1) = 2, f(0.5) = 1, f'(0.5) = a. The second condition for accepting the Newton step is
2 \< |a|
while the resulting point is 0.5 - 2/a. Now let a = 4.1, so that x\_3 \approx 0.01. If f(x\_3) \< 0, (x\_3, x\_2) will be the new bracket, and is about 99% of the old one.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [September 22, 2022, 3:24pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/8 "2022-09-22T15:24:46Z")
Following the notation in the 3rd edition of the book, if I am not mistaken, the next point computed seems to be 0.378:
```julia
x1, x2 = 0.0, 1.0
fl, fh = -1.0, 2.0
xl, xh = x1, x2
rts = 0.5 * (x1 + x2) # = 0.5
dxold = abs(x2-x1)
dx = dxold
f = 0.5
df = 4.1
(((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0) || (abs(2.0*f) > abs(dxold*df)) # false
# so take the Newton step:
dxold = dx;
dx = f/df
temp = rts
rts -= dx # new rts = 0.378
```
---
### Author:  [@ctkelley](https://discourse.julialang.org/u/ctkelley)
#### Post date: [September 22, 2022, 5:09pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/9 "2022-09-22T17:09:54Z")
I vaguely remember W. Kahan giving a talk about how he got this idea into the solver (the **solve** button) for the HP35 programmable calculator. This was from the late 1970s and I do not know if he published it or not. I think the opening gambit was bisection and then the solver switched to Newton when the intervals got small.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [September 22, 2022, 6:48pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/10 "2022-09-22T18:48:09Z")
> [@ctkelley](#):
>
> W. Kahan giving a talk about how he got this idea into the solver (the **solve** button) for the HP35 programmable calculator. This was from the late 1970s and I do not know if he published it or not.
[This article by W. Kahan](https://www.keesvandersanden.nl/calculators/hp_journals/HP_Journal_7912_Personal_Calculator_has_Key_to_Solve_Any_Equation.pdf) about the HP-34C solver seems to be based on the secant method.
---
### Author:  [@ctkelley](https://discourse.julialang.org/u/ctkelley)
#### Post date: [September 22, 2022, 6:50pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/11 "2022-09-22T18:50:59Z")
That’s the one. Thanks for doing the error correction Newton → secant from my fading memory.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [September 22, 2022, 7:32pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/12 "2022-09-22T19:32:00Z")
That is implemented in the default Order0 method of find\_zero in Roots. There is also the possibility of other hybrid methods, switching to a bracketing method if one is identified.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [September 23, 2022, 2:13pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/13 "2022-09-23T14:13:14Z")
> [@rafael.guerra](#):
>
> if I am not mistaken, the next point computed seems to be 0.378
Good point, thanks for working through the numbers (which I should have done before posting the example 😊 ).
Try f(0.5) = -0.5 and f´(0.5) = -(1 + \epsilon), which gets you arbitrarily close to 0, so the new bracket is (\approx 0, 1). Eg
```julia
f = -0.5
df = -1.1
```
gives `rts = 0.04545...`
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [September 23, 2022, 3:18pm UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/14 "2022-09-23T15:18:16Z")
I forget if the convergence analysis is in detail, but in the [LithBoonkkampIJzermanBracket](https://github.com/JuliaMath/Roots.jl/blob/master/src/Derivative/lith.jl#L288) method of `find_zero` from `Roots` there is an implementation of this kind of algorithm and a reference to follow up on.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [September 26, 2022, 9:58am UTC](https://discourse.julialang.org/t/references-and-implementation-for-safe-hybrid-univariate-newtons-method/87638/15 "2022-09-26T09:58:34Z")
Thanks, the article ([non-paywall version here](https://arxiv.org/pdf/1702.03174.pdf)) was extremely informative. This is the kind of answer I was looking for, and I am delighted to learn that it is already implemented in the excellent Roots.jl.