Re: Exit codes for fzero

[Jaroslav Hajek]

On Wed, Feb 10, 2010 at 7:22 PM, Rik <address@hidden> wrote:
> Jaroslav,
>
> The fzero routine you contributed is very good, but I have a question about
> the exit code for certain pathological functions.
>
> I am looking at the function 1/(x - pi) for the range [0, 5].  The zero,
> such as it is, is at pi.  On either side of this the function is basically
> -/+ infinity.
>
> Running the following code I get an exit code of 1 which seems to indicate
> that the algorithm has converged successfully.
>
> [x,fval,info,output] = fzero(@(x) 1./(x-pi), 3)
> x =  3.1416
> fval =  7.5060e+14
> info =  1
> output =
> {
>  iterations =  85
>  funcCount =  87
>  bracket =
>
>     3.1416   3.1416
>
>  bracketf =
>
>    -1.1259e+15   7.5060e+14
>
> }
>
> When I run with fsolve I get a return code of zero indicating that there
> was a problem during the computation of the solution.
>
>  [X, FVEC, INFO, OUTPUT, FJAC] = fsolve (@(x) 1./(x-pi), 3)
> X = -2746.1
> FVEC = -3.6373e-04
> INFO = 0
> OUTPUT =
> {
>  iterations =  401
>  successful =  344
>  funcCount =  519
> }
>
> FJAC = -1.3230e-07
>
> Earlier versions of Octave were even more explicit.  From the 3.0.X branch
> I get "error: fsolve: number of function calls exceeded limit" when running
> fsolve.
>
> My question is whether fzero could be modified to indicate correctly when
> there was a problem with the algorithm.  I suspect in my case that fzero
> switched to a bisection algorithm and kept adding one bit of precision with
> each function evaluation while simultaneously flipping above and below the
> value of pi.  Eventually it collected the 53 bits of precision for IEEE
> floating point numbers and then returned.
>
> One thought I had would be a single sanity check at the end of fzero to
> verify that fval is somewhere near zero.  That would easily catch this case
> where the end value of 7x10^14 cannot be confused with zero.  Is there
> another better way to get an error exit code when the algorithm has failed?
>
> --Rik
>
>

Hi Rik,
your explanation is correct - the algorithm falls back to bisection
when inverse interpolation fails. Basically, the function is assumed
to be continuous, so based on that assumption the algorithm
incorrectly deduces that a root is bracketed (as if the -Inf and Inf
were connected).
This is where the actual bracketing interval may be useful. After
fsolve "converges", you can check out the final bracketing, and if the
jump in values is excessive, consider it to be a discontinuity. I
don't think it's much useful to put such a check directly into fzero;
the best way to react is probably problem-dependent.

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz