0
$\begingroup$

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?

Improve this question
$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$
    – J. M.'s missing motivation ♦
    Oct 15, 2018 at 16:36

2 Answers 2

Reset to default
4
$\begingroup$

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.2`20;
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.

Improve this answer
$\endgroup$
1
$\begingroup$

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 *)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.