MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

an example... Re: Why Mathematica does not issue a warning when the calculations

  • To: mathgroup at
  • Subject: [mg117803] an example... Re: Why Mathematica does not issue a warning when the calculations
  • From: Daniel Lichtblau <danl at>
  • Date: Thu, 31 Mar 2011 04:06:18 -0500 (EST)

Some background from recent notes.

>> In an example of this sort a program COULD, in principle, raise the 
>> precision of the input argument. That will only work if that input is 
>> exactly representable in the underlying arithmetic.

> Here is a fundamental disagreement.  When a human types 0.1 as a number 
> to be represented as a hardware float (binary) it is converted to binary 
> and has a certain bit representation, B.  That now is the number that 
> the computer has in its electronic memory.  B is close to, but not 
> exactly 0.1 decimal.  It is closer to 0.10000000000000000555  except 
> rounded to "machine precision". it is 0.1.
> People of a certain sophistication will realize that not all decimal 
> numbers have exact binary representations. Just the fraction 1/3 does 
> not have an exact decimal representation.  We must live with some 
> version of this radix stuff and hope we can explain it when necessary.
> [...]
>> For decimals, in programs that use base 10 behind the scenes, this can 
>> be done. If the input is decimal and the representation is binary and 
>> the decimal is not exactly representable (e.g. 1.3 in binary) then 
>> you'd be in trouble.

> Only with people who don't understand that 1.3 is not exactly 
> representable in binary, and if they really want that, they could use 
> 13/10.
> ....

>> Though it is necessary if we want Mathematica to emulate fixed 
>> precision arithmetic in any sensible way. I do not regard manufacture 
>> of digits or bits, that is, raising precision artificially, as 
>> "sensible", unless one has a priori a justification for doing so.

> I, on the other hand, see no reason to assume that some number in the 
> computer, is not a perfectly valid number, exactly what the user wants 
> to compute with, and is accurate to all digits and can be extended by 
> zeros if need be.
> [...]
> Executive summary of the paragraph above:
>  You couldn't possibly mean THAT number, which is so close to 1/10.  So 
> I will refuse to allow you to do arithmetic with it.
> [...]

It was here that I responded, in an earlier note, that I would provide 
an example. An example of what? Well, of something.

There is actually a story to this example. I'll get to that later.

Here is my example (not originally mine, but that's part of the 
story...). We are given some algebraic equations with approximate 
coefficients, and we wish to find solutions. So in Mathematica this is 
an NSolve problem.

Here is the setup.

x0 = .23; x1 = .119; x2 = .182; x3 = .282; x4 = .1049; Kc = \
34.4; Kr = 1.45; Kch = Kr;
eq1 = Kc (x0 - eh - ec)*(x1 - ec +
       ech)*(1 + 2*ec + 2*eh)^2 == (x2 + eh + 2*ec - ech)^2*(x3 -
        3*eh + 2*ec + ech)^2;
r21 = (x4 - x2 - 2*ec)/2 - 3*(x1 - ec - ech)/(2 Kr);
eq2 = eh ==
    r21 + Sqrt[r21^2 + Kr*(1 - Kr) (x4 - ech) (x2 + 2*ec - ech)];
r31 = (Kch (x2 + x4 + 2*ec) + x3 + x1 + 3*eh + ec)/(2*(Kch - 1));
r32 = (Kch (x2 + eh + 2*ec)*(x4 - eh) - (x3 + 3*eh + 2 ec)*(x1 -
         ec))/(Kch - 1);
eq3 = ech == r31 + Sqrt[r31^2 - r32];
vars = {ec, eh, ech};

Here are the actual equations. I use Mathematica InputForm because (1) 
it copy-pastes nicely from the Mathematica front end, and (2) it shows 
the input, as Mathematica sees it, in its full glory (the input's glory, 
not Mathematica's...).

eqns = {eq1, eq2, eq3}

Out[254]= {34.4*(0.119 - ec + ech)*(0.23 - ec -
      eh)*(1 + 2*ec + 2*eh)^2 ==
      (0.282 + 2*ec + ech - 3*eh)^2*(0.182 + 2*ec - ech + eh)^2,
    eh == (1/2)*(-0.0771 - 2*ec) +
        Sqrt[((1/2)*(-0.0771 - 2*ec) -
         1.0344827586206897*(0.119 - ec - ech))^2 -
            0.6525*(0.1049 - ech)*(0.182 + 2*ec - ech)] -
    1.0344827586206897*(0.119 - ec - ech),
    ech == 1.1111111111111112*(0.40099999999999997 + ec +
       1.45*(0.2869 + 2*ec) + 3*eh) +
     1.234567901234568*(0.40099999999999997 + ec +
          1.45*(0.2869 + 2*ec) + 3*eh)^2 -

      2.2222222222222223*(1.45*(0.1049 - eh)*(0.182 + 2*ec + eh) -
                 (0.119 - ec)*(0.282 + 2*ec + 3*eh))]}

Notice that some of these numbers do not appear to be in any way "nice". 
That is, they are long decimals with no repeating parts. One of them 
comes about simply by taking a reciprocal of a "nice" decimal.

In[269]:= (2 Kr) // InputForm

In[268]:= 1/(2 Kr) // InputForm

In a sense this actually is a nice number. Specifically, Rationalize 
will recognize it as a rational with a small denominator. But I'm not 
sure it is fair to use Rationalize as a measure of "niceness" in the 
context of this discussion. For our purposes, the sensible measure is 
hhis. Do we know what we have in such a way that will allow us to raise 
precision by appending zeroes without messing it up?

Obviously other decimal numbers are "nice" even if long. For example, 
2.2222222222222223 is clearly an approximation, with error in the last 
place (from arithmetic, which, as I keep saying, happens...), of 

Now we use NSolve, and get a result. I'll show it again in InputForm.

In[255]:= Timing[NSolve[{eq1, eq2, eq3}, {ec, eh, ech}]]

Out[255]= {0.102983, {{ec -> -15.1618, eh -> 4.85832,
    ech -> -10.528}, {ec -> -0.728199 - 0.00417776 I,
    eh -> 0.614256 - 0.13956 I,
    ech -> 0.287064 - 1.04778 I}, {ec -> -0.728199 + 0.00417776 I,
    eh -> 0.614256 + 0.13956 I,
    ech -> 0.287064 + 1.04778 I}, {ec -> -0.504525, eh -> 0.110101,
    ech -> -0.39278}}}

Out[255]= {0.1029829999999663, {{ec -> -15.161779300885101,
    eh -> 4.858321317991649,
        ech -> -10.527968401345873}, {ec -> -0.7281987191080787 -
        eh -> 0.6142563132407307 - 0.1395597138750204*I,
        ech -> 0.287064092184311 - 1.0477819258510883*I},
      {ec -> -0.7281987191080787 + 0.004177758076267676*I,
        eh -> 0.6142563132407307 + 0.1395597138750204*I,
        ech ->
     0.287064092184311 +
      1.0477819258510883*I}, {ec -> -0.5045254112449927,
        eh -> 0.11010134840801741, ech -> -0.3927799738963131}}}

It looks to my view that none of the numbers are "nice" decimals. I 
realize neither one of us has stated that they should be. Here is what I 
am getting at. I do not see how a view of arithmetic that involves 
throwing in binary or decimal zeros to increase precision will do the 
least bit of good for working with a problem like this, in terms of 
improving the quality of the result. (Mathematica will raise the 
precision in this manner, but for a different reason, noted below.)

I'm not sure this contradicts your view of how numerics should be 
handled. But it certainly is not a case where your arithmetic model will 
be of much use. Moreover, for the implementation we have, it would not 
work at all. The reason there is a technicality. We use Groebner bases 
as a first serious step (after basic preprocessing, that is), and this 
cannot be done with fixed precision. Also it requires a modest initial 
precision, which is the reason we do in fact raise it from what is in 
the input. But I make no pretense that that will improve the quality of 
the result, only the likelihood that we actually get a result.

You might prefer to rationalize and use an exact solver. Fine, but slow. 
It has not completed after several minutes thus far.

A variation on NSolve will make an exact computation run in around three 

In[272]:= Timing[
  exsols = NSolve[Rationalize[eqns], vars,
     WorkingPrecision -> Infinity];]
Out[272]= {168.239, Null}

But you won't want to look at the result. Also N at machine precision 
will give a bad outcome (plug into equations, residuals of one will be 

End of aside]

Or you can rationalize and use NSolve. But this will still numericize, 
and that will take us into computations that will require significance 
arithmetic. (Again, that may be an implementation detail. But it is a 
serious one.)

I guess this is a bit open-ended. But here is my question. Does your 
view of numerics play well with an example like this one? I do not refer 
only to use of raised precision (which in point of fact NSolve will do 
anyway), but also to the underlying arithmetic. What will it mean in 
regards to results, if it can be made to produce results? Will they be 
more accurate?


Here is the story behind the example. Once upon a time, around 
2002-2004, there were some series of posts from a person who went 
variously by "Peter Kosta", "Anonymous Coward", and "Andy Kowar". He 
also posted occasionally under his actual name (and he had been an 
employee here, circa 1992-1997).

He made claims similar to yours, in regard to considering .1 to mean, in 
essence, 1/10. He made them loud and often and about as rudely as 
possible. Unlike yourself, he was entirely familiar with numerical 
analysis, so some of the rants made no sense (he took the arguments 
places you'd be scared to go...)

The irony is that, a few years prior when he was still a reasonably 
pleasant former colleague, he had sent me the above NSolve example 
(under his real name, of course). He was delighted with the results. 
never asked whether I was regarding the decimals as "exact", or having 
fuzz, or whatever.

I don't know, maybe it's my present mood. But I think the example is 
both relevant to the recent thread, and funny given who brought it to my 
attention and subsequent history with same person.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: read CSV
  • Next by Date: Re: "set" data structure in Mathematica? (speeding up graph traversal function)
  • Previous by thread: Re: was: All we have is integers.
  • Next by thread: Re: Plotting problem !