3
$\begingroup$

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
Improve this question
$\endgroup$
9
  • 6
    $\begingroup$ 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 $\endgroup$
    – Bill
    Commented May 16, 2020 at 2:33
  • 5
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented May 16, 2020 at 2:35
  • 3
    $\begingroup$ @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. $\endgroup$
    – MarcoB
    Commented May 16, 2020 at 4:23
  • 2
    $\begingroup$ 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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 16, 2020 at 5:02
  • 2
    $\begingroup$ 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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 18, 2020 at 5:11

0

Reset to default