Using Mathematica 10.1, depending on whether I write a numerical constant as 0.6x10^-15, or 6x10^-16, I get very different answers from the following code to solve a cubic equation. The first form leads to a presumed, 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?
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 *)
Presumed 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