Thread Subject: x-> 0 limit problem in MATLAB?

Hi,
I recently tried to evaluate the function:

f=(x.*exp(x)-exp(x)+1)./x.^2

for values very close to zero. Typically, whenever the argument is
zero, I replace it with eps to avoid the singularity in the
denominator.

It can be shown that the limit of f for x->0 is 0.5. Indeed, for small
values of x, f correctly approaches 0.5.

E.g:
x=1e-6;f=(x.*exp(x)-exp(x)+1)./x.^2

f =

    0.5000

However, something strange happens when I make x even smaller:

x=1e-7;f=(x.*exp(x)-exp(x)+1)./x.^2

f =

    0.5107

And ultimately:

x=eps;f=(x.*exp(x)-exp(x)+1)./x.^2

f =

     0

Which is clearly wrong

Even worse, when I slightly rewrite the same function as:

f=(x.*exp(x)-(exp(x)-1))./x.^2

I get:

x=eps;f=(x.*exp(x)-(exp(x)-1))./x.^2

f =

     1

There clearly is a problem with approaching the limit for small values
of x.
Is there any way to rewrite the function that would avoid this issue?

The MATLAB version number is 6.1.0.450 (R12.1)


Thanks,

Sigurd

This is a problem with floating point math.

f(10^-7) = -99999990000000*exp(1/10000000)+100000000000000
~ -9.999999999999949999996666666542*10^13 + 100000000000000

You are subtracting two very large numbers that differ only
slightly. Floating point doesn't handle this very well.
For this kind of problem you either have to define the
function to be the limit within a certain tolerance or use
Maple.

sigurd.ss@gmail.com wrote in message <42f8a1c3-6680-4f98-
aa85-86066016d616@p25g2000hsf.googlegroups.com>...
> Hi,
> I recently tried to evaluate the function:
>
> f=(x.*exp(x)-exp(x)+1)./x.^2
> ...... (snip) ........
> There clearly is a problem with approaching the limit for small values
> of x.
> Is there any way to rewrite the function that would avoid this issue?
> Sigurd
------------
  If your f(x) is something you need to compute frequently with small values
of x, there is a way around your difficulties, but it requires some considerable
effort.

  If you expand your particular f(x) in an infinite power series of powers of x,
you get

 f(x) = 1/2! + 2/3!*x + 3/4!*x^2 + 4/5!*x^3 + 5/6!*x^4 + .....

This shows that you can use appropriate polynomials to accurately
approximate f(x) for all small values of x, and your original form for larger
x's, thereby avoiding computational problems. However, as you can imagine,
it requires a careful error analysis as to how high a degree polynomial is
necessary and at what point the transition in form is to be made.

  As Matt has pointed out, the problems in your original form of f(x) with very
small values of x are caused by the subtraction of two sizable quantities
whose difference is so small that you have lost a substantial amount of
accuracy. After all, all computers have only a finite number of digits to use in
expressing their numbers, and such a subtraction can suddenly create very
much larger errors than are present in ordinary computation. This is not just
a problem with matlab. It affects all numerical computation methods. In
matlab's case their binary floating point numbers have 53 bits of significand
for an accuracy of one part in 2^53, or about 16 decimal places, which is a
very high accuracy in most respects. But if you subtract two numbers which
are, say, near 1, and yet differ only out in their 16th place, your answer
becomes completely obscured by the consequent round off errors.

Roger

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in
message <g18737$8jm$1@fred.mathworks.com>...
> sigurd.ss@gmail.com wrote in message <42f8a1c3-6680-4f98-
> aa85-86066016d616@p25g2000hsf.googlegroups.com>...
> > Hi,
> > I recently tried to evaluate the function:
> >
> > f=(x.*exp(x)-exp(x)+1)./x.^2
> > ...... (snip) ........
> > There clearly is a problem with approaching the limit for small values
> > of x.
> > Is there any way to rewrite the function that would avoid this issue?
> > Sigurd
> ------------
> If your f(x) is something you need to compute frequently with small values
> of x, there is a way around your difficulties, but it requires some
considerable
> effort.
>
> If you expand your particular f(x) in an infinite power series of powers of
x,
> you get
>
> f(x) = 1/2! + 2/3!*x + 3/4!*x^2 + 4/5!*x^3 + 5/6!*x^4 + .....
>
> This shows that you can use appropriate polynomials to accurately
> approximate f(x) for all small values of x, and your original form for larger
> x's, thereby avoiding computational problems. However, as you can
imagine,
> it requires a careful error analysis as to how high a degree polynomial is
> necessary and at what point the transition in form is to be made.
>
> As Matt has pointed out, the problems in your original form of f(x) with
very
> small values of x are caused by the subtraction of two sizable quantities
> whose difference is so small that you have lost a substantial amount of
> accuracy. After all, all computers have only a finite number of digits to use
in
> expressing their numbers, and such a subtraction can suddenly create very
> much larger errors than are present in ordinary computation. This is not
just
> a problem with matlab. It affects all numerical computation methods. In
> matlab's case their binary floating point numbers have 53 bits of
significand
> for an accuracy of one part in 2^53, or about 16 decimal places, which is a
> very high accuracy in most respects. But if you subtract two numbers
which
> are, say, near 1, and yet differ only out in their 16th place, your answer
> becomes completely obscured by the consequent round off errors.

Just to expand on Roger's comments...

This is called subtractive cancellation.
In your expression, think about what
happens.

(x.*exp(x)-exp(x)+1)./x.^2

if you go back to the Taylor series, then
you should see it. Since the Taylor series
for exp(x) is just

1+x+x^2/2+...+x^n/factorial(n)+...

if we substitute this into the expression,
we see that not only does the constant
term drop off, but so does the first order
term. What remains is only the second
order term and above. Then you divide
by x^2 at the very end to yield the series
that Roger reports.

So anything under sqrt(eps) in this
expression will produce pure garbage,
and even as we approach sqrt(eps),
expect it to be nasty. Try it out.

fun = @(x) (x.*exp(x)-exp(x)+1)./x.^2;
format long g

fun(1e-2)
ans =
         0.503345866736948

fun(1e-3)
ans =
         0.500333458330893
fun(1e-4)
ans =
         0.500033325856464

fun(1e-5)
ans =
         0.500004482262284

fun(1e-6)
ans =
         0.500044450291171

fun(1e-7)
ans =
         0.510702591327572

fun(1e-8)
ans =
          1.11022302462516

fun(1e-9)
ans =
     0

See what happened? As we get near the
point where the subtractive cancellation
starts to dominate, things go haywire.

Just for kicks, I had to throw this function
at my own utility for extracting general
limits.

[lim,err] = limest(fun,0,'MethodOrder',5)
lim =
         0.500000000165239
err =
       9.2459501732486e-10

It gets the correct result, considerably
more accurate than any simple function
value, but even so, the error estimate is
still not much smaller than sqrt(eps).

John

Just a minor comment on John's reply:

fun = @(x) (x.*exp(x)-exp(x)+1)./x.^2;

might not make much sense to the OP since he is using ver 6.1!

"Matt Fig" <spamanon@yahoo.com> wrote in message
<g195rl$3q4$1@fred.mathworks.com>...
> Just a minor comment on John's reply:
>
> fun = @(x) (x.*exp(x)-exp(x)+1)./x.^2;
>
> might not make much sense to the OP since he is using ver 6.1!

Oops. Correct. Use an inline function instead.

fun = inline('(x.*exp(x)-exp(x)+1)./x.^2','x');

John

sigurd.ss@gmail.com wrote in message
> It can be shown that the limit of f for x->0 is 0.5.
Indeed, for small
> values of x, f correctly approaches 0.5.
Hello Sigurd:
I tried to figure out if you need a general solution to
this, or just a particular solution for this expression.
In any case, you can use the L'Hopital rule in cases like this.
http://en.wikipedia.org/wiki/L'H%C3%B4pital's_rule
So the limit you are trying to calculate can be
appropriately evaluated as
(derivative of the numerator)/(derivative of the denominator)
which in this case simply leads to
exp(x)/2
which does not suffer to "catastrophic cancellation" for any
x, but is valid when x is close to zero.

The best answer is strongly connected with your needs.
Another solution, provided you know the shape of the leading
term of the Taylor expansion, is to use Richardson
Extrapolation
(http://en.wikipedia.org/wiki/Richardson_extrapolation)
which can be safely applied using values far from the
catastrophic area.
Why are you further interested in finding a limit you
already know? If you describe better the context of your
problem maybe we can provide other choices.
Regards
Carlos

This message is intentionally left blank

"carlos lopez" <clv2clv_00000000_@adinet.com.uy> wrote in message
<g1d32e$l46$1@fred.mathworks.com>...
> Hello Sigurd:
> I tried to figure out if you need a general solution to
> this, or just a particular solution for this expression.
> In any case, you can use the L'Hopital rule in cases like this.
> http://en.wikipedia.org/wiki/L'H%C3%B4pital's_rule
> So the limit you are trying to calculate can be
> appropriately evaluated as
> (derivative of the numerator)/(derivative of the denominator)
> which in this case simply leads to
> exp(x)/2
> which does not suffer to "catastrophic cancellation" for any
> x, but is valid when x is close to zero.
>
> The best answer is strongly connected with your needs.
> Another solution, provided you know the shape of the leading
> term of the Taylor expansion, is to use Richardson
> Extrapolation
> (http://en.wikipedia.org/wiki/Richardson_extrapolation)
> which can be safely applied using values far from the
> catastrophic area.
> Why are you further interested in finding a limit you
> already know? If you describe better the context of your
> problem maybe we can provide other choices.
> Regards
> Carlos
--------------
Hello Carlos,

  Sigurd stated that he wanted to "evaluate the function ... for values very
close to zero" which is asking somewhat more than simply finding its limit as
x approaches 0. He specifically asked for a way to "rewrite the function" so as
to avoid the trouble near 0. That is the reason I showed him a power series
expansion about 0. Otherwise there would have been no point in doing so.
By the way, that expansion shows that exp(x)/2 would not be a very good
approximation near x = 0. It comes into zero with the wrong slope.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
> Hello Carlos,
>
> Sigurd stated that he wanted to "evaluate the function
... for values very
> close to zero" which is asking somewhat more than simply
finding its limit as
> x approaches 0. He specifically asked for a way to
"rewrite the function" so as
> to avoid the trouble near 0. That is the reason I showed
him a power series
> expansion about 0. Otherwise there would have been no
point in doing so.
> By the way, that expansion shows that exp(x)/2 would not
be a very good
> approximation near x = 0. It comes into zero with the
wrong slope.
>
> Roger Stafford
As usual you are right...
Indeed this problem suffers from cancellation twice:
.- when evaluating exp(x)-1, which has a reminder of the
order of x
.- when substracting x*exp(x) from the above, because it
also cancels out the x term
For the former I feel that he can achieve a modest
improvement (a factor of two, i.e. one bit) by using the
following identity:
(1-exp(x))=(1-exp(2*x))/(1+exp(x))
because the denominator does not suffer from cancellation
and the numerator uses a larger exponent (twice). The
identity can be used more than once, producing for example
(1-exp(x))=(1-exp(4*x))/(1+exp(x))/(1+exp(2x))
keeping the cancellation effect only in the numerator, but
every time with a larger exponent (which remediates modestly
the problem).
Notice that the terms are all powers of exp(x), so there is
no need for extensive calculation of exponentials.
If Sigurd needs a simpler solution, he might want to
consider using higher-than-double precision in the
calculations. The Multiprecision toolbox (by Ben Barrowes)
achieves that whithout manipulating the expression but just
changing the machine epsilon at will.
I could not figure out which is the context where he needs
this for a particular function; maybe he es seeking for a
general solution for other functions. In any case, your
approach using Taylor expansion is always valid, but
requires a little effort to do and it has always an exact
truncation term. Richardson extrapolation has also the same
drawback, because it requires to know the order of the first
non-zero term of the expansion.
Regards
Carlos

"carlos lopez" <clv2clv_00000000_@adinet.com.uy> wrote in message
<g1e875$dro$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid>
> > Hello Carlos,
> >
> > Sigurd stated that he wanted to "evaluate the function
> ... for values very
> > close to zero" which is asking somewhat more than simply
> finding its limit as
> > x approaches 0. He specifically asked for a way to
> "rewrite the function" so as
> > to avoid the trouble near 0. That is the reason I showed
> him a power series
> > expansion about 0. Otherwise there would have been no
> point in doing so.
> > By the way, that expansion shows that exp(x)/2 would not
> be a very good
> > approximation near x = 0. It comes into zero with the
> wrong slope.
> >
> > Roger Stafford
> As usual you are right...
> Indeed this problem suffers from cancellation twice:
> .- when evaluating exp(x)-1, which has a reminder of the
> order of x
> .- when substracting x*exp(x) from the above, because it
> also cancels out the x term
> For the former I feel that he can achieve a modest
> improvement (a factor of two, i.e. one bit) by using the
> following identity:
> (1-exp(x))=(1-exp(2*x))/(1+exp(x))
> because the denominator does not suffer from cancellation
> and the numerator uses a larger exponent (twice). The
> identity can be used more than once, producing for example
> (1-exp(x))=(1-exp(4*x))/(1+exp(x))/(1+exp(2x))
> keeping the cancellation effect only in the numerator, but
> every time with a larger exponent (which remediates modestly
> the problem).
> Notice that the terms are all powers of exp(x), so there is
> no need for extensive calculation of exponentials.
> If Sigurd needs a simpler solution, he might want to
> consider using higher-than-double precision in the
> calculations. The Multiprecision toolbox (by Ben Barrowes)
> achieves that whithout manipulating the expression but just
> changing the machine epsilon at will.
> I could not figure out which is the context where he needs
> this for a particular function; maybe he es seeking for a
> general solution for other functions. In any case, your
> approach using Taylor expansion is always valid, but
> requires a little effort to do and it has always an exact
> truncation term. Richardson extrapolation has also the same
> drawback, because it requires to know the order of the first
> non-zero term of the expansion.
> Regards
> Carlos

I decided to submit my limest code to the FEX.
It might not appear until tomorrow given that
today is a holiday. It uses a family of polynomial
extrapolants to build a good approximation to
the limit.

fun = @(x) (x.*exp(x)-exp(x)+1)./x.^2;
[lim,err] = limest(fun,0)

lim =
         0.500000000094929
err =
      1.56957823156003e-09

For a little more accuracy and a bit more
work, a higher order method is available.

[lim,err] = limest(fun,0,'methodorder',8,'step',1.25);

lim
lim =
         0.500000000002992
 
John

On May 26, 5:17 am, "John D'Errico" <woodch...@rochester.rr.com>
wrote:
> "carlos lopez" <clv2clv_000000...@adinet.com.uy> wrote in message
>
> <g1e875$dr...@fred.mathworks.com>...
>
>
>
> > "Roger Stafford" <ellieandrogerxy...@mindspring.com.invalid>
> > > Hello Carlos,
>
> > > Sigurd stated that he wanted to "evaluate the function
> > ... for values very
> > > close to zero" which is asking somewhat more than simply
> > finding itslimitas
> > > x approaches 0. He specifically asked for a way to
> > "rewrite the function" so as
> > > to avoid the trouble near 0. That is the reason I showed
> > him a power series
> > > expansion about 0. Otherwise there would have been no
> > point in doing so.
> > > By the way, that expansion shows that exp(x)/2 would not
> > be a very good
> > > approximation near x = 0. It comes into zero with the
> > wrong slope.
>
> > > Roger Stafford
> > As usual you are right...
> > Indeed this problem suffers from cancellation twice:
> > .- when evaluating exp(x)-1, which has a reminder of the
> > order of x
> > .- when substracting x*exp(x) from the above, because it
> > also cancels out the x term
> > For the former I feel that he can achieve a modest
> > improvement (a factor of two, i.e. one bit) by using the
> > following identity:
> > (1-exp(x))=(1-exp(2*x))/(1+exp(x))
> > because the denominator does not suffer from cancellation
> > and the numerator uses a larger exponent (twice). The
> > identity can be used more than once, producing for example
> > (1-exp(x))=(1-exp(4*x))/(1+exp(x))/(1+exp(2x))
> > keeping the cancellation effect only in the numerator, but
> > every time with a larger exponent (which remediates modestly
> > the problem).
> > Notice that the terms are all powers of exp(x), so there is
> > no need for extensive calculation of exponentials.
> > If Sigurd needs a simpler solution, he might want to
> > consider using higher-than-double precision in the
> > calculations. The Multiprecision toolbox (by Ben Barrowes)
> > achieves that whithout manipulating the expression but just
> > changing the machine epsilon at will.
> > I could not figure out which is the context where he needs
> > this for a particular function; maybe he es seeking for a
> > general solution for other functions. In any case, your
> > approach using Taylor expansion is always valid, but
> > requires a little effort to do and it has always an exact
> > truncation term. Richardson extrapolation has also the same
> > drawback, because it requires to know the order of the first
> > non-zero term of the expansion.
> > Regards
> > Carlos
>
> I decided to submit my limest code to the FEX.
> It might not appear until tomorrow given that
> today is a holiday. It uses a family of polynomial
> extrapolants to build a good approximation to
> thelimit.
>
> fun = @(x) (x.*exp(x)-exp(x)+1)./x.^2;
> [lim,err] = limest(fun,0)
>
> lim =
> 0.500000000094929
> err =
> 1.56957823156003e-09
>
> For a little more accuracy and a bit more
> work, a higher order method is available.
>
> [lim,err] = limest(fun,0,'methodorder',8,'step',1.25);
>
> lim
> lim =
> 0.500000000002992
>
> John

All,
thanks for the comments so far.

I realize this is an issue of numerical instability, but I'm just
surprised MATLAB would give no hint or warning that the result may be
unreliable (as it does in other cases). I thought when evaluating a
function at the built-in value "eps", results would be consistently
reliable - apparently not. For other functions (e.g. sin(x)/x) I
always got a good approximation of the actual limit by working this
way.
I could use the Taylor expansion for small x, but I'm not specifically
interested in small values only. x=0 (or eps) is just one of the
values and that's where I happened to notice some weird behavior when
I plotted the function.

I'll certainly lake a look at the "limest" function.

regards,

Sigurd

sigurd.ss@gmail.com wrote:
> Hi,
> I recently tried to evaluate the function:
>
> f=(x.*exp(x)-exp(x)+1)./x.^2

Sigurd, you don't say how accurately you need this function, or how small your
arguments really need to be. Everything that everyone else has said is right,
but you might also look at the following very simple modification to your
expression, using EXPM1:

x = logspace(0,-14,10000);
limit = .5;

% original expression
f1 = (x.*exp(x)-exp(x)+1)./x.^2;
subplot(3,1,1), loglog(x,f1-limit,'-');

% use expm1, a bit better
f2 = (x.*exp(x)-expm1(x))./x.^2;
subplot(3,1,2), loglog(x,f2-limit,'-');

% use truncated taylor series for very small values
f3 = (x.*exp(x)-expm1(x))./x.^2;
small = (x<1e-5); % choose cutoff to trade off between roundoff noise and bias
f3(small) = (1/2) + (1/3)*x(small) + (3/8)*x(small).^2;
subplot(3,1,3), loglog(x,f3-limit,'-');

Peter Perkins <Peter.PerkinsRemoveThis@mathworks.com> wrote in message
<g1hha4$ghe$1@fred.mathworks.com>...
> .........
> % use truncated taylor series for very small values
> f3 = (x.*exp(x)-expm1(x))./x.^2;
> small = (x<1e-5); % choose cutoff to trade off between roundoff noise and
bias
> f3(small) = (1/2) + (1/3)*x(small) + (3/8)*x(small).^2;
> subplot(3,1,3), loglog(x,f3-limit,'-');
---------------
Hello Peter,

  Minor correction: The coefficient for x^2 should be 1/8, not 3/8, in the Taylor
expansion of Sigurd's function.

Roger Stafford

Roger Stafford wrote:

> Minor correction: The coefficient for x^2 should be 1/8, not 3/8, in the Taylor
> expansion of Sigurd's function.

Yes, thanks (I should have just read your first post instead of doing the
cipherin' myself).

Peter Perkins <Peter.PerkinsRemoveThis@mathworks.com> wrote
> % use truncated taylor series for very small values
> f3 = (x.*exp(x)-expm1(x))./x.^2;
> small = (x<1e-5); % choose cutoff to trade off between
roundoff noise and bias
> f3(small) = (1/2) + (1/3)*x(small) + (3/8)*x(small).^2;
I might suggest a small improvement: use the derivative
values at x=small instead of x=0 in order to have higher
order continuity (up to the second derivative)
Regards
Carlos

On May 27, 12:47 pm, sigurd...@gmail.com wrote:
-----SNIP
> All,
> thanks for the comments so far.
>
> I realize this is an issue of numerical instability, but I'm just
> surprised MATLAB would give no hint or warning that the result may be
> unreliable (as it does in other cases). I thought when evaluating a
> function at the built-in value "eps", results would be consistently
> reliable - apparently not. For other functions (e.g. sin(x)/x) I
> always got a good approximation of the actual limit by working this
> way. I could use the Taylor expansion for small x,

You don't seem to be getting the point.

MATLAB yields 1 + z = 1 for z < eps.

Therefore, as x --> 0

sin(x) --> x -x^3/6
         = x*(1-x^2/6)
         = x for x < sqrt(6*eps) = 3.6500e-008.

Therefore, MATLAB should yield

(x-sin(x))./x.^3 = 0 for x < sqrt(6*eps)

instead of 1/6.

Consider

 >> format short e
 x = 10.^(-(1:20)');

 result = [ x, x.^3/6, x - x.*(1-x.^2/6), x- (x-x.^3/6), x - sin(x) ]

result =

  1.0000e-001 1.6667e-004 1.6667e-004 1.6667e-004 1.6658e-004
  1.0000e-002 1.6667e-007 1.6667e-007 1.6667e-007 1.6667e-007
  1.0000e-003 1.6667e-010 1.6667e-010 1.6667e-010 1.6667e-010
  1.0000e-004 1.6667e-013 1.6667e-013 1.6667e-013 1.6667e-013
  1.0000e-005 1.6667e-016 1.6667e-016 1.6667e-016 1.6667e-016
  1.0000e-006 1.6667e-019 1.6665e-019 1.6665e-019 1.6665e-019
  1.0000e-007 1.6667e-022 1.7205e-022 1.7205e-022 1.7205e-022
  1.0000e-008 1.6667e-025 0 0 0
  1.0000e-009 1.6667e-028 0 0 0
  1.0000e-010 1.6667e-031 0 0 0
  1.0000e-011 1.6667e-034 0 0 0
  1.0000e-012 1.6667e-037 0 0 0
  1.0000e-013 1.6667e-040 0 0 0
  1.0000e-014 1.6667e-043 0 0 0
  1.0000e-015 1.6667e-046 0 0 0
  1.0000e-016 1.6667e-049 0 0 0
  1.0000e-017 1.6667e-052 0 0 0
  1.0000e-018 1.6667e-055 0 0 0
  1.0000e-019 1.6667e-058 0 0 0
  1.0000e-020 1.6667e-061 0 0 0


  Therefore blind use of series expansions are not the answer.


For your original function, it is easier to see what is happening
by factoring out exp(x) to obtain

f = exp(x).*( exp(-x) -1 + x )./x.^2

Since the expression in parentheses --> 0.5*x.^2
as x --> 0, consider

format short e
x = 10.^(-(1:20)');
result = [ x, 0.5*x.^2, (1 - x + 0.5*x.^2)-(1 - x), exp(-x)-(1-x)]

>>
result =

  1.0000e-001 5.0000e-003 5.0000e-003 4.8374e-003
  1.0000e-002 5.0000e-005 5.0000e-005 4.9834e-005
  1.0000e-003 5.0000e-007 5.0000e-007 4.9983e-007
  1.0000e-004 5.0000e-009 5.0000e-009 4.9998e-009
  1.0000e-005 5.0000e-011 5.0000e-011 5.0000e-011
  1.0000e-006 5.0000e-013 5.0004e-013 5.0004e-013
  1.0000e-007 5.0000e-015 4.9960e-015 4.9960e-015
  1.0000e-008 5.0000e-017 0 1.1102e-016
  1.0000e-009 5.0000e-019 0 0
  1.0000e-010 5.0000e-021 0 0
  1.0000e-011 5.0000e-023 0 0
  1.0000e-012 5.0000e-025 0 0
  1.0000e-013 5.0000e-027 0 0
  1.0000e-014 5.0000e-029 0 0
  1.0000e-015 5.0000e-031 0 0
  1.0000e-016 5.0000e-033 0 0
  1.0000e-017 5.0000e-035 0 0
  1.0000e-018 5.0000e-037 0 0
  1.0000e-019 5.0000e-039 0 0
  1.0000e-020 5.0000e-041 0 0

  Hope this helps.

  Greg

sigurd.ss@gmail.com wrote in message
> I realize this is an issue of numerical instability, but
I'm just
> surprised MATLAB would give no hint or warning that the
result may be
> unreliable (as it does in other cases). I thought when
evaluating a
> function at the built-in value "eps", results would be
consistently
> reliable - apparently not.
Hello Sigurd:
I think that MATLAB does not provide any warning because
IEEE748 might not request so.
If you want to consider ways to extend MATLAB to do so, you
need to modify the behavior of the addition and substraction
operation, requesting a comparison. For example, you have to
create a directory @double below any directory in the path,
and store there a file minus.m which will look similar to this:
function c=minus(a,b)
%Minus operation between matrices a and b with warning about
%shift out and catastrophic cancellation
%UNTESTED VERSION 20080528

%check where a has non-zero elements
q=find(a);
if any(b(q)./a(q) <= eps)
   %Issue a warning, but perform the expected operation
   %for backward compatibility
   warning('MATLAB:ShiftOut','Shift-out or catastrophic
cancellation in progress')
end
c=builtin('minus',a,b);

Something similar should be performed for the plus.m file. I
haven't worried about checking the sign of the terms (it
should be done) but hope it illustrates the idea.
Solutions like these come with a price: every time there is
an addition or substraction a quotient and a comparison
should be carried out. The overhead is too heavy to be
considered in normal cases, and that's why it has not been
implemented widely.
In any case you might want to enable/disable this feature at
your will.
Once you create the directory and have the files there, to
enable it you should issue once the command:
rehash toolboxcache
If you want to disable this operation with warning, rename
the directory to something like @double.inactive, and rehash
again (there might be other simpler alternatives).
I expect to hear other opinions.
Regards
Carlos

"carlos lopez" <clv2clv_00000000_@adinet.com.uy> wrote in message
<g1jk6a$joc$1@fred.mathworks.com>...
> sigurd.ss@gmail.com wrote in message
> > I realize this is an issue of numerical instability, but
> I'm just
> > surprised MATLAB would give no hint or warning that the
> result may be
> > unreliable (as it does in other cases). I thought when
> evaluating a
> > function at the built-in value "eps", results would be
> consistently
> > reliable - apparently not.
> Hello Sigurd:
> I think that MATLAB does not provide any warning because
> IEEE748 might not request so.

But do you think that most users of Matlab
would want to see a warning message come
out EVERY time they do something simple?

x = (1:10) - cos(2*pi);

would result in something like

WARNING! You may be committing an act of
subtractive cancellation here. Are you really
positive that you know what you are doing?
Perhaps you should call our free help line for
the numerically challenged: 1-800-num-anal.

(Gosh, I hate to think what that number does
reach. ;-) I hope it is for numerical analysis.)

The point is, you would not want warning
messages for every possible theoretical
numerical infraction, at least unless you
are a lawyer. In that case, the lawyers would
have the Mathworks generate a long warning
message after EVERY line of code. Some of
the messages would tell us of possible
subtractive cancellation, some would tell us
that our call to inv is potentially suspect, or
that all calls to an adaptive quadrature
routine may result in an incorrect result,
or that calling an optimizer may not
converge unless we have adequate starting
values.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in
> But do you think that most users of Matlab
> would want to see a warning message come
> out EVERY time they do something simple?
>
> x = (1:10) - cos(2*pi);
>
> would result in something like
>
> WARNING! You may be committing an act of
> subtractive cancellation here. Are you really
> positive that you know what you are doing?
> Perhaps you should call our free help line for
> the numerically challenged: 1-800-num-anal.
>
> (Gosh, I hate to think what that number does
> reach. ;-) I hope it is for numerical analysis.)
>
> The point is, you would not want warning
> messages for every possible theoretical
> numerical infraction, at least unless you
> are a lawyer. In that case, the lawyers would
> have the Mathworks generate a long warning
> message after EVERY line of code. Some of
> the messages would tell us of possible
> subtractive cancellation, some would tell us
> that our call to inv is potentially suspect, or
> that all calls to an adaptive quadrature
> routine may result in an incorrect result,
> or that calling an optimizer may not
> converge unless we have adequate starting
> values.
>
> John
:-) You provided all the necessary arguments to explain why
this is NOT the way it is implemented (neither in Matlab nor
anywhere else I am aware of); hope that Sigurd will read
your answer.
Anyway, Sigurd asked for a "warning tool", so I tried to
suggest one. It might be an interesting idea, at least for
some narrow purposes... like locating the offending line in
a long code.
Regards
Carlos
PD: I always love to look at your comments.