Correctly evaluate exponential function for limit values

Hi all,

I have a function, of the form

whereby g(x) and h(x) are similar functions that I have expressions for, and x is a 2D meshgrid of values. For certain x values, g(x) = 0 and/or h(x) = 0. At these points, the half with the exponent function equal to 0 should evaluate to N (from Taylor Series expansion i.e. if both g(x) and h(x) = 0 for same x value). However, MATLAB doesn't know this limit and so evaluates the function as NaN for that value of x.

Is there an elegant way to for MATLAB to correctly evaluate this function at these points (other than just including many if statements inside the function)?

Thank you in advance.

Answer 1

[Torsten]

function result = f(x,g,h,N)
  result = N^2*ones(size(x));
  K = g(x).*h(x) ~= 0;
  if (any (K(:)))
    result(K) = your expression for f;
  end
end

[Dog_Biscuit]

Thank you :)

[James Tursa]

You may want more elaborate checking than this. E.g.,

If only one of g(x) or h(x) is 0, then you will not get for a limiting result. You will get N times the other terms.

If either g(x) or h(x) is very near 0 but not exactly 0, then you still might get a 0/0 (NaN) result because exp(stuff) could evaluate to be exactly 1 even though stuff is nonzero. In this case you could use a truncated Taylor series to evaluate the 1-exp(stuff) as (-stuff) or perhaps (-stuff-stuff^2/2) etc. E.g.,

exp(1e-20)

ans = 1

So 1-exp(1e-20) will evaluate to be exactly 0 even though the argument to exp( ) is nonzero. But using a truncated Taylor series for these cases will easily get you the 1e-20 result that is needed to avoid the NaN. The code would be something like this with an appropriate tolerance:

if( abs(stuff) < tolerance )
    result = -stuff - stuff^2/2;
else
    result = 1 - exp(stuff);
end

Bottom line is how robust do you want this code to be? If you want it to be bulletproof and return limiting values and avoid NaN in these cases then you will need to check the exp( ) arguments individually and have more elaborate if-then-else function evaluation code.

[Torsten]

You are right - the solution given is not correct.

Should be replaced by

function result_gh = f(x,g,h,N)
  result_g = N*ones(size(x));
  result_h = result_g;
  K_g = g(x)~=0;
  if (any(K_g(:))
    result_g(K_g) = the quotient with g(x)
  end
  K_h = h(x)~= 0;
  if (any(K_h(:))
    result_h(K_h) = the quotient with h(x)  
  end
  result_gh = result_g.*result_h;
end