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

**Re: read CSV**

**Re: "set" data structure in Mathematica? (speeding up graph traversal function)**

**Re: was: All we have is integers.**

**Re: Plotting problem !**