Thread Subject: Help! Numerically handle removable singularities: a more clear

Hi all,

I am helping my friend on this issue.

We have some numerical computation which involves a lot of "log(g)*g"
and we have to evaluate these for g on the interval [0, 1],
numerically... The reason that 0 is included is that because the
singularity at 0 is removable.


In Maple:

> x := 0; log(x)*x;
%;
                                      0
Error, (in ln) numeric exception: division by zero

In Matlab:


>> x=0; log(x)*x

ans =

   NaN


Our procedure is that we first derive very complicated symbolic
expressions in Maple and then translate into Matlab to get numerical
calculation.

How do handle the above "NaN" in Matlab? This is just one intermediate
step using "log(x)*x" (this expression is from Maple symbolic
results).


Of course we could take a limit at x=0 and handle that case in Matlab
separately, let say when |x| < 1e-5, we switch to this limiting
result; and when |x|>=1e-5, we still use the numerical approach.


But this creates a discontinuity at |x|=1e-5 in numerical results in
Matlab.


Could anybody shed some lights on us about how to handle such a
situation?


Thanks a lot!

On Sat, 14 Nov 2009 13:52:56 -0800 (PST), Luna Moon
<lunamoonmoon@gmail.com> wrote:

>Hi all,
>
>I am helping my friend on this issue.
>
>We have some numerical computation which involves a lot of "log(g)*g"
>and we have to evaluate these for g on the interval [0, 1],
>numerically... The reason that 0 is included is that because the
>singularity at 0 is removable.
  

http://ideas.repec.org/c/boc/bocode/x031701.html

http://fmwww.bc.edu/repec/bocode/s/slog.ox

A.L.

Luna Moon wrote:
> Hi all,
>
> I am helping my friend on this issue.
>
> We have some numerical computation which involves a lot of "log(g)*g"
> and we have to evaluate these for g on the interval [0, 1],
> numerically... The reason that 0 is included is that because the
> singularity at 0 is removable.
>
>
> In Maple:
>
>> x := 0; log(x)*x;
> %;
> 0
> Error, (in ln) numeric exception: division by zero
>
> In Matlab:
>
>
>>> x=0; log(x)*x
>
> ans =
>
> NaN
>
>
> Our procedure is that we first derive very complicated symbolic
> expressions in Maple and then translate into Matlab to get numerical
> calculation.
>
> How do handle the above "NaN" in Matlab? This is just one intermediate
> step using "log(x)*x" (this expression is from Maple symbolic
> results).
>
>
> Of course we could take a limit at x=0 and handle that case in Matlab
> separately, let say when |x| < 1e-5, we switch to this limiting
> result; and when |x|>=1e-5, we still use the numerical approach.
>
>
> But this creates a discontinuity at |x|=1e-5 in numerical results in
> Matlab.
>
>
> Could anybody shed some lights on us about how to handle such a
> situation?
>
>
> Thanks a lot!

log(x)*x when x = 0 is undefined because log(x) is undefined

limit(log(x)*x,x=0) = 0

Generally if you are simply doing a numerical computation you could handle
it as a side case such as

if (x == 0) \\ or |x| <= 0.00000001 or something similar
return 0;

or

// hack that is not generally a good way
if (x==0) \\ or |x| <= 0.00000001
 x = 0.00001
// do comp

If your doing some advanced numerical method then generally you will find
some numerial appoximation to log(x)*x in the first place where the
discontinuity is removed such as a taylor series


log(x)*x ~= -sum((1 - x)^k/k)*x
(not a good method)

or by using a faster converging series such as some found at

http://en.wikipedia.org/wiki/Logarithm

On 2009-11-14 17:52:56 -0400, Luna Moon <lunamoonmoon@gmail.com> said:

> Hi all,
>
> I am helping my friend on this issue.
>
> We have some numerical computation which involves a lot of "log(g)*g"
> and we have to evaluate these for g on the interval [0, 1],
> numerically... The reason that 0 is included is that because the
> singularity at 0 is removable.

Just as there is a more common name for sin(x)/x as sinc(x) which
is fairly common in signal processing there may be a common name
for x*log(x) which shows up in entropy calculations. There may even
be known approximations for x*log(x) as there is for sin(x)/x. Google
is your friend as I do not know of a ready source.

Try defining a new function xln(x) = x * ln(x) in Maple and let it chew
on the result. At least it will be easier to know where the simplifiction
will occur.

> In Maple:
>
>> x := 0; log(x)*x;
> %;
> 0
> Error, (in ln) numeric exception: division by zero
>
> In Matlab:
>
>
>>> x=0; log(x)*x
>
> ans =
>
> NaN
>
>
> Our procedure is that we first derive very complicated symbolic
> expressions in Maple and then translate into Matlab to get numerical
> calculation.
>
> How do handle the above "NaN" in Matlab? This is just one intermediate
> step using "log(x)*x" (this expression is from Maple symbolic
> results).
>
>
> Of course we could take a limit at x=0 and handle that case in Matlab
> separately, let say when |x| < 1e-5, we switch to this limiting
> result; and when |x|>=1e-5, we still use the numerical approach.
>
>
> But this creates a discontinuity at |x|=1e-5 in numerical results in
> Matlab.
>
>
> Could anybody shed some lights on us about how to handle such a
> situation?
>
>
> Thanks a lot!

BTW, in maple, for such things, I tend to use the piecewise function

xlog := x->piecewise(x=0, 0, x*log(x));

which should return 0 when x is 0 else x*log(x).

I suggest in MATLAB:
x=0+eps;
ie,using eps to avoid the situation where some is divided by zero

On Nov 15, 4:19 am, Gordon Sande <g.sa...@worldnet.att.net> wrote:
> On 2009-11-14 17:52:56 -0400, Luna Moon <lunamoonm...@gmail.com> said:
>
> > Hi all,
>
> > I am helping my friend on this issue.
>
> > We have some numerical computation which involves a lot of "log(g)*g"
> > and we have to evaluate these for g on the interval [0, 1],
> > numerically... The reason that 0 is included is that because the
> > singularity at 0 is removable.
>
> Just as there is a more common name for sin(x)/x as sinc(x) which
> is fairly common in signal processing there may be a common name
> for x*log(x) which shows up in entropy calculations. There may even
> be known approximations for x*log(x) as there is for sin(x)/x. Google
> is your friend as I do not know of a ready source.
>
> Try defining a new function xln(x) = x * ln(x) in Maple and let it chew
> on the result. At least it will be easier to know where the simplifiction
> will occur.
>


There is a big difference between sin(x)/x and x*log(x); the first is
analytic
(removable singularity at 0) but the second is not. `xln` is
continuous in
[0,oo) but is not differentiable at 0, so it will generate many
problems.

Mate

Hi,

On Nov 15, 2:19 am, Gordon Sande <g.sa...@worldnet.att.net> wrote:

> Just as there is a more common name for sin(x)/x as sinc(x) which
> is fairly common in signal processing there may be a common name
> for x*log(x) which shows up in entropy calculations. There may even
> be known approximations for x*log(x) as there is for sin(x)/x. Google
> is your friend as I do not know of a ready source.

A good idea.

Another idea, not necessarily good and its applicability depends a lot
on the remaining part of your problem.
Do a variable substitution exp(y)=x, hence x*ln(x) becomes exp(y)*ln
(exp(y)) which is exp(y)*y, defined on the entire real line. The
substitution is only really useful in case you can do it on all of
your x variables.

Sebastian

Gordon Sande wrote:
 > On 2009-11-14 17:52:56 -0400, Luna Moon <lunamoonmoon@gmail.com> said:
 >

 >> We have some numerical computation which involves a lot of "log(g)*g"
 >> and we have to evaluate these for g on the interval [0, 1],
 >> numerically... The reason that 0 is included is that because the
 >> singularity at 0 is removable.
 >
 > Just as there is a more common name for sin(x)/x as sinc(x) which
 > is fairly common in signal processing there may be a common name
 > for x*log(x) which shows up in entropy calculations. There may even
 > be known approximations for x*log(x) as there is for sin(x)/x. Google
 > is your friend as I do not know of a ready source.

xlog(x) has somewhat nastier behaviour at the origin than sinc(x).

QAGS in GNU SL will probably handle this integral directly. They give a
similar tricky one on their examples page. Never used their code in
anger myself.

http://www.gnu.org/software/gsl/manual/html_node/Numerical-integration-examples.html

Regards,
Martin Brown