# equal expressions give different numerical result

I have two expressions, and they are really equal:

aa = -32*pc1T1*(CAp*(CB*(-(CBp*m*SA) + CA*(E1 + E2)*SBp) +
SB*(-(CA*CBp*(E1 + E2)) + m*SA*SBp)) +
SAp*(-((E1 + E2)*SA*(CBp*SB - CB*SBp)) +
CA*m*(CB*CBp - SB*SBp)));
bb = -32 CA CAp CB E1 pc1T1 SBp -
32 CA CAp CB E2 pc1T1 SBp + 32 CA CAp CBp E1 pc1T1 SB +
32 CA CAp CBp E2 pc1T1 SB - 32 CA CB CBp m pc1T1 SAp +
32 CA m pc1T1 SAp SB SBp + 32 CAp CB CBp m pc1T1 SA -
32 CAp m pc1T1 SA SB SBp - 32 CB E1 pc1T1 SA SAp SBp -
32 CB E2 pc1T1 SA SAp SBp + 32 CBp E1 pc1T1 SA SAp SB +
32 CBp E2 pc1T1 SA SAp SB;


To see they are equal:

aa - bb // Simplify


This gives:

0.


However, if I evaluate these two expressions in numerical values, they produce very different result:

m=0.2;
aa /. {CA -> 5.050247469193957, CB -> 5.050247469193957,
CAp -> 5.050247469193957, CBp -> 5.050247469193957,
SA -> 4.950252468319164, SB -> 4.950252468319164,
SAp -> 4.950252468319164, SBp -> 4.950252468319164, pc1T1 -> 3,
E1 -> 5.834380858325929, E2 -> 5.834380858325929}
bb /. {CA -> 5.050247469193957, CB -> 5.050247469193957,
CAp -> 5.050247469193957, CBp -> 5.050247469193957,
SA -> 4.950252468319164, SB -> 4.950252468319164,
SAp -> 4.950252468319164, SBp -> 4.950252468319164, pc1T1 -> 3,
E1 -> 5.834380858325929, E2 -> 5.834380858325929}


Result is like the following and stable:

3.41061*10^-13
-2.32831*10^-10


One is plus sign and the other is minus sign. This is really strange. Does anyone know how to fix this?

• Just because two expressions are algebraically equivalent does not mean that they will also be numerically equivalent, when inexact arithmetic is involved. Part of the business of numerical analysis is the design and modification of formulae to make them numerically stable. – J. M.'s technical difficulties Oct 15 '18 at 16:36

Both aa and bb vanish for those values of the parameters. The discrepancy is just numerical error. Use a higher precision to get the correct answer:

m = 0.220;
aa /. SetPrecision[{CA -> 5.050247469193957,
CB -> 5.050247469193957, CAp -> 5.050247469193957,
CBp -> 5.050247469193957, SA -> 4.950252468319164,
SB -> 4.950252468319164, SAp -> 4.950252468319164,
SBp -> 4.950252468319164, pc1T1 -> 3, E1 -> 5.834380858325929,
E2 -> 5.834380858325929}, 20]
bb /. SetPrecision[{CA -> 5.050247469193957,
CB -> 5.050247469193957, CAp -> 5.050247469193957,
CBp -> 5.050247469193957, SA -> 4.950252468319164,
SB -> 4.950252468319164, SAp -> 4.950252468319164,
SBp -> 4.950252468319164, pc1T1 -> 3, E1 -> 5.834380858325929,
E2 -> 5.834380858325929}, 20]

(* 0.*10^-14 *)
(* 0.*10^-13 *)


In fact, you can establish the vanishing of these expressions symbolically:

Clear[m]
FullSimplify[
aa /. {CA -> CB, CAp -> CB, CBp -> CB} /. {SA -> SB, SAp -> SB, SBp -> SB}
]
(* 0 *)


and a similar calculation for bb.

It is a precision issue

aa = -32*pc1T1*(CAp*(CB*(-(CBp*m*SA) + CA*(E1 + E2)*SBp) +
SB*(-(CA*CBp*(E1 + E2)) + m*SA*SBp)) +
SAp*(-((E1 + E2)*SA*(CBp*SB - CB*SBp)) + CA*m*(CB*CBp - SB*SBp)));

bb = -32 CA CAp CB E1 pc1T1 SBp - 32 CA CAp CB E2 pc1T1 SBp +
32 CA CAp CBp E1 pc1T1 SB + 32 CA CAp CBp E2 pc1T1 SB -
32 CA CB CBp m pc1T1 SAp + 32 CA m pc1T1 SAp SB SBp +
32 CAp CB CBp m pc1T1 SA - 32 CAp m pc1T1 SA SB SBp -
32 CB E1 pc1T1 SA SAp SBp - 32 CB E2 pc1T1 SA SAp SBp +
32 CBp E1 pc1T1 SA SAp SB + 32 CBp E2 pc1T1 SA SAp SB;

aa == bb // Simplify

(* True *)

param = {m -> 0.2, CA -> 5.050247469193957, CB -> 5.050247469193957,
CAp -> 5.050247469193957, CBp -> 5.050247469193957,
SA -> 4.950252468319164, SB -> 4.950252468319164,
SAp -> 4.950252468319164, SBp -> 4.950252468319164, pc1T1 -> 3,
E1 -> 5.834380858325929, E2 -> 5.834380858325929};

{aa, bb} /. param

(* {3.41061*10^-13, -2.32831*10^-10} *)

Subtract @@ %

2.33172*10^-10


However, if expressions are simplified prior to applying numeric values

Simplify[{aa, bb}] /. param

(* {8.73115*10^-11, 8.73115*10^-11} *)

Subtract @@ %

(* 0. *)


or

FullSimplify[{aa, bb}] /. param

(* {-8.55209*10^-11, -8.55209*10^-11} *)

Subtract @@ %

(* 0. *)


Alternatively, if all approximate numbers are converted to exact values, both expressions are seen to be identically zero

{aa, bb} /. Rationalize[param, 0]

(* {0, 0} *)


Or using RootApproximant

{aa, bb} /. (param /. x_Real :> RootApproximant[x]) // RootReduce

(* {0, 0} *)


EDIT: Using arbitrary precision

{aa, bb} /. SetPrecision[param, 15]

(* {0.*10^-9, 0.*10^-8} *)

Subtract @@ %

(* 0.*10^-8 *)
`