3

I have the following expression:

     eqs3=-2.750863054772633*^7*Subscript[x, 18] + 197943.17573209677*Sqrt[1.03*^7 - 1.*Subscript[x, 19]] + 
 2.750863054772633*^7*Subscript[x, 19] + 
 (-95931.14816252001 + (1.2622519495068421*^9*
     (-0.0001568174846625767*Sqrt[1.03*^7 - 1.*Subscript[x, 18]] + 0.021793296146996537*
       Subscript[x, 18] + 0.000076*Sqrt[1.03*^7 + Subscript[x, 18]] - 
      0.021793296146996537*Subscript[x, 19]))/Sqrt[1.03*^7 + Subscript[x, 18]])*
  Sqrt[1.03*^7 + Subscript[x, 19]]

I'm going to apply several rules, as follows.

eqs3 /. 
 Thread[{Subscript[x, 19], Subscript[x, 18]} -> -8673600.700]

(*0*)

eqs3 /. 
 Thread[{Subscript[x, 18], Subscript[x, 19]} -> -8673600.700]

(*0*)

eqs3 /. Subscript[x, 18] -> -8673600.700 /. 
 Subscript[x, 19] -> -8673600.700

(*0*)

So far so good, if it wasn't that...

eqs3 /. Subscript[x, 19] -> -8673600.700 /. 
 Subscript[x, 18] -> -8673600.700

(*0.00842667*)

Why is so??? This is very much disappointing.

gwr
  • 13,452
  • 2
  • 47
  • 78
Mirko Aveta
  • 2,192
  • 11
  • 26
  • 1
    With the definition you give, eqs3[[8]] does not make much sense?! – gwr Feb 28 '17 at 15:19
  • Sorry I'll correct this. – Mirko Aveta Feb 28 '17 at 15:20
  • Exactly as @Mr.Wizard shows, I can also not reproduce your results. (I get the same as Mr.Wizard does) – gwr Feb 28 '17 at 15:24
  • I think this is a matter of precision on the coefficients of the expression. When copying the expression some of the coefficients must have been chopped. I'll try and sort this out by putting the bits of code in InputForm. @gwr – Mirko Aveta Feb 28 '17 at 15:28

2 Answers2

4

I believe you are operating in the range of numeric noise for the given computation, similar to Funny behaviour when plotting a polynomial of high degree and large coefficients. As astutely illustrated by Michael E2 in that Q&A setting false accuracy (i.e. Rationalize) doesn't produce a meaningful result. If we instead use arbitrary precision, which tracks precision, we see that your results are all equally bad:

eqsTracked = SetPrecision[eqs3, $MachinePrecision];
val = SetPrecision[-8673600.700, $MachinePrecision];

r1 = eqsTracked /. Subscript[x, 18 | 19] -> val
r2 = eqsTracked /. Subscript[x, 18] -> val /. Subscript[x, 19] -> val
r3 = eqsTracked /. Subscript[x, 19] -> val /. Subscript[x, 18] -> val

Precision /@ {r1, r2, r3}
0.*10^-1

0.*10^-1

0.*10^-1

{0, 0, 0}

So every result has zero digits of Precision.

Here is another way to look at this. If we convert every machine precision number to an explicit Interval equivalent of exact arithmetic, using Michael's formula, we find:

machineToInterval = 
  c_Real?MachineNumberQ :>
    Interval[ {1 - 2^-53, 1 + 2^-53} * SetPrecision[c, ∞] ];

val = -8673600.700 /. machineToInterval;

eqs3 /. machineToInterval /. Subscript[x, 18 | 19] -> val // N
Interval[{-0.21192, 0.21192}]

So at most we can say that your result is somewhere in the interval $[-0.21192, 0.21192]$

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • Please note that I have edited my question copying the expression in InputForm. – Mirko Aveta Feb 28 '17 at 15:30
  • Would you please explain further more the functioning og machineToInterval? I'm kind of puzzled. How does it manage to define the interval? What would be the mathematical meaning of this? Thanks for your answer, I have the feeling I'm learning something important. – Mirko Aveta Feb 28 '17 at 21:57
  • 1
    @MirkoAveta Assuming I didn't screw it up it should be just as Michael E2 describes in this answer. 2^-54 comes from the number of bits used in the binary representation of machine numbers. See another answer again by Michael E2. – Mr.Wizard Feb 28 '17 at 22:01
  • May ask a further question? I faced this problem when attempting to solve a very complicated system of 20 nonlinear equations. One of the variables after several substitutions turned out to be eqs3. I was wondering, if I could proceed with the intervals, as I'm doing now and what happens when I apply machineToInterval to more variables? How do these intervals act together? – Mirko Aveta Feb 28 '17 at 22:02
  • 1
    @MirkoAveta It may be best to post a new Question on that topic. Not all operators will work with Interval but the basic arithmetic ones do, as do a fair number of numeric functions, e.g. Through[{Min, Max, Floor, Sin, Abs, UnitStep}[Interval[{-3.5, 5}]]]. – Mr.Wizard Feb 28 '17 at 22:15
  • 1
    @MirkoAveta I just realized from Michael's analysis I should have used 2^-53 for the maximum error rather than supposed mean error. I'll update accordingly. – Mr.Wizard Feb 28 '17 at 22:46
2

This seems to be an issue with numerical precision. Note that the problem does not arise when you work with exact coefficients (e.g. using Rationalize):

eqs3r = Rationalize[ eqs3, 0];

With[ { c = Rationalize[-8673600.700] },
  SameQ[
    (eqs3r /. Subscript[x, 19] -> c /. Subscript[x, 18] -> c)
    ,
    (eqs3r /. Subscript[x, 18] -> c /. Subscript[x, 19] -> c)
  ]
]

True

gwr
  • 13,452
  • 2
  • 47
  • 78