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

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…

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?

I thank you also for this last comment

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?

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

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.

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.

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!

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> 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.

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. But that package makes use of 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.