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
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
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.
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
format only deals with how many digits are displayed as output to the console. It doesn't change the underlying number - tmpearce 2012-04-04 01:11