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

Specific Domains Numerics
, ,
1

I am trying to increase precision for a numeric solution using Roots.jl using ArbFloat (of ArbNumerics.jl).
But, when I decrease the tolerance, the outcome does not change.
I keep comparing it to the numerical solution given by Mathematica for a precision of 100 bits.

using Roots, ArbNumerics

setprecision(ArbFloat, bits=1000)
# define the functions
f(x) = elliptic_e(x) - 1.3
fp(x) = (elliptic_e(x) - elliptic_k(x))/(2*x)
f_df_for_problem = (f,fp)
prob = ZeroProblem(f_df_for_problem, ArbFloat.((0., 1.)))

# the Mathematica solution:
mat_sol = ArbFloat(0.5970995261535437979056290157285245402529726426867049646646376914313944654874531279615647217356979859)

# checking the tolerances:
tol_arr = ArbFloat.([1e-5, 1e-10, 1e-15, 1e-20, 1e-25, 1e-30])
sol_diff = zeros(ArbFloat{}, 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) - mat_sol
end
diff(sol_diff)
# 5-element Vector{ArbFloat{1024}}:
# -1.101830987403286005066785079340465546840396709080396701806064417093954217433957889770394422941865029374686504172338753893173456687130137027492372054459748648498572128908751908720015249250176160412383446695027162276445486575283091025261380542647988174481266056856722571444701223523899719468795845530662e-6
#  0
# -4.890235796098201158338650651658624761209184554431080843729011388461680467110056575333933006390200023438226437943815995258787281268074457177754306036569641977736232991446959605210152128331054802662009701025783145520192069461248260365300342805569393266846881813550215442269927455543559205690749057359874e-18
#  0
#  0

How can I increase the precision of this example?

2

I suspect that the parser parses the number as a Float64 and passes it to ArbFloat constructor. You probably want that to be in a string or something.

2 Likes
3

If I change the code to ignore mat_sol

tol_arr = ArbFloat.([1e-5, 1e-10, 1e-15, 1e-20, 1e-25, 1e-30])
sol_diff = zeros(ArbFloat{}, 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

I get the same thing

diff(sol_diff )
# 5-element Vector{ArbFloat{1024}}:
#  -1.101830987403286005066785079340465546840396709080396701806064417093954217433957889770394422941865029374686504172338753893173456687130137027492372054459748648498572128908751908720015249250176160412383446695027162276445486575283091025261380542647988174481266056856722571444701223523899719468795845530662e-6
#   0
# -4.890235796098201158338650651658624761209184554431080843729011388461680467110056575333933006390200023438226437943815995258787281268074457177754306036569641977736232991446959605210152128331054802662009701025783145520192069461248260365300342805569393266846881813550215442269927455543559205690749057359874e-18
#   0
#   0
4

f(x) = elliptic_e(x) - 1.3 this is the Float64(1.3). You might want to use 13//10 (an exact rational that will be promoted to the appropriate type)

1 Like
5

dlakelan and Oscar_Smith are giving you good advice.

using ArbNumerics
setprecision(ArbFloat, bits=1000) # you do not need 1000 bits, I expect.

# the Mathematica solution:
mat_sol_str = "0.5970995261535437979056290157285245402529726426867049646646376914313944654874531279615647217356979859"
mat_sol = ArbFloat(mat_sol_str)

# to keep (ensure, force) your functions remain good with ArbNumeric processing
# again, perhaps unnecessary -- maybe helpful as you investigate
thirteen_tenths = ArbFloat(13.0) / ArbFloat(10.0)
f(x) = elliptic_e(ArbFloat(x)) - thirteen_tenths
fp(x) = fpx(ArbFloat(x))
fpx(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)

etc

1 Like
6

Thanks @JeffreySarnoff @Oscar_Smith and @dlakelan!

I tried to make sure that everything is ArbFloat{1024}, but it does not seem to work. I get the same results as before. If you have more ideas I’d love to hear!

setprecision(ArbFloat, bits=1000)
# both of these options do not work:
thirteen_tenths = ArbFloat(13.0) / ArbFloat(10.0)
thirteen_tenths = 13//10
# trying to set the functions propely
f(x) = elliptic_e(x) - thirteen_tenths
# fp(x) = (elliptic_e(x) - elliptic_k(x))/(2*x)
fp(x) = fpx(ArbFloat(x))
fpx(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)
f_df_for_problem = (f,fp)

prob = ZeroProblem(f_df_for_problem, ArbFloat.((0., 1.)))

tol_arr = ArbFloat.([1e-5, 1e-10, 1e-15, 1e-20, 1e-25, 1e-30])
sol_diff = zeros(ArbFloat{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
7

using Wolfram Alpha
with f(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)
f(0.0) is indeterminate
f(1.0) is complex infinity

1 Like
8

good point!
But I even when I change the boundaries the issue persists:

setprecision(ArbFloat, bits=1000)
# thirteen_tenths = ArbFloat(13.0) / ArbFloat(10.0)
f(x) = elliptic_e(x) - 13//10
# fp(x) = (elliptic_e(x) - elliptic_k(x))/(2*x)
fp(x) = fpx(ArbFloat(x))
fpx(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)
f_df_for_problem = (f,fp)

# taking proper boundaries
prob = ZeroProblem(f_df_for_problem, ArbFloat.((2//10, 7//10)))

tol_arr = ArbFloat.([1e-5, 1e-10, 1e-15, 1e-20, 1e-25, 1e-30])
sol_diff = zeros(ArbFloat{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

diff(sol_diff)
# 5-element Vector{ArbFloat{1024}}:
#   0
#  -1.263691682891297874032500101802607475645278384356269074460873730882216789672564541790403977384163012574771779868569275528853410011912960808494838703525069307787221051902872795671316331610025642341529374557913699031175136942578446693471510936037011017742278852572178006157009005170082418898873784590014816036e-13
#   0
#   0
#   0
9

do you have a rough (less high precision) answer to what you are solving?

[It is unlikely that the ArbNumerics functions are working improperly,
and there is the possibility that something in the way you have structured
the task has overlooked a computational detail which is behind
“The issue that persists”.]

I have no idea what you are trying to determine, nor why.
What we need is a very-much-stripped-down exemplar
(for this purpose there is no need to test various tolerances
nor any need to iterate over whathaveyou, the solvers
do that for you.) with some explaination.

1 Like
10

What I am giving here is just some minimal example. My goal is to invert a function with very high precision. The problem is that I did not manage yet to increase the precision. I was thinking that the way to go is to reduce the tolerances. But I see that reducing it does not affect the outcome. That is why I gave this simpler example.

I agree that I am probably doing something wrong.

The point is not to say that ArbNumerics has a bug, but that I do not know how to work with Roots.jl and ArbNumerics.jl to increase the precision of the inverted function.

11

Are you trying to invert f(x) = (elliptic_e(x) - elliptic_k(x)) / (2*x)?
If not, write the formula with elliptic functions that you want to invert … in math rather than in Julia

1 Like
12

I am trying to invert this function
f(x, \lambda)= \frac{4 \sqrt{2} }{3} \sqrt{-q_1^2} \left( \lambda^2 E_e \left( \frac{q_2^2}{q_1^2} \right) + \sqrt{4 x + \lambda^{4} } E_k \left( \frac{q_2^2}{q_1^2} \right) \right)
where
q_1^2 = \lambda^2 - \sqrt{4 x + \lambda^4} < 0
q_2^2 = \lambda^2 + \sqrt{4 x + \lambda^4} >0
\lambda > 0 and -\lambda^4/4 < x<0.
And E_k(m), E_e(m) are the complete elliptic integrals of the first and second kind.
I consider \lambda as a parameter, and I want to solve for x(f, \lambda = \lambda_0).

13

Check out this example

setprecision(ArbFloat, bits=1000)

function get_root(y::ArbReal)
    f(x) = elliptic_e(x) - y
    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,fp)
    f_df_for_problem = (f_arb,fp_arb)

    prob = ZeroProblem(f_df_for_problem, ArbReal.((1//10, 9//10)))
    solve(prob, Roots.LithBoonkkampIJzermanBracket())
end

and then

sol_of_root = get_root(ArbReal(13//10))

typeof(sol_of_root)
# ArbFloat{1000}

So the solve function returns an ArbFloat even though all variables are ArbReal. Could it be messing up with the precision?

14

no – ArbReals are are intervals, consider ArbFloats the midpoints of such intervals.
Those intervals are generally tight enough. That is not a problem.
If you prefer, just use ArbReal throughout. Then the result will include the true result (assuredly).

1 Like
15

Many root-finding schemes will use a Newton method when close to the answer. Since Newton has quadratic convergence, it will stop at the same iteration for a range of tolerance values which includes progressively more terms in a series like the one you are using. A simple Newton scheme converges nicely for your example problem, so Arb[Numerics] seems to be doing fine here. (However it deviates from your Mathematica result at 256 bits - is that the precision they use?) It is still possible that Roots.jl is not fully type-consistent, but you need to change your test to show that.

2 Likes
16

What do you mean? I thought I was using ArbReal throughout…

17

OH – I had not seen the revision.

18

Thank you for the math. Now for some baby steps for me:

Before inverting the function, what is it doing – to what end.

How does x go from a variable (an argument to f(x,) that constrains bounds) to a function ( to solve for x(f,) )?

1 Like
19

I believe @Ralph_Smith is correct in why you get confusing results (with the ArbFloat example). 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.

As a starting point let me compute the root to high precision (using my personal ArbExtras.jl package).

using SpecialFunctions, Arblib, ArbExtras

setprecision(Arb, 1024)

# Note that ellipe doesn't support Arb arguments so we go through the complex version
f(x::Arb) = real(ellipe(Acb(x))) - 13 // 10
df(x::Arb) = real(ellipe(Acb(x)) - ellipk(Acb(x)))/(2*x)

r0 = setball(Arb, 0.597, 1e-3) # Some initial enclosure of the root

r = ArbExtras.refine_root(f, r0, df = df, verbose = true)

This gives the output

[ Info: enclosure: [0.60 +/- 4.01e-3]
[ Info: enclosure: [0.59710 +/- 4.44e-6]
[ Info: enclosure: [0.5970995262 +/- 7.95e-11]
[ Info: enclosure: [0.5970995261535437979056 +/- 4.12e-23]
[ Info: enclosure: [0.5970995261535437979056290157285245402529726426867 +/- 2.12e-50]
[ Info: enclosure: [0.597099526153543797905629015728524540252972642686704964664637691431394465487456318730673086362723578116701 +/- 5.30e-106]
[ Info: enclosure: [0.59709952615354379790562901572852454025297264268670496466463769143139446548745631873067308636272357811670124731820824719418900915095793285799996819422478713629021168865105270107934991811685519209186640529643684042652259 +/- 9.09e-219]
[ Info: enclosure: [0.597099526153543797905629015728524540252972642686704964664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561743 +/- 8.78e-307]
[ Info: enclosure: [0.597099526153543797905629015728524540252972642686704964664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561743 +/- 8.50e-307]
[ Info: diameter only improved from (637534209 * 2^-1047) to (620756993 * 2^-1047) - stopping early
[0.597099526153543797905629015728524540252972642686704964664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561743 +/- 8.50e-307]

We can see that the Mathematica solution is only accurate to about 80 digits.

julia> r - Arb("0.5970995261535437979056290157285245402529726426867049646646376914313944654874531279615647217356979859")
[3.190769108364627025592216701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561743e-78 +/- 8.53e-307]

This is also the same solution you code finds

setprecision(ArbFloat, bits=1000)
f(x) = elliptic_e(x) - 13 // 10
df(x) = (elliptic_e(x) - elliptic_k(x))/(2*x)
f_df_for_problem = (f, df)

prob = ZeroProblem(f_df_for_problem, ArbFloat.((0., 1.)))

tol_arr = ArbFloat.([1e-5, 1e-10, 1e-15, 1e-20, 1e-25, 1e-30])
sol_arr = map(tol_arr) do tol
    solve(prob, Roots.LithBoonkkampIJzermanBracket(), rtol = tol)
end

# Very high precision enclosure of root
r = ArbFloat("0.597099526153543797905629015728524540252972642686704964664637691431394465487456318730673086362723578116701247318208247194189009150957932857999968194224787136290211688651052701079349918116855192091866405296436840426522586446300241603087130298569879618160424746435212787624616426552989938655802225368867561743")

# Show error for roots
sol_arr .- r

This gives the output

6-element Vector{ArbFloat{1024}}:
 1.101830987408181890667755211574514881826117882386650477156856095255801747838866539359810422200848823088037092438226759563013846778015117614144511412083260196085434790806715349605824142279781614219020667403865081804905862617288346518126978284084381903510854053651780045544266684946154934721956960778503e-6
 4.890235796098269886094746750943470118688768023424265459778149270207280703880883637756928902897501946940702325383547627861579105626711762075959504672976726128831360057337088609386672293334896682530707328643737555808770885624950154112424482385206495461376641200642920707794266798310707653185660986592683e-18
 4.890235796098269886094746750943470118688768023424265459778149270207280703880883637756928902897501946940702325383547627861579105626711762075959504672976726128831360057337088609386672293334896682530707328643737555808770885624950154112424482385206495461376641200642920707794266798310707653185660986592683e-18
 1.245373578596568724871677098026871781915477713895901244737175565618205629160447838989352550644230462223862825279454214426092289169636327266777080919747716291758468599027274227774655613574232951647168308530825712622544584789621990300686490037973155307065048439975236284273010927759324897508834604413348e-52
 1.245373578596568724871677098026871781915477713895901244737175565618205629160447838989352550644230462223862825279454214426092289169636327266777080919747716291758468599027274227774655613574232951647168308530825712622544584789621990300686490037973155307065048439975236284273010927759324897508834604413348e-52
 1.245373578596568724871677098026871781915477713895901244737175565618205629160447838989352550644230462223862825279454214426092289169636327266777080919747716291758468599027274227774655613574232951647168308530825712622544584789621990300686490037973155307065048439975236284273010927759324897508834604413348e-52

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.

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.

2 Likes
20

Thanks for the valuable comment!
Indeed the Mathematica result was not accurate to 100 digits (I copied the wrong thing). This is the Mathematica result up to 100 digits:

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

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