How to cancel floating point factors?

Question

I am currently having problems with some floating points.

I have a function, which gives as an intermediate result (for example)

-(10000. - 10000. a) E^(-16.4157 kp^2 x^0.277)/((-1. + a) a x^0.252525 xj^4.47047)

It gets passed through into another function, which will fill in the value a = 1. Then Mathematica reacts with a 1/Sqrt[0] error. However, as you can see, the factor (10000. - 10000. a) should cancel with the factor (-1.+ a) in the denominator to give 10000, the result is thus well-defined.

Do you have any idea how to let Mathematica cancel these factors out? I have tried Simplify and FullSimplify with Assumptions -> a!=1 but it doesn't work. I cannot change much to the function itself, as it is an intermediate result (and it is just an example for this post; other intermediate results also occur, sometimes with the same problem, sometimes they work fine).

Floating point numbers are different than rational numbers. Since floating point numbers don't really follow the rules of algebra, you'll find that Simplify and the like will not do much at all wit them.

Run the calculations with increased precision using SetPrecision. (still essentially the same problem)

I think there are two options:

Use the function Rationalize to turn the floating point numbers into rational numbers and then use FullSimplify. (I'm not a big fan of this because it ignores the fact that the numbers of a limited precision) — Searke, Mar 01 '12 at 19:47
@David That is right. Addition/multiplication is not even associative. — Szabolcs, Mar 01 '12 at 19:57
@David for example, (1 + $MaxMachineNumber) == $MaxMachineNumber evaluates to True — acl, Mar 01 '12 at 20:02
@David try 1.*10^16 + 1. - 1.*10^16 which returns different answers depending on where 1. is, however -10.^16 + (10.^16 + 1.) and -10.^16 + (1. + 10.^16) return the same incorrect answer. This indicates that addition is not associative, but commutative. Also, given g[a_, b_, c_] := a^2 - 2 a b + b^2 - c^2, g[10.^10, 10.^10, 10.] == 0 implying that mma is doing speculative processing as the first three terms simplify to $(a - b)^2$ and should cancel. — rcollyer, Mar 01 '12 at 20:43
Even more odd, 1. + (1.1 $MachineEpsilon) === 1. evaluates to True. In any case, point taken, thanks everyone. — David, Mar 01 '12 at 20:47
A great example is that if you use Floating point numbers in DSolve , you will find that its ability to solve differential equations becomes greatly reduced because of all this. They are difficult to reason about algebraically. — Searke, Mar 01 '12 at 21:28
Wow, a lot new stuff I didn't realize existed! I did the trick by trying to postpone in general my numerics (keeping the calculation symbolic until the last step), and adding Rationalize and Cancel for those situations where I cannot keep the calculation symbolic. I couldn't implement Limit, but I'm happy that I know it now. Thanks a lot everyone! — freddieknets, Mar 01 '12 at 23:32

score 9 · Accepted Answer · answered Mar 01 '12 at 19:53

9

The problem is that symbolic computation doesn't mix well with floating point. If you use Rationalize first to convert to an exact representation, Simplify (or the faster Cancel) will be able to do it's job.

answered Mar 01 '12 at 19:53

Szabolcs

234,956
30
623
1,263

David · Answer 2 · 2012-03-01T19:53:16.127

5

Use Limit, Mathematica's implementation of $\lim$:

$\displaystyle\lim_{x\to a}f(x) = $ Limit[f[x], x -> a]

Applied to your equation, this results in

Limit[
    -(((10000. - 10000. a) E^(-16.4157 kp^2 x^0.277)/((-1. + a) a x^0.252525 xj^4.47047))),
    a -> 1
]

(10000. E^(-16.4157 kp^2 x^0.277))/(x^0.252525 xj^4.47047)

The error you are talking about is not caused by floating point arithmetic, it's how Mathematica (or any other language for that matter) handles computations in general. If you substitute $a=1$ then you've got a division by zero, and once this occurs you're out of the domain of, well, mathematics; therefore, you can't cancel the two factors anymore. (Actually, your function is not defined at $a=1$, it has a removable singularity there.)

edited Mar 01 '12 at 19:53

answered Mar 01 '12 at 19:45

David

14,911
6
51
81

The function is well defined at a = 1. The singularity is artificial, and has more to do with the order of evaluation than a being a property of the function, unlike in $\sin(x)/x$. Since mma evaluates the terms in isolation from one another, i.e. $(a - 1)^{-1}$ is evaluated separately from $1/a$, etc., this is likely to occur. – rcollyer Mar 01 '12 at 20:08
1

$\frac xx$ stands for the product of x and its multiplicative inverse, which does not exist for $x=0$. It's not well-defined or even defined at that point, although the limit exists of course. – David Mar 01 '12 at 20:10
On the surface, I'd disagree with you, as it is plain that $\frac xx$ simplifies to $1$ prior to the limit being taken. But, upon further review, it appears that removable singularities may in fact correspond to how we evaluate things, as opposed to being intrinsic. After all, $\frac\sin{x} = 1 - x^2/3! + x^4/5!\ \cdots$, implying that the singularity has more to due with the evaluation process than actually being present. – rcollyer Mar 01 '12 at 20:34
It's the same thing for $\sin(x)/x$, not defined at $x=0$. It's of course very convenient to always assume that the removable singularity is removed (by defining the additional point using the limit). Mathematically what you do when canceling fractions is choosing another element of the equivalence class of two numbers, but if the denominator is zero you don't have a rational number in the first place, meaning no arithmetic is defined. (Also see Wikipedia, Constructing the rationals) – David Mar 01 '12 at 20:41
2

This, of course, reminds me why I hated my analysis course. – rcollyer Mar 01 '12 at 20:44

score 5 · Answer 3 · answered Mar 01 '12 at 21:42

Well...

[Warning: This act is performed by trained professionals. Do not try it at home.]

In[24]:= expr = -(10000. - 10000. a) E^(-16.4157 kp^2 x^0.277)/
  ((-1. + a) a x^0.252525 xj^4.47047);

In[25]:= Cancel[expr]
Out[25]//InputForm= 10000./(a*E^(16.4157*kp^2*x^0.277)*x^0.252525*xj^4.47047)

Okay, I cheated. I'm using the development version of Mathematica.

I will add some emphasis to remarks others have made. Doing computer algebra with approximate numbers is almost never a safe bet. Even to the extent that it might be supported, things can go wrong. And do. On a remarkably frequent basis.

How to cancel floating point factors?

3 Answers3

Linked

Related