calculating a function in matlab with very small values

Question

I am making a function in matlab to compute the following function:

for this function we have:

This is my implementation in matlab of the function:

function [b]= exponential(e)
%b = ? 
b= (exp (e) -1)/e;

When I test the function with very small values, indeed the limit of the function is to 1, but when the number is very small (eg 1*e-20) the limit goes to zero? what is the explanation to this phenomenon?. Am I doing something wrong?.

x= 10e-1 , f (x)= 1.0517

x= 10e-5 , f (x)=  1.0000

x= 10e-10 , f (x)=  1.0000

x= 10e-20 , f (x)=  0

Answer 1

The problem is that exp(x) is approx 1+x, but is being evaluated as 1 due to the 1 being indistinguishable from 1+x in the floating point representation. There is a MATLAB function expm1(x) (which is exp(x)-1 implemented for small x) that avoids the problem and works well for small arguments:

>> x=1e-100;expm1(x)/x
ans =
     1

Answer 2

I had to try my LIMEST tool out on it. As with any adaptive tool, it can be fooled, but it is usually pretty good.

fun = @(x) (exp(x) - 1)./x;

As you can see, fun has problems at zero.

fun(0)
ans =
   NaN

although if we evaluate fun near zero, we see it is near 1.

format long g
fun(1e-5)
ans =
          1.00000500000696

LIMEST succeeds, and even is able to provide an error estimate of the limit.

[lim,err] = limest(fun,0,'methodorder',3)

lim =
     1

err =
      2.50668568491927e-15

LIMEST uses a sequence of polynomial approximations, coupled with an adaptive Richardson extrapolation to generate both a limit estimate and a measure of the uncertainty on that limit.

So what problem are you seeing? The failure you saw is simple subtractive cancellation error. Thus, look at the value of

exp(1e-20)
ans =
     1

Even with format long g, the double precision value of exp(1e-20) is simply too close to 1 that when we subtract off 1, the result is an exact zero. Divide that by any non-zero value, and we get zero. Of course, when x is actually zero then we have a 0/0 condition, so NaN resulted when I tried that.

Lets see what happens in high precision. I'll use my HPF tool for that computation, and work in 64 decimal digits.

DefaultNumberOfDigits 64
exp(hpf('1e-20'))
ans =
1.000000000000000000010000000000000000000050000000000000000000166

See that when we sutract off 1, the difference between 1 and the exponential value is less than eps(1), so MATLAB must produce a zero value.

exp(hpf('1e-20')) - 1
ans =
1.000000000000000000005000000000000000000016666666666670000000000e-20

The question unasked is how I would choose to generate that function accurately in MATLAB. Clearly, you don't want to use brute force and define fun as I have, as you lose a great deal of accuracy. I'd probably just expand the Taylor series in a limited region around zero, and use fun as above for x significantly different from zero.

Answer 3

I think it has to do with the precision of your numbers. In short, the default precision for MATLAB is 5 digits. You can extend it to 15 with format long. Check out this article for more information about MATLAB precision