CODEWITHSUNDEEP
matlab
Ka Wa Yip
Ka Wa Yip

Reputation: 2993

MATLAB: very small number divided by very small number

I have a very small number divided by very small number and MATLAB gives me Inf (when x is too small) or NaN (when x is zero). But the correct limit value should be 2. How to solve this problem? To get j as 2 for too small x or x = 0?

x = 0.000001 

>> x = 0.000001
x =
   1.0000e-06
>> j = 2*x*(1+(-1)*exp(-x))^(-1)*exp(-x)
j =
    2.0000

x = 0.01

>> x = 0.01
x =
    0.0100
>> j = 2*x*(1+(-1)*exp(-x))^(-1)*exp(-x)
j =
    1.9900

x = 1e-19

>> x = 1e-19
x =
   1.0000e-19
>> j = 2*x*(1+(-1)*exp(-x))^(-1)*exp(-x)
j =
   Inf

x = 0

>> x = 0
x =
     0
>> j = 2*x*(1+(-1)*exp(-x))^(-1)*exp(-x)
j =
   NaN

Upvotes: 2

Views: 1367

Answers (2)

nirvana-msu
nirvana-msu

Reputation: 4077

The fact that you get NaN for x=0 is correct because you cannot divide by zero. The limit is only valid when x->0, not at x=0.

It does work for up to x ~= 10^-12, the problem is just that you quickly start loosing precision when the exponent approaches 1 due to specifics of floating point arithmetics:

f = @(x) 2*x*exp(-x)/(1-exp(-x));
pow = (-1)*(1:19);
X = 10.^pow; % 10^-1 to 10^-19
[pow', arrayfun(f, X)']

ans =
   -1.0000    1.9017
   -2.0000    1.9900
   -3.0000    1.9990
   -4.0000    1.9999
   -5.0000    2.0000
   -6.0000    2.0000
   -7.0000    2.0000
   -8.0000    2.0000
   -9.0000    2.0000
  -10.0000    2.0000
  -11.0000    2.0000
  -12.0000    2.0000
  -13.0000    1.9994
  -14.0000    2.0016
  -15.0000    2.0016
  -16.0000    1.8014
  -17.0000       Inf
  -18.0000       Inf
  -19.0000       Inf

One way to get around that would be to use symbolic math, as suggested in the other answer. If you want to do it numerically, you would have to find the trade-off between evaluating at the point too far from limiting point where the limit is not yet reached, and evaluating too close to the limiting point where numerical errors overwhelm the result.

There's a great Adaptive numerical limit (and residue) estimation submission on FEX that tries to achieve this. The approach it takes is somewhat similar to Richardson extrapolation method. In particular:

  1. The limit is evaluated at a sequence of geometrically spaced points approaching the limiting point.
  2. Polynomial is used to fit the resulting sequence of values. The trick is to find the balance of the points chosen in step 1. The closer they are to the limiting point, the more distorted the polynomial fit becomes.
  3. Take the constant term of the polynomial model as the limit.

In your particular case that utility works like a charm without any need to tweak the sequence of points. Note that it expects your function to be vectorized, i.e. being able to evaluate result at a vector of input points, so use element-wise operators when defining it:

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

lim =
    2.0000
err =
   5.5161e-12

Upvotes: 2

rahnema1
rahnema1

Reputation: 15837

you may use symbolic math toolbox :

syms x;
limit(2*x*(1+(-1)*exp(-x))^(-1)*exp(-x), 0)

Upvotes: 1

Related Questions