3

Using Mathematica 10.1

Depending on whether I write a numerical constant as 0.6 * 10^-15, or 6 * 10^-16, I get very different answers from the following code to solve a cubic equation. The first form leads to a nonphysical, imaginary root, and the second to a physical real root. When printed, the numbers are not different.

I know now to avoid the first form, but why is the answer wrong?

P0 = 2*10^-6; E0 = 1*10^-5; Z0 = 3*10^-6;
KdP = 8*10^-16; (* assumed protein Kd *) 
KdE = 0.6*10^-15;  
N[KdE]
a := KdP + KdE + P0 + E0 - Z0;
b := KdE*(P0 - Z0) + KdP*(E0 - Z0) + KdP*KdE;
c := -KdP*KdE*Z0;
theta := ArcCos[(-2*a^3 + 9*a*b - 27*c)/(2 Sqrt[(a^2 - 3*b)^3])];
Z := -a/3 + (2/3)*Cos[theta/3]*Sqrt[a^2 - 3 b];
PZ := P0*Z/(Z + KdP);
EZ := E0*Z/(Z + KdE);
N[Z + EZ + PZ, 3] (* Zn conserved? *)
N[PZ/P0, 6] (*fraction Zn bound to P *)
N[EZ/E0, 6] (*same for EDTA *)

Incorrect result from above:

6.*10^-16 
0.0000119979 + 2.08196*10^-7 I 
0.999686 + 0.0219132 I 
0.999855 + 0.016437 I

Correct result if I write KdE = 6 * 10^-16

6.*10^-16
3.00*10^-6
0.207304
0.258539
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
jremington
  • 31
  • 2
  • 6
    0.6x10^-15 has a decimal point in it and is thus an approximate floating point number. 6x10^-16 has no decimal point and thus is an exact rational number. Every other constant in your code is an exact number and thus all the calculations will be exact if you don't have a decimal in KdE, down to the point where you use N to display a decimal version of the result, but those decimals are not further used in calculations. Often Mathematica chooses completely different algorithms when dealing with approximate numbers than it does when dealing with exact numbers. That likely explains your code – Bill May 16 '20 at 02:33
  • 5
    To further confirm that it has to do with exact vs. Approximate numbers, notice that, if you use 6. *10^-16 instead (note the period!), which is a machine precision version, you get the same unphysical result. – MarcoB May 16 '20 at 02:35
  • Thanks, but that is a bit discouraging! – jremington May 16 '20 at 02:49
  • 3
    @jremington Why is it discouraging? You should always be aware of whether you are working with machine-precision numbers, or using the arbitrary precision system. There can be VERY large differences in speed between the two as well. – MarcoB May 16 '20 at 04:23
  • 2
    Also, you seem to be using Cardano's formula manually. Why not use Solve[] or NSolve[] instead? But if you insist on explicitly solving a cubic, there are numerical pitfalls to be aware of; see this. – J. M.'s missing motivation May 16 '20 at 05:02
  • You would also get the expected results using arbitrary-precision, KdE = 0.6\16*10^-15;` – Bob Hanlon May 16 '20 at 05:05
  • @MarcoB: it is discouraging, because as an educator, I think about students encountering this for the first time. The user is not informed of the switch, or the difference in approaches, so it is a trap set for the unsuspecting. It is not helpful to state that one should always understand all of the details the tool one is using, but thank you anyway. – jremington May 17 '20 at 21:37
  • 2
    It is common that a lot of the classical formulae for computation (even the simple quadratic formula!!) are symbolically good, but numerically unreliable. Someone using these formulae with pen and paper can easily adjust, but the computer does not have faculties like that. Better for your students to know that while they are still students. – J. M.'s missing motivation May 18 '20 at 05:11
  • 1
    Someone once told me: what you want to do is one thing; what you tell the computer to do is another; what the computer does is yet something else. When using a computer to do math, one needs to be aware that the computer uses a different math than the user expects. E.g., look here. Moreover, the computer doesn't even multiply numbers in the way we're taught in school - it uses a whole different, more efficient algorithm. And, in spirit of J.M.'s comment, as an educator you have the chance to educate the students how does a computer work. – corey979 May 18 '20 at 08:09

0 Answers0