MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Zero testing

  • To: mathgroup at smc.vnet.net
  • Subject: [mg53133] Re: [mg53106] Zero testing
  • From: Andrzej Kozlowski <akoz at mimuw.edu.pl>
  • Date: Fri, 24 Dec 2004 05:59:39 -0500 (EST)
  • References: <200412231259.HAA21206@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On 23 Dec 2004, at 21:59, Maxim wrote:

> By 'zero testing' I mean the problem of deciding whether a given 
> numerical
> quantity is exactly equal to zero. Often Mathematica doesn't make
> necessary checks at all:
>
> In[1]:=
> a = Pi/4; b = I*Log[-(-1)^(3/4)];
> Limit[(a - b)/x, x -> 0]
> Integrate[1/((x - a)*(x - b)), x]
> Integrate[1/(x - b), {x, 0, 1}]
>
> Out[2]=
> DirectedInfinity[Pi - (4*I)*Log[-(-1)^(3/4)]]
>
> Out[3]=
> (-4*((-2*I)*ArcTan[Log[-(-1)^(3/4)]/x] - 2*Log[Pi - 4*x] + Log[16*(x^2
> + Log[-(-1)^(3/4)]^2)]))/(2*Pi - (8*I)*Log[-(-1)^(3/4)])
>
> Out[4]=
> ((-I)*Pi - (2*I)*ArcTan[Log[-(-1)^(3/4)]] - Log[Log[-(-1)^(3/4)]^2]
> + Log[1 + Log[-(-1)^(3/4)]^2])/2
>
> The constants a and b are equal to each other. We can see that Out[2] 
> is
> incorrect because it is equal to ComplexInfinity and the limit of 0/x 
> is
> 0; Out[3] is also equal to ComplexInfinity (the denominator is zero) 
> and
> therefore wrong as well; finally, Out[4] is incorrect because the 
> integral
> doesn't converge.

This problem amounts to the difference between Simplify and 
FullSimplify. FullSimplify is needed to determine that a and b are 
equal, just Simplify cannot achieve this:

In[1]:=
a = Pi/4; b = I*Log[-(-1)^(3/4)];

In[2]:=
Simplify[a - b]

Out[2]=
(1/4)*(Pi - 4*I*Log[-(-1)^(3/4)])

In[3]:=
FullSimplify[a - b]

Out[3]=
0

Functions like integrate use Simplify automatically but not 
FullSimplify. Why not? The answer you provided yourself (well, almost) 
in an earlier posting: FullSimplify generally takes far too long to be 
used automatically. Even using Simplify leads to problems, as you have 
noted. Are you then arguing that FullSimplify ought to be used 
automatically? I know it ould b enice if one were able to eat the cake 
and still have it ...

If you apply FullSimplify to the argument in the above examples you 
will get the answer you stated as correct. In all cases, that is, 
except in one when it is you and not Mathematica that blundered. The 
case in question is:


> Integrate[1/((x - a)*(x - b)), x]


Are you seriously claiming that

Integrate[1/((x - Pi/4)^2), x]

is ComplexInfinity? At least Mathematica knows better than that:


Integrate[FullSimplify[1/((x - a)*(x - b))], x]

4/(Pi - 4*x)

>
> A variation of the same problem is when we're numerically evaluating a
> discontinuous function for an argument close to a discontinuity point. 
> In
> those cases Mathematica does perform necessary checks but handles more
> complicated cases incorrectly, especially when the branch cuts are
> concerned:
>
> In[5]:=
> Sign[Im[Sqrt[-1 - 10^-18*I]]]
>
> Out[5]=
> 1


The problem here is that Mathematica has great difficulty in 
determining Im[Sqrt[-1 - 10^-18*I] is -1. Algebraically. You ned ot use


Sign[FullSimplify[ComplexExpand[Im[Sqrt[-1-10^-18*I]]]]]

-1


This I agree to be a problem that ought to be fixed. It seems to me 
that numerical methods should be used since:


Sign[N[Im[Sqrt[-1-10^-18*I]]]]

-1

works easily and reliably. On the other hand it could be argued that 
this sort of thing should be approached numerically anyway.
>
> In[6]:=
> ((1 + I*(-167594143/78256779 + Pi))^4)^(1/2)
>
> Out[6]=
> (1 + I*(-167594143/78256779 + Pi))^2
>
> In fact In[5] is equal to -1 and In[6] is equal to -(1
> + I*(-167594143/78256779 + Pi))^2. In the last case we do not even 
> apply
> any functions, the error is in the automatic simplifications. Sometimes
> setting a high value for $MinPrecision helps, but here it doesn't have 
> any
> effect either.
>

This one is indeed a problem. Even this

FullSimplify[ComplexExpand[
     (1 + I*(-167594143/
          78256779 + Pi))^4]]^
   (1/2)


-(((167594143 + 78256779*I) -
      78256779*Pi)^2/
    6124123459454841)

gives the wrong answer ( the sign is wrong). To get the right one you 
need to use high precision arithmetic and Unevaluated:


N[Unevaluated[
    ((1 + I*(-167594143/
          78256779 + Pi))^4)^
     (1/2)], 20]


1.5635873239815087691232`4.043\
607143112669*^-16 -
   2.00000000000000015635873239\
81508708003`20.15051499783199*I

The only way that I have found of getting the correct exact answer is:


FullSimplify[ComplexExpand[((1 + I*(-u + Pi))^4)^
      (1/2)] /. u -> 167594143/78256779]

((167594143 + 78256779*I) - 78256779*Pi)^2/
   6124123459454841

In[54]:=
N[%, 20]

Out[54]=
1.5635873239815087691362`4.043607143112669*^-16 -
   2.00000000000000015635873239815087005589`20.15051499783\
1987*I

It is curious that this method works while the one using Unevaluated 
does not:


FullSimplify[ComplexExpand[((1 + I*(-u + Pi))^4)^
      (1/2)] /. u -> 167594143/78256779]
FullSimplify[ComplexExpand[Unevaluated[
     ((1 + I*(-167594143/78256779 + Pi))^4)^(1/2)]]]


((167594143 + 78256779*I) - 78256779*Pi)^2/
   6124123459454841


-(((167594143 + 78256779*I) - 78256779*Pi)^2/
    6124123459454841)



While I think this Mathematica's dealing wiht this situation could be 
improved (particularly the last example) I find it  difficult to see 
how such problems could be avoided in general. The problem seems to be 
very subtle. On the one hand, one can't expect the Mathematica 
evaluator to automatically apply powerful transformation s such as 
ComplexExpand or FullSimplify. On the other hand, certain 
transformations have to be automatically performed by the evaluator 
otherwise Mathematica would become unbearably tedious to use (one could 
leave every expression unevaluated until Simplify or FullSimplify or 
possibly a weaker function were applied to it but it would still 
require the user to judge when to apply which).

As far as I can tell all this means that the user has to remain 
vigilant at all times and always think carefully about the answers 
returned by the program.


Andrzej Kozlowski
Chiba, Japan
http://www.akikoz.net/~andrzej/
http://www.mimuw.edu.pl/~akoz/


  • References:
  • Prev by Date: Re: Need your help to solve a PDE.
  • Next by Date: Re: Re: Re: Mathematica language issues
  • Previous by thread: Zero testing
  • Next by thread: Re: Zero testing