# How to increase precision of solution using Roots.jl and ArbNumerics.jl

**URL:** <https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850>\
**Category:** Numerics\
**Tags:** question, roots, precision\
**Created:** [August 16, 2023, 2:55am UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850 "2023-08-16T02:55:24Z")\
**Posts on this page:** 9\
**Page:** 2

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [August 17, 2023, 2:52pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/21 "2023-08-17T14:52:45Z")

</div>

> [@JeffreySarnoff](#):
>
> Before inverting the function, what is it doing – to what end.

I am inverting it, because I need to solve some differential equations.

> [@JeffreySarnoff](#):
>
> How does x go from a variable (an argument to f(x,_) that constrains bounds) to a function ( to solve for x(f,_) )?

I have some value of f, say f\_0, and I need to find the x that solves f(x) = f\_0. So basically, I am trying to implement the root-finding function to give me this inverse. I am not sure this answers your question…

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [August 17, 2023, 3:12pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/22 "2023-08-17T15:12:54Z")

</div>

> [@joeldahne](#):
>
> The Newton iterations are converging quadratically but you are only increasing the tolerance linearly, so most of the time it does the same number of iterations. Indeed the Mathematica result is also not accurate to 100 digits.

Thanks! you are indeed correct as I wrote

> [@roi.holtzman](#):
>
> Indeed the Mathematica result was not accurate to 100 digits (I copied the wrong thing). This is the Mathematica result up to 100 digits:
> 
> ```julia
> mat_sol = 5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781
> 
> ```
> 
> I still do not understand why I am getting that the deviation from the mathematical result is so big compared to the tolerance I am using
> 
> ```julia
> setprecision(ArbReal, bits=2000)
> 
> mat_sol = ArbReal(0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781)
> 
> f(x) = elliptic_e(x) - 13//10
> fp(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)
> 
> f_arb(x) = f(ArbReal(x))
> fp_arb(x) = fp(ArbReal(x))
> f_df_for_problem = (f_arb,fp_arb)
> 
> prob = ZeroProblem(f_df_for_problem, ArbFloat.((2//10, 7//10)))
> 
> tol_arr = ArbReal.([1e-10, 1e-20, 1e-40, 1e-80, 1e-160, 1e-320])
> sol_diff = zeros(ArbReal{1024}, length(tol_arr))
> for (i, tol) in enumerate(tol_arr)
> sol_diff[i] = solve(prob, Roots.LithBoonkkampIJzermanBracket(), atol=tol, rtol=tol, xatol=tol, xrtol=tol)
> end
> 
> sol_diff .- mat_sol
> # 6-element Vector{ArbReal{1024}}:
> # 1.2638671375145267156858622267668263146486194631103066439239457343927831910566533896936293116052288778029257199931559187656131795963980324568386040031356265232081227664948940326357833243023160430254889004060107464311764950603434757473331346283872256753547221429720703375142021224453030934691368662734391262496e-13
> # 1.7545462322884165336212496421883900334107875403756946307200351056640138408884790322533422106586522815394012458664323675976958448507164834376529961055721542090171459202123696446699269229040068395952584809704740000135811776502905386162369235021465761244329039989233135719311727522067456789461205007606883446998e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426553039569708383033944653456941061581000481656896e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561742567542124464703805e-17
> 
> ```
> 
> ```julia
> 
> ```
> 
> ```julia
> mat_sol = 5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781
> 
> ```

> [@joeldahne](#):
>
> Notice how the error jumps from `1e-18` to `1e-52` between one step. Then it doesn’t improve further due to already satisfying the required tolerance.

I fixed my example with the actual Mathematica proper result, but I do not get the same results as you. I do not understand why the deviation in my case is so large compared to yours. Do you?

> [@roi.holtzman](#):
>
> ```julia
> setprecision(ArbReal, bits=2000)
> 
> mat_sol = ArbReal(0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781)
> 
> f(x) = elliptic_e(x) - 13//10
> fp(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)
> 
> f_arb(x) = f(ArbReal(x))
> fp_arb(x) = fp(ArbReal(x))
> f_df_for_problem = (f_arb,fp_arb)
> 
> prob = ZeroProblem(f_df_for_problem, ArbFloat.((2//10, 7//10)))
> 
> tol_arr = ArbReal.([1e-10, 1e-20, 1e-40, 1e-80, 1e-160, 1e-320])
> sol_diff = zeros(ArbReal{1024}, length(tol_arr))
> for (i, tol) in enumerate(tol_arr)
> sol_diff[i] = solve(prob, Roots.LithBoonkkampIJzermanBracket(), atol=tol, rtol=tol, xatol=tol, xrtol=tol)
> end
> sol_diff .- mat_sol
> # 6-element Vector{ArbReal{1024}}:
> # 1.2638671375145267156858622267668263146486194631103066439239457343927831910566533896936293116052288778029257199931559187656131795963980324568386040031356265232081227664948940326357833243023160430254889004060107464311764950603434757473331346283872256753547221429720703375142021224453030934691368662734391262496e-13
> # 1.7545462322884165336212496421883900334107875403756946307200351056640138408884790322533422106586522815394012458664323675976958448507164834376529961055721542090171459202123696446699269229040068395952584809704740000135811776502905386162369235021465761244329039989233135719311727522067456789461205007606883446998e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578117626680882423863314123512122397648188809999514616284468318022816285430076900650010820250236575985366472186173954302239428583719428163347192246596616159009013210278691731167307196736105564879687806472283979634912581903772354e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426553039569708383033944653456941061581000481656896e-17
> # 1.7545462322884165336079884568956236214664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561742567542124464703805e-17
> 
> ```

I thank you also for this last comment

> [@joeldahne](#):
>
> Regarding why the version with `ArbReal` returns `ArbFloat` I would guess it is because they use the type `float(ArbReal)` for the computations, which is `ArbFloat`. For most Julia types `float(T)` gives you a **more** capable type numerically so it is very commonly used is numerical code. For `ArbReal` it instead gives you a **less** capable type, which is a bit unfortunate.

I do not exactly understand what is more capable vs less capable. I guess it would be nice if the root-finding function would return an `ArbReal` with the error estimation from the root-finding algorithm.

In my other functions, I try to use `ArbReal`. Do you think since this is “less” capable, it is a bad choice?  
If I were to keep everything `ArbReal`, would you just convert the value that is returned by `solve` to `ArbReal`?

---

<div class="post-metadata">

**Author:** ![joeldahne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joeldahne/32/45861_2.png) [@joeldahne](https://discourse.julialang.org/u/joeldahne)\
**Post date:** [August 17, 2023, 5:30pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/23 "2023-08-17T17:30:58Z")

</div>

> I fixed my example with the actual Mathematica proper result, but I do not get the same results as you. I do not understand why the deviation in my case is so large compared to yours. Do you?

You are missing the citation marks around the value you give to `ArbReal`. If you replace it with

```julia
mat_sol = ArbReal("0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781")

```

it works fine!

> I do not exactly understand what is more capable vs less capable. I guess it would be nice if the root-finding function would return an `ArbReal` with the error estimation from the root-finding algorithm.

With more capable I mean that `ArbReal` is more capable than `ArbFloat`, it can do anything `ArbFloat` can do and at the same time it keeps track of the error. To say a bit more about this: The promotion system in Julia typically returns the more “general” of the input types, for example `promote_type(Int32, Int64) = Int64` and `promote_type(Float64, BigFloat) = BigFloat`. In base Julia I believe it always holds that `promote_type(T, float(T)) = float(T)`. But for `ArbReal` this is not the case, then `float(ArbReal) = ArbFloat` but `promote_type(ArbReal, ArbFloat) = ArbReal`.

In your case I believe you can just as well use `ArbFloat`. The only benefit with `ArbReal` is that you would notice if you lose too much precision in the numerical calculations, but it wont tell you anything about the error bound for the root anyway. If you want to get error bounds for the root you have to use some other algorithm. For example the interval Newton method, which is what the `refine_root` function I used in my previous post makes use of.

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [August 17, 2023, 6:15pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/24 "2023-08-17T18:15:21Z")

</div>

> [@joeldahne](#):
>
> You are missing the citation marks around the value you give to `ArbReal`.

This is the same error I mentioned early on. To understand what’s going on here, the Julia Parser reads the number as a Float64 thereby ignoring all the excess digits, then passes that to ArbReal … If you put the quote marks, the parser reads it as a string and passes the string to ArbReal which converts it to a full precision value. You simply can’t give a naked string of digits to the parser and get anything other than a native float type. The parser doesn’t know anything about Arb types.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [August 17, 2023, 7:16pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/25 "2023-08-17T19:16:39Z")

</div>

Thanks so much for the answer! It is very clear.

> [@joeldahne](#):
>
> If you want to get error bounds for the root you have to use some other algorithm. For example the interval Newton method, which is what the `refine_root` function I used in my previous post makes use of.

It would be great to get the error. This `refine_root` method is part of your package `ArbExtras` right? Do you know if there is a similar analog in `Roots`?

And thank you both @joeldahne and @dlakelan for pointing out the issue with the string parsing!

---

<div class="post-metadata">

**Author:** ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)\
**Post date:** [August 17, 2023, 7:40pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/26 "2023-08-17T19:40:29Z")

</div>

Thanks for this informative response. `Roots` uses `float` for the input values as it seems like the right generic for different number types.

When I try to solve this with `A42` which might be recommended as there is a known bracket, the convergence is quite quick for smooth enough functions, and the tolerance is like `2abs(r)*eps(one(r))` between iterates, I get an answer that is off from `mat_sol`, but yields an exact zero for `f`:

```julia
julia> setprecision(ArbFloat, bits=1000)
1000

julia> f(x) = elliptic_e(x) - 13 // 10
f (generic function with 1 method)

julia> r = find_zero(f, ArbReal.( (0.0, 1.0) ), A42())
0.5970995261535437979056290157285245402529726426867049646646376914313944654874563187306730863627235781167012473182082471941890091509579328579999681942247871362902116886510527010793499181168551920918664052964368404265225864463002416030871302985698796181604247464352127876246164265529899386558022253688676

julia> f(r)
0

julia> r - mat_sol
-2.496873060882422448783512389107976765165905935558302891285924378533353483003135167996405765830182705071755548296901317683611453726541415832228563754139160817061255022802266605200253564787212375383573447010061344197774631132438257432457875535296194985152788292500142950589819644532161239691396597893509e-78

```

I didn’t check with other zero-finding methods, and am not sure this happens with ArbReals which should be guaranteed intervals (as I understand it), but over floating point numbers some functions can have 0, 1 or more “zeros” near a mathematical zero.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [August 17, 2023, 10:05pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/27 "2023-08-17T22:05:47Z")

</div>

> [@j\_verzani](#):
>
> `Roots` uses `float` for the input values as it seems like the right generic for different number types.

Yes, that seems reasonable. Is there a way to get back from `Roots` the error estimate of the found solution?

---

<div class="post-metadata">

**Author:** ![joeldahne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joeldahne/32/45861_2.png) [@joeldahne](https://discourse.julialang.org/u/joeldahne)\
**Post date:** [August 18, 2023, 7:29am UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/28 "2023-08-18T07:29:42Z")

</div>

> [@roi.holtzman](#):
>
> It would be great to get the error. This `refine_root` method is part of your package `ArbExtras` right? Do you know if there is a similar analog in `Roots`?

The `refine_root` is indeed part of my `ArbExtras.jl` package, which uses `Arblib.jl` instead of `ArbNumerics.jl` for the interface to Arb. The code is fairly short and self-contained, so you could take a look at it.

In general you need to use different algorithms to get rigorous error bounds and I think that would be outside the scope of `Roots.jl`. I’m only aware of one package in the general repository that does this, [IntervalRootFinding.jl](https://github.com/JuliaIntervals/IntervalRootFinding.jl). But that package makes use of [`IntervalArithmetic.jl`](https://github.com/JuliaIntervals/IntervalArithmetic.jl) for the arithmetic and that doesn’t support the elliptic functions that you use. So if you want rigorous error bounds you will either have to use `ArbExtras.jl`, which is maybe not optimal since it is not in the general repository, or implement your own version of the interval Newton method, which is in principle a very simple method but might be hard if you are not used to working with rigorous numerics.

---

<div class="post-metadata">

**Author:** ![roi.holtzman](https://avatars.discourse-cdn.com/v4/letter/r/f05b48/32.png) [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)\
**Post date:** [August 18, 2023, 8:59pm UTC](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850/29 "2023-08-18T20:59:06Z")

</div>

Thanks so much for the information! I will have a look at `refine_root`.

[Previous page](https://discourse.julialang.org/t/how-to-increase-precision-of-solution-using-roots-jl-and-arbnumerics-jl/102850.md?page=1)