Numerical voodoo in mathematica? [closed]

So, I was doing some calculations in Mathematica, and came across the following paradox! (I may be exaggerating, but...) I wanted to solve the following equation:

$3.33333*10^{-15} - 3.33333*10^{-15} \exp(-t/5000000) = 10^{-122}$.

I used NSolve to solve this, and I got $t = 0.009951$.

To check the solution, when I substituted this back into the expression, I got: $1.59926*10^{-32}$ which is obviously not correct!

Anyone have an idea of what is going on?

• NSolve gives a machine-precision number while regular calculations give an arbitrary-precision number. This causes a rounding error and a fallacious answer especially when the expected answer is very small. – JungHwan Min Jul 27 '16 at 21:17
• Okay. Thanks for this. So, which is the fallacious answer? Is it the result from solving the equation of t = 0.009951, or the substituted check? – Thomas Moore Jul 27 '16 at 21:24
• I vote to close this question as unclear. It could be made more clear by including the code that you used to solve the equation. When I input the obvious : Solve[3.33333*10^-15 - 3.33333*10^-15*Exp[-t/5000000]==10^-122,t]  I get 0.0. – Mark McClure Jul 28 '16 at 11:56

To get exact output, use exact input (and please input your question using plain text Mathematica code that one can easily copy).

Clear[t];
x = Rationalize[3.33333];
eq = x*10^(-15) - x*Exp[-t/5000000]*10^(-15) == 10^(-122);
sol = t /. First@Solve[eq, t];
sol = sol /. C[1] -> 0


eq /. t -> sol


t should be positive. Since

$MaxExtraPrecision = 1000; N[Log[333333000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000000/ 3333329999999999999999999999999999999999999999999999999999999999999999999999\ 99999999999999999999999999999999], 1000]  and this is little over zero, and it is being multiplied by 5000000 since t was so the final t should be little over exact zero. N[5000000 Log[ 333333000000000000000000000000000000000000000000000000000000000000000000000\ 000000000000000000000000000000000/ 333332999999999999999999999999999999999999999999999999999999999999999999999\ 999999999999999999999999999999999], 1000]  • Hi. Thanks for this. So, if I take C[1] = 0 in this calculation, t works out to be 0.009951, which is similar to what I got with NSolve. So, is t = 0.009951 the correct answer, or is it wrong? – Thomas Moore Jul 27 '16 at 21:23 • @ThomasMoore You are mistaken, taking C[1]=0 Yields -3.45388*10^8 – Feyre Jul 27 '16 at 21:25 • @Feyre I think there is a typo in the Eq that Nasser posted. There should be a 10^(-15) in the first term, not 10^(15). – Thomas Moore Jul 27 '16 at 21:26 • @ThomasMoore will correct the typo now. It does not affect the process. (You see, if you pasted your code in plain text, I would not had to type it all over and make mistake :) – Nasser Jul 27 '16 at 21:28 • Okay, will remember that for next time. – Thomas Moore Jul 27 '16 at 21:28 The problems at hand comes from using inadequate numerical approximation. The equation behind is quite simple and it can be solved exactly, we put 10/3 instead of 3.33333 (if usual mathematical meaning is assumed use 333333/100000 instead of 10/3) : t0 = t /. Solve[10/3 (1 - Exp[-(t/5000000)]) 10^-15 == 10^-122, t, Reals]//First  5000000 (108 Log[2] + 108 Log[5] - Log[23] - Log[83] - Log[137401213] - Log[3207480541] - Log[92762470571] - Log[128134821352182563257005775275977598174025940104318328372829269\ 39321359468931])  Numerically it can be well estimated with Block[{$MaxExtraPrecision = 200}, N[t0, 112]]

1.50000000000000000000000000000000000000000000000000000000000000000000\
0000000000000000000000000000000000000002250*10^-101


The symbolic result is the first thing we can rely on (in general, there are certain exceptions like bugs). If no exact approach is available one should carefully consider problems case by case, nevertheless one can be quite sure Mathematica provides appropriate tools to harness.

When substituting back, the computation that occurs is subtraction $10^{-15}$ (with 16 significant figures, hence 31 correct digits in the fractional part) from something very similarly. At the end you still have about 31 correct digits in the fractional part, which is why you are not able to obtain $10^{-122}$ whose first significant figure is the 122nd digit of the fractional part.

• Okay, I see. So, is there any way of confirming that t = 0.009951 is the correct answer? – Thomas Moore Jul 27 '16 at 21:25

"$3.33333 \times 10^{−15}$" is not a sharp number. It represents the interval $(3.333325 \times 10^{-15}, 3.333335 \times 10^{-15})$. Replacing your unsharp numbers with the appropriate intervals, ...

Reduce[
Interval[{3.333325*^-15, 3.333335*^-15}] -
Interval[{3.333325*^-15, 3.333335*^-15}] Exp[-t/5000000] ==
Interval[{0.5*10^-122, 1.5*10^-122}],
t, Reals]

(*  t == Interval[{-15.,15.}]  *)


So, given the (poor) precision of your inputs, the answer could justifiably anything between $-15$ and $15$. If you want a narrower interval, you must provide more precise, or exact numbers in your equation. For example,

Reduce[
Interval[{3.333333333325*^-15, 3.333333333335*^-15}] -
Interval[{3.333333333325*^-15, 3.333333333335*^-15}] Exp[-t/5000000] ==
Interval[{0.5*10^-122, 1.5*10^-122}],
t, Reals]

(*  t == Interval[{-0.0000150041, 0.0000150047}]  *)


The above assumes that the two uses of Interval[{3.333325*^-15, 3.333335*^-15}] are independent. If not, we could do more work, but we would get the same interval for this problem (because the two uses of this interval factor out to a single usage).

• Did you know that applying Interval[] directly to an inexact number automatically generates the corresponding interval? – J. M.'s discontentment Jul 28 '16 at 18:42
• @J.M. : I did not. I do now. Thanks. – Eric Towers Jul 29 '16 at 1:47

As others have said, the issue with NSolve is loss of precision. By the way, well done for at least back-substituting your answer to test it -- good idea!

You can get a correct answer with pencil and paper, and this is a good skill to practise. Let x = (t/5e6). Divide everything by 3.33333e-15 and rearrange: exp(-x) = 1 - 1e-107/3.33333. Now it should be clear that 0 < x << 1. So, a first-order approximation of exp suggests that x ~ 1e-107/3.33333, and this shows that a second-order approximation is unnecessary if you only want 100 or so correct digits. So t ~ 5e-101/3.33333. For 6-digit accuracy, we could regard 3.33333 ~ 10/3 and the answer is t ~ 1.50000e-101.

The pencil-and-paper method has avoided loss of precision because we reached an equation (1-x) = (1-small), and we cancelled out those two "large" 1s exactly. (OK, this assumes that the first 3.33333e-15 exactly equals the second. This assumption is not necessarily justified. The answer will be very sensitive to the precision of the inputs in this equation.)