an example... Re: Why Mathematica does not issue a warning when the calculations
- To: mathgroup at smc.vnet.net
- Subject: [mg117803] an example... Re: Why Mathematica does not issue a warning when the calculations
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Thu, 31 Mar 2011 04:06:18 -0500 (EST)
Some background from recent notes.
(DL)
>> 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.
(RJF)
> 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.
> [...]
(DL)
>> 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.
(RJF)
> 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.
> ....
(DL)
>> 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.
(RJF)
> 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) +
Sqrt[
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
Out[269]//InputForm=
2.9
In[268]:= 1/(2 Kr) // InputForm
Out[268]//InputForm=
0.3448275862068966
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
2.222...(repeating).
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 -
0.004177758076267676*I,
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.
[Aside:
A variation on NSolve will make an exact computation run in around three
minutes.
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
horrible).
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