Query status of the solution when using Roots

Specific Domains Numerics
, , ,
1

Edited to make it more specific
I would like to query status of the solution when using the Roots package, before the function returns an answer. I can see in the source code of “find_zero.jl” that it has some functionality state.convergence_failed = true and also a function assess_convergence(method::Any, state::UnivariateZeroState{T,S}, options) where {T,S} , but there is no documentation on how to use it. I have given a MWE and my attempts at resolving it below.

My questions is:
How to query the status of the solution before it returns a value? e.g. I would like to use the Schroder method first, and if it fails to converge then it should switch to Halley method. If the roots are complex or the algorithm fails to converge then it should return NaN.

using Roots
function root_func(c;guess=0.01,max_evals=100)
    t = collect(0:length(c)-1)
    f(x) = sum(@views c[i]/(1+x)^t[i] for i in 1:length(c))
    #find derivatives of f(x)
    df(x) = sum(@views -t[i]*c[i]/(1+x)^(t[i]+1) for i in 1:length(c))
    ddf(x) = sum(@views (t[i]^2+t[i])*c[i]/(1+x)^(t[i]+2) for i in 1:length(c))
    return find_zero((f,df,ddf),guess,Roots.Schroder();maxevals=max_evals)
end

Example usage:
root_func([-100,10,100,1,1]

My attempts at checking status of the solution before switching methods or returning a value are given here:

using Roots
function root_func(c;guess=0.01,max_evals=100)
    t = collect(0:length(c)-1)
    f(x) = sum(@views c[i]/(1+x)^t[i] for i in 1:length(c))
    #find derivatives of f(x)
    df(x) = sum(@views -t[i]*c[i]/(1+x)^(t[i]+1) for i in 1:length(c))
    ddf(x) = sum(@views (t[i]^2+t[i])*c[i]/(1+x)^(t[i]+2) for i in 1:length(c))
    r = find_zero((f,df,ddf),guess,Roots.Schroder();maxevals=max_evals)
    if r.convergence_failed = false
       return r
    end
    r = find_zero((f,df,ddf),guess,Roots.Halley();maxevals=max_evals)
    if r.convergence_failed = false
       return r
     else return NaN
   end
end
2

The part of the code to modify would be the method here https://github.com/JuliaMath/Roots.jl/blob/master/src/find_zero.jl#L796 which implements an algorithm that checks if a bracket has occurred in which case the method is switched. For your interest, you would call update_state, check the conditions you wanted and if you want to switch, you would switch the method, update the state and continue with update_state until convergence or failure to converge. One caveat, I forget how much support there is for complex zeros, if I recall correctly that is limited to a few methods.

3

Thanks. Do I necessarily have to modify the source code of “find_zero.jl” for this? This is difficult for me as I am new to Julia.

All I want from my function is to switch from Schroder to Halley when it fails to converge rather than returning an error. And if Halley’s method also fails to converge than it should just return NaN rather than an error message.

4

This could be done along these lines (using a internal form of find_zero which returns NaN instead of an error):

using Roots
using ForwardDiff
D(f) = x -> ForwardDiff.derivative(f, float(x))

f(x) = (x^2+1)*(x-1)*(x-2)
a = 1.5

F = Roots.callable_function((f, D(f), D(D(f))))
M = Roots.Schroder()
state = Roots.init_state(M, F, a)
options = Roots.init_options(M, state)

val = find_zero(M, F, options, state)

if isnan(val)
    M = Roots.Halley()
    state = Roots.init_state(M, F, a)
    val = find_zero(M, F, options, state)
end

val

In this case, I’m not sure what functions will not satisfy Schroder, but will satisfy Halley, but presumably you know of one or two.

1 Like
5

Thank you! This solves my problem.

Schroder’s method fails for root_func([10,10] but Halley’s method would work using the code in the original post, though I am not sure if this a good example