# Finding root numerically for a function defined in a finite range
**URL:** https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729
**Category:** Numerics
**Created:** [August 11, 2023, 6:12pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729 "2023-08-11T18:12:39Z")
**Posts on this page:** 17
**Page:** 1
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 11, 2023, 6:12pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/1 "2023-08-11T18:12:39Z")
I have a function f(x) defined in range (-a, \infty). f(x) is actually a piecewise function f\_1(x) for (-a,0), and f\_2(x) for (0, \infty).
For these two functions f\_1, f\_2 the 2nd derivative diverges at the boundary x=0.
I want to invert to find the function f^{-1}.
For any value I am interested in, I am trying to solve for the root using [Roots.jl](https://juliamath.github.io/Roots.jl/stable/). I tried to use Newton’s method and other bracketing methods.
However, I get into trouble. Either I get a solution that is not accurate to my precision level (and I cannot get better precision by changing the tolerance). Or sometimes, I get some solutions that do not make sense.
I would love to get some ideas on how to deal with this.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [August 11, 2023, 6:21pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/2 "2023-08-11T18:21:34Z")
Can you open an issue on `Roots.jl` and hopefully share the functions?
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 11, 2023, 6:25pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/3 "2023-08-11T18:25:58Z")
Sure. My functions are a bit complicated, but I hope I make it clear enough.
Also, let me mention that these include elliptic integrals, so I have been using `ArbNumerics` to compute them. I assume that this is fine for `Roots.jl`, right?
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [August 11, 2023, 6:49pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/4 "2023-08-11T18:49:20Z")
It might be an issue checking for convergence, but let’s see.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 11, 2023, 8:55pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/5 "2023-08-11T20:55:52Z")
In the region where my function is not defined, should I return `NaN`?
I think this also messes up things.
(I am working on an github issue in the meantime)
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [August 11, 2023, 9:57pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/6 "2023-08-11T21:57:34Z")
BIsection will work if you return Inf or -Inf (assuming some sign is known), but otherwise the algorithms tend to stop when an NaN is encountered as they propagate to the next iterate.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 11, 2023, 10:01pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/7 "2023-08-11T22:01:01Z")
I see, and what about Newton? Should I also use `Inf`?
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [August 11, 2023, 10:10pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/8 "2023-08-11T22:10:19Z")
You are likely to have the same issue, though perhaps the algorithm never sees these points if you have reasonable initial estimates.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 11, 2023, 10:19pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/9 "2023-08-11T22:19:13Z")
I see. I am not sure what to do then. Is there a way in Roots.jl to keep Newton’s method above or below some threshold?
When I try the bisection method it also does not work so I am not sure how to solve it.
BTW, I just opened an [issue](https://github.com/JuliaMath/Roots.jl/issues/393).
---
### Author:  [@longemen3000](https://discourse.julialang.org/u/longemen3000)
#### Post date: [August 12, 2023, 1:19am UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/10 "2023-08-12T01:19:39Z")
> I see. I am not sure what to do then. Is there a way in Roots.jl to keep Newton’s method above or below some threshold?
Actually, yes, but is not called newton anymore, try this instead:
```julia
f(x, p) = x^5 - x - p #function
fp(x, p) = 5x^4 - 1 #derivative
function fdf1(x, p=1) #function and derivative, if you can calculate those two at the same time
fx = f(x, p)
∂fx = fp(x)
return fx, fx / ∂fx
end
fdf2 = (f,fp) #another way to specify function and derivative
#we specify the bracketed problems
prob1 = ZeroProblem(fdf1, (1, 2))
prob2 = ZeroProblem(fdf2, (1, 2))
#and then we solve those
solve(prob1, Roots.LithBoonkkampIJzermanBracket(),1.0) # p = 1.0
```
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 13, 2023, 10:25pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/11 "2023-08-13T22:25:29Z")
Thanks! This is very useful!
My goal is now to increase the precision of the root that this finds. Do you have any suggestions on how to do it?
I am trying to use `ArbFloat` but there seems to be a [bug](https://github.com/JuliaMath/Roots.jl/issues/393) with it so I cannot yet test if this would improve the precision.
Do you have more ideas on how to improve the precision? Maybe adding the 2nd derivative?
---
### Author:  [@longemen3000](https://discourse.julialang.org/u/longemen3000)
#### Post date: [August 14, 2023, 5:50am UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/12 "2023-08-14T05:50:22Z")
Hmm, you could start with a bracketed step on reduced precision, and then switch to a [LithBoonkkampIJzerman](https://juliamath.github.io/Roots.jl/stable/reference/#Roots.LithBoonkkampIJzerman) method with high order in `S`.They even have this as a note:
> For the larger values of `S`, the expressions to compute the next value get quite involved. The higher convergence rate is likely only to be of help for finding solutions to high precision.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [August 14, 2023, 11:59am UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/13 "2023-08-14T11:59:10Z")
There was a missing promotion with the `LithBoonkkampIJzerman` method that just got addressed. But, the basic tradeoff is adding the derivative (or more points of memory) will speed convergence, but precision is likely going to be greater with a bracketing method. So it might be better to switch this suggestion and start with a higher-order method and finish with a bracketing method, should one be identified along the way.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 14, 2023, 2:44pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/14 "2023-08-14T14:44:06Z")
Thanks @j_verzani, it seems like there is no error now. However, I do not see any effect of decreasing the tolerance on the result.
This is my function
```julia
function E_of_I_dw_above(I::ArbFloat, λ::ArbFloat; Emax=ArbFloat(5), maxiters=20, atol=eps(), rtol=eps(), xatol=eps(), xrtol=eps())
Emin = 0.0
I, λ, Emin, Emax = promote(I, λ, Emin, Emax)
I, atol, rtol, xatol, xrtol = promote(I, atol, rtol, xatol, xrtol)
# making sure that Emax is large enough
while I_dw_above(Emax, λ) < I
# Emax is too small
Emax = 2*Emax
end
EdE_pair_for_prob = (E -> I_dw_above(E, λ) - I, E -> dIdE_dw_above(E, λ))
E_prob = ZeroProblem(EdE_pair_for_prob, (Emin, Emax))
solve(E_prob, Roots.LithBoonkkampIJzermanBracket(), maxiters=maxiters, atol=atol, rtol=rtol, xatol=xatol, xrtol=xrtol)
end
E_of_I_dw_above(I, λ; Emax=10.0, maxiters=20, atol=eps(), rtol=eps(), xatol=eps(), xrtol=eps()) = E_of_I_dw_above(ArbFloat(I), ArbFloat(λ); Emax=Emax, maxiters=maxiters, atol=atol, rtol=rtol, xatol=xatol, xrtol=xrtol)
```
And I get the following:
```julia
E_of_I_dw_above(0.5, 0.1)
# 0.02888908568904820499569305844526
E_of_I_dw_above(0.5, 0.1, atol=ArbFloat(1e-60), rtol=ArbFloat(1e-60), xatol=ArbFloat(1e-60), xrtol=ArbFloat(1e-60))
# 0.0288890856890482049956930584451
bits = 400
E_of_I_dw_above(ArbFloat(0.5, bits=bits), ArbFloat(0.1, bits=bits))
# 0.02888908568904820499569305844469
E_of_I_dw_above(ArbFloat(0.5, bits=bits), ArbFloat(0.1, bits=bits), atol=ArbFloat(1e-60), rtol=ArbFloat(1e-60), xatol=ArbFloat(1e-60), xrtol=ArbFloat(1e-60))
# 0.0288890856890482049956930584451
```
What am I doing wrong?
About your advise regarding the accuracy. I am now using a bracketing method, `LithBoonkkampIJzermanBracket`. Does this mean that I cannot improve the accuracy?
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 14, 2023, 2:49pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/15 "2023-08-14T14:49:05Z")
Thanks a lot!
I didn’t understand from the docs what is `S` and what is `D`. Could you explain what these are?
I thought that a bracketing method would be better since I am interested in finding values near the boundary of where my function is defined. Near the boundary, I have that the first derivative vanishes, and the 2nd derivative diverges. And I guess that would be hard for the algorithm.
---
### Author:  [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [August 14, 2023, 5:24pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/16 "2023-08-14T17:24:58Z")
For that method, the authors use S for the number of past steps remembered (S=1 for Newton, S=2 for secant) and D for the number of derivatives used (S=1 for Newton, S=0 for secant). Basically, increasing either increases the convergence rate at the expense of either more function calls or a more complicated update step or both.
As for what offers more precision, that is about the convergence criteria. For that, by default the bracketing methods – when available – generally have tighter criteria, as they only depend on the past x values getting close to each other and that is chosen tightly, as the convergence is guaranteed up to floating point accuracy. (Which for Arb floats might be very accurate, if my understanding is correct).
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [August 14, 2023, 5:28pm UTC](https://discourse.julialang.org/t/finding-root-numerically-for-a-function-defined-in-a-finite-range/102729/17 "2023-08-14T17:28:23Z")
Thanks again @j_verzani!
That explanation is very clear.
I guess I will stick with the bracketing method `LithBoonkkampIJzermanBracket` since I care most about accuracy. I will try to figure out why the accuracy is not improving when I set smaller tolerances using `ArbFloat` for the solver.