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