[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Exit codes for fzero
From: |
Rik |
Subject: |
Re: Exit codes for fzero |
Date: |
Tue, 16 Feb 2010 17:57:53 -0800 |
Jaroslav Hajek wrote:
>>
>> For a patch, I considered a simple level test, such as fval <
>> SOME_TOLERANCE or (fb - fa) < SOME_TOLERANCE. These would catch one type
>> of singularity where the right- or left-hand limit was undefined and going
>> to infinity. However, it wouldn't catch a simple jump discontinuity such
>> as a step function from -1 to +1. In that case both limits would exist and
>> both would be well-bounded but the two limits do not converge on the same
>> value. Instead, my patch looks at the numerically calculated slope
>> between the last two computed points of fzero. If the slope is extremely
>> high this *may* indicate a discontinuity of either the infinite or bounded
>> sort.
>>
>> A draft of the patch is attached. If it looks good and I get no other
>> comments I will commit it in a few days.
>>
>
> I have some comments. This looks good, but I think you should more
> carefully consider what an "extremely high" slope is.
> In particular, I don't like the absolute constant 1e6, because it
> makes the test scaling-dependent. Yes, absolute scaling independence
> is impossible due to limited FP range, but 1e6 is way too small. I
> think it would be better if the slope of the final bracketing were
> compared to the slope of the initial bracketing.
There is some scale dependence in the patch. I used a comparison to
max(1e6, 0.5/tolX) for the determination of "extremely steep". This breaks
the comparison into a fixed zone of 1e6 and a scale dependent zone.
My goal was to catch a jump discontinuity of 1 unit or a delta_y of 1. The
delta_x that fzero arrives at is basically the tolerance in x that the user
has specified. The slope for any jump discontinuity of this sort is
therefore 1/delta_x which is why I set one of the bounds as 0.5/tolX. This
scales with the user's requirements. For the default tolerance of 1e-8
this is 5E7. I don't know if this is as large as you want, but if it is
increased it won't catch the jump discontinuity. For comparison, a
function diverging to +/- infinity such as tan(x) or 1/x has a slope of
2e16 at the default tolerance.
This test can be fooled of course. The following straight line with large
enough slope does it.
[x, fval, info, out] = fzero (@(x) 1e8*x-pi, [-1,1], struct("TolX",1e-8))
which returns an info of -2 with the new patch. This seems okay to me
because the test was only meant to be a warning about singularities rather
than a definitive judgment. Moreover, there is a rather straightforward
test to clarify the matter. I can keep calling fzero with tighter
tolerances and if fzero continues to report the solution as a singularity
it very likely is. Using that technique
[x, fval, info, out] = fzero (@(x) 1e8*x-pi, [-1,1], struct("TolX",1e-9))
properly returns info=1 indicating that it has converged on a true zero.
The problem with the above setting is that low tolerances result in totally
unacceptable results. If I'm only looking for the rough location of a zero
I might use a tolX of 1e-2 which would result in a maximum slope of just 50
before classifying the result as a singularity. I can think of lots of
functions that have slopes above 50. Thus, for sloppy tolerances I put a
bound of 1e6 before classifying something as a singularity. This bound
won't catch tan(x) at a tolx of 1e-2, but just catches it at a tolx of
1e-3. I'll happily change this if you a better scale-dependent algorithm.
>
> Btw., Matlab 2007 misclassifies for instance this simple example:
>
>>> [x, fval, info, out] = fzero (@(x) sign(x) + (x==0), [-5,4])
Interesting, the new patched fzero catches this function.
--Rik
- Exit codes for fzero, Rik, 2010/02/10
- Re: Exit codes for fzero, Jaroslav Hajek, 2010/02/10
- Re: Exit codes for fzero, Rik, 2010/02/10
- Re: Exit codes for fzero, Jaroslav Hajek, 2010/02/11
- Re: Exit codes for fzero, Rik, 2010/02/15
- Re: Exit codes for fzero, Jaroslav Hajek, 2010/02/16
- Re: Exit codes for fzero,
Rik <=
- Re: Exit codes for fzero, Jaroslav Hajek, 2010/02/17
- Re: Exit codes for fzero, Rik, 2010/02/18
- Re: Exit codes for fzero, Jaroslav Hajek, 2010/02/19
Exit codes for fzero, John W. Eaton, 2010/02/10