Finding root numerically for a function defined in a finite range

Specific Domains Numerics
1

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

2

Can you open an issue on Roots.jl and hopefully share the functions?

1 Like
3

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?

4

It might be an issue checking for convergence, but let’s see.

1 Like
5

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)

6

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.

7

I see, and what about Newton? Should I also use Inf?

8

You are likely to have the same issue, though perhaps the algorithm never sees these points if you have reasonable initial estimates.

9

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.

10

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:

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
1 Like
11

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

12

Hmm, you could start with a bracketed step on reduced precision, and then switch to a 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.

1 Like
13

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.

3 Likes
14

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

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:

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?

15

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.

16

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

2 Likes
17

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.