A12 input = 5e+5x1+2e-1x2; StringJoin[Characters[input] //. e -> *10^] 5*10^+5x1+2*10^-1x2 ToExpression[%] 500000*x1 + x2/5 Bobby Treat -----Original Message----- should be transcribed into 5 10^5 x1 + 0.2 x2. Chuck ==== Neither SetAccuracy[expr,n] nor SetPrecisions[expr,n] modify expr. These functions modify the prinout not the internal representation. So, the first computation of f is done with approximate numbers and doesn't result in a correct answer due to approximate arithmetic. By assigning a rational expression to each of the variables, you have made them exact numbers and Mathematica responds with an exact solution. >Could someone explain what is going on here, please? > >In[1]:= a = 77617.; b = 33096.; > >In[2]:= SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; > >In[4]:= f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >5.5*b^8 + a/(2*b) > >In[5]:= SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; > >In[6]:= f > >Out[6]= -1.1805916207174113*^21 > >In[7]:= a = 77617; b = 33096; > >In[8]:= g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >(55/10)*b^8 + a/(2*b) > >In[9]:= g > >Out[9]= -(54767/66192) > >In[10]:= N[%] > >Out[10]= -0.8273960599468214 > > >PK > ==== SetAccuracy. However, I still don't understand why the order in which we set the accuracies for f, a, and b matters. In[1]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; a = SetAccuracy[77617., Infinity]; b = SetAccuracy[33096., Infinity]; In[4]:= f Out[4]= -(54767/66192) In[5]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity] Out[5]= 1180591620717411303424 Similarily: In[1]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 50]; a = SetAccuracy[77617., 100]; b = SetAccuracy[33096., 100]; In[4]:= f Out[4]= -0.8273960599468212641107299556`11.4133 In[5]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; Out[5]= 1.180591620717411303424`121.0721*^21 -PK ==== Hallo, I have the problem, that I want to determine the numerical solution of a double integral, like e.g. the following: << Statistics`ContinuousDistributions` ndist = NormalDistribution[0, 1]; N[Sigma^2/(Abs[Mu_1^2 - Mu_2^2]) Integrate[(Integrate[1/y_1 1/(y - y_1) PDF[ndist, d*Log[y_1] + c*Log[y - y_1] + f]* PDF[ndist, Sigma*(c*Log[y_1] + d*Log[y - y_1] + g)], {y_1,0,y}]), {y,0,Lambda}]] whereby Lambda := 3; Mu_1 := 4; Mu_2 := 5; Sigma := 0.25; a = 1/2 (Mu_1^2 + Mu_2^2/Sigma^2) b = 1/2 (Mu_2^2 + Mu_1^2/Sigma^2) c = -Mu_2/(Mu_1^2 - Mu_2^2) d = -Mu_1/(Mu_2^2 - Mu_1^2) f = d*a + c*b g = c*a + d*b Sorry for the poor code... The problem is, that Mathematica doesn't give me a result (after waiting 2 hours I turned my machine off). Thus, the question is, if there is still a possibility of solving such complicated expressions... TIA, Sven. ==== I am trying to solve a system of simultaneous equations with 26 variables and 14 equations (the 12 free variables can be any of the 26 from the eqns.. preferably ones that will minimize the computation time for the other 14). These equation are not linearly related.. the highest degree in any one eqn is degree 4 i believe.. and there are some cross terms in the equations but not every equation depends on every variable.. (some are actually rather simple eqns). Any ideas on how to get started with this using mathematica (any ideas for algorithms).. Anything will be helpful.. I can be reached at ngupta2@seas.upenn.edu Many thanks, Nachi ==== I should add that the solution is over natural numbers.. this will probably make a big difference.. Nachi > I am trying to solve a system of simultaneous equations with 26 variables > and 14 equations (the 12 free variables can be any of the 26 from the eqns.. > preferably ones that will minimize the computation time for the other 14). > These equation are not linearly related.. the highest degree in any one eqn > is degree 4 i believe.. and there are some cross terms in the equations but > not every equation depends on every variable.. (some are actually rather > simple eqns). Any ideas on how to get started with this using > mathematica (any ideas for algorithms).. > > Anything will be helpful.. > > I can be reached at ngupta2@seas.upenn.edu > > Many thanks, > > Nachi > > > > ==== This may seem like a trivial issue, but I find it very frustrating. I use variety of intellectual domains. In every package (JBuilder, KMail, Emacs, XEmacs, Mozilla, xterm, Konsole, etc.) the keyboar mapping is a bit differnt from the other. There are certain idioms which I find to be fairly invaraint between these different packages. I tend to use this common subset more than the package specific idioms. Switching from one package to the next can be a very disorienting experience. It can be even trickier to try and copy and paste from one to the next. I also use 'special' characters in certain domains, ç,?,.81,?,[CapitalYAcute],[CapitalThorn], ß, ¯, etc. Add to all of this, that I run beta code for just about everything. The function of my keyboard changed like the weather. I have spent hours trying to figure out why I can no longer type '.9a'. This doesn't even address the problems of switching between character encodings, or keyboard compose modes. The last thing I want to start doing is messing with the key mappings in my user environment. I want to control the way my keyboard works with Mathematica from within the Mathematica runtime environment. That is, the configuration should be loaded when Mathematica loads, and should not impact the rest of my X environment. If I have come across as a bit jaded regarding this issue, there are reasons. There is a history. I don't find keyboard configuration issues interesting. I want my fine keyboard to just work, the way I want it to Shift+point movement = select text. Ctrl+Insert = copy Shift+Insert = paste Ctrl+c = copy Ctrl+x = cut Ctrl+v = paste Shift+End = select to end of line. etc. Yes, I said I use XEmacs, and Emacs. Yes (X)Emacs is different, but adding yet another alteration with Mathematica is just too much. Is there a way around this? STH ==== >FrontEndExecute[{ > FrontEnd`NotebookWrite[FrontEnd`SelectedNotebook[], > [LeftDoubleBracket][RightDoubleBracket],After]}] FrontEndExecute[{ FrontEnd`NotebookWrite[FrontEnd`SelectedNotebook[], [LeftDoubleBracket][RightDoubleBracket],After], FrontEnd`SelectionMove[FrontEnd`SelectedNotebook[], Previous, Character]}] or FrontEndExecute[{ FrontEnd`NotebookWrite[FrontEnd`SelectedNotebook[], [LeftDoubleBracket][SelectionPlaceholder][RightDoubleBracket], Placeholder]}] -------------------------------------------------------------- Omega Consulting The final answer to your Mathematica needs Spend less time searching and more time finding. http://www.wz.com/internet/Mathematica.html ==== I'm trying to address the special issues related to using Mathematica on the would not if I used a Windows system. These are typically not all that big if a problem _once I figure out what's going on_. What I hope to do is collect all such matters and document them effectively in something like a haven't discussed here http://66.92.149.152/proprietary/com/wri/index.html I'm interested in hearing what you have to say. Of particular interest are the issues faced by a person who is not familiar with the technical aspects would help make this less painful? STH ==== particular, the y-axis label is typically rotated by 90deg so that it reads up the y-axis. This works fine on macs and windoze, but under linux the label runs up the y-axis; however, the letters are rotated so that they are the same orientation as that for the x-axis. Printouts of the NB are fine, but it looks really stupid when using a video projector to teach or give a talk. I have mentioned this before, and the stock answer is that ... it is Mathematica. I personally don't care why the label looks peculiar, just that it does and that there is no easy workaround. -- Kevin J. McCann Joint Center for Earth Systems Technology (JCET) Department of Physics UMBC Baltimore MD 21250 > I'm trying to address the special issues related to using Mathematica on the > would not if I used a Windows system. These are typically not all that big > if a problem _once I figure out what's going on_. What I hope to do is > collect all such matters and document them effectively in something like a which I > haven't discussed here http://66.92.149.152/proprietary/com/wri/index.html > I'm interested in hearing what you have to say. Of particular interest are > the issues faced by a person who is not familiar with the technical aspects > would help make this less painful? > > STH > Reply-To: jmt@dxdydz.net ==== See the MathLink API for C, or the JLink API for java. JLink is easier to use, MathLink is a lower level but native interface. JLink is built on MathLink. As far as I know, the perl API Math::ematica has not been upgraded to Mathematica 4. > > Is it possible to get a document of the description of how to > interface with the kernel? Kind of what should an interface say to the > kernel and how to connect to it. > > > Paulo > > ==== Is it possible to get a document of the description of how to interface with the kernel? Kind of what should an interface say to the kernel and how to connect to it. Paulo ==== for your extensive answer but I still have some doubts about convergence of the following integral (m,n integrers>=0) W[m_,n_]:=Integrate[BesselJ[m, x]*BesselJ[n, x], {x, 0, Infinity}] for wich Mathematica gives the close form W[m_,n_]:= -Cos[(m-n)Pi/2]/(2 Pi)* ( 2 EulerGamma + Log[4] + PolyGamma[0, 1/2(1 + m - n)] + PolyGamma[0, 1/2(1 - m + n)] + 2PolyGamma[0, 1/2(1 + m + n)] ) You say this integral is convergent to 1/2 for m=0 and n=1. Also Mathematica agrees to you since for m>=0 W[m,m+1]=1/2 W[m,m+3]=-1/2 Numerically we have NIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, 0, Infinity}] NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy.... 0.597973 NIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, 0, Infinity}, Method -> Oscillatory] NIntegrate::ploss : .... 0.5 So I define also the corresponding numeric definition NW[m_, n_] := NIntegrate[BesselJ[m, x]*BesselJ[n, x], {x, 0, Infinity}, Method -> Oscillatory] THEORY The integral is the critical case of Weber-Schafheitlin integral (see Watson book on Bessel function p.402, or Ryzhik-Gradshteyn 6.574(2)). According to this theory WS[m_,n_,p_]:=Integrate[BesselJ[m, x]*BesselJ[n, x] x^-p, {x, 0, Infinity}] = A/B where A=Gamma[p]*Gamma[(n+m-p+1)/2] B=2^p Gamma[(n-m+p+1)/2]Gamma[(n+m+p+1)/2]Gamma[(m-n+p+1)/2] By the presence of Gamma[p] in numerator A, in the case p=0 as in W[m,n] all these integrals are divergent since Gamma[0]=Infinity. The integral exist if m+n+1 > p > 0. ASYMPTOTICS The Watson theory is in conflict with Mathematica and your notes according which the asyntotic trend 1/x of the integrand in W[m,n] is enough for convergernce. I divide the integral in two parts Wasy[m_,n_,a_]=NIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, 0, a]+ NIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, a, Infinity}] and if a>>1 I use asyntotic expansion of Bessel function in the second integral so that I can write Wasy[m_,n_,a_]= int1[m,n,a]+int2[m,n,a] where int1[m_,n_,a_]:=NIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, 0, a]+ int2[m_,n_,a_]:=(2/Pi)Integrate[Cos[x-(2m+1)Pi/4]*Cos[x-(2n+1)Pi/4], {x, a, Infinity}] The first integral is a quite normal finite integral. The second (int2) is singular and according to Mathematica 4.1 int2[m_, n_, a_] := -(1/Pi)*Log[a]*Cos[1/2(m - n)Pi]*]Log[a] + (1/Pi)*CosIntegral[2 a]*Sin[1/2(m+n)Pi] + 1/(2*Pi)*Cos[1/2(m+n)Pi]*(Pi-SinIntegral[2*a]) RESULTS m=1;n=0;a=20.; WS[m,n,0]=divergent W[m,n]=1/2 NW[m,n]=0.5 Wasy[m,n,a]=.49816 m=2;n=0;a=20.; WS[m,n,0]=divergent W[m,n]=0.427599 NW[m,n]=-2.43818 Wasy[m,n,a]=-1.48052 m=3;n=1;a=20.; WS[m,n,0]=divergent W[m,n]=0.639806 NW[m,n]=-2.31957 Wasy[m,n,a]=-1.26822 m=4;n=0;a=20.; WS[m,n,0]=divergent W[m,n]=-.852012 NW[m,n]=1.45786 Wasy[m,n,a]=1.06835 The cases W[m,m+1],W[m,m+3] well agrre with the numerical counterpart. Other case are doubtfully. I think the main problem is the convergence of this kind of integrals. Any suggestion will be well considerd. Robert Roberto Brambilla CESI Via Rubattino 54 20134 Milano tel +39.02.2125.5875 fax +39.02.2125.5492 rlbrambilla@cesi.it ==== On Sun, 29 Sep 2002 09:35:41, in the message Re: A Bessel integral, VB>> The expression for W[m_,n_] returned by Mathematica is wrong. VB>> VB>> To prove, just substitute m = n = 0 which is exactly what you had done VB>> VB>> and observe that the output you had had VB>> VB>> W[0,0]=-(2 EulerGamma + Log[4] + 4 PolyGamma[0, 1/2])/(2 Pi) VB>> VB>> = 0.84564 VB>> VB>> was incorrect. The correct answer is 1/2. ^^^^^^^^^^^^^^^^^^^^^^^^^^ TB> W[0,0]diverges. Mathematica gets that wrong. (That my terrible bug shows how it is dangerous to do several posting to the MathGroup before sending them ;-) Why sure, you are right, the integral Integrate[BesselJ[0, z]^2, {z, 0, Infinity}] diverges because the integrand is bounded everywhere over the integration region and decays at z -> Infinity as Cos[Pi/4 - z]^2/z + o(z), that is as In[1] := Expand[TrigExpand[Cos[Pi/4 - z]^2/z]] // InputForm Out[1] = 1/(2*z) + (Cos[z]*Sin[z])/z which means that the integral Integrate[BesselJ[0, z]^2, {z, 0, x}] diverges logarithmically in x. By the way, the main term of Expand[Normal[Series[BesselJ[0, z], {z, Infinity, 1}]]^2] is (2*Cos[Pi/4 - z]^2)/(Pi*z) which conveys the suggestion that we should try it, too. This reveals us another integral which Mathematica 4.1 fails to calculate In[1] := Integrate[Cos[Pi/4 - z]^2/z, {z, 1, Infinity}] // N Out[1]= -0.0173083 Out[2]= 4.1 for Microsoft Windows (November 2, 2000) but Mathematica 4.2 handles correctly In[1] := Integrate[Cos[Pi/4 - z]^2/z, {z, 1, Infinity}] Out[1] = Integrate::idiv: Integral of... does not converge on {1, Infinity). Out[2]= 4.2 for Microsoft Windows (February 28, 2002) Even simpler, In[1] := Integrate[Cos[z]^2/z, {z, 1, Infinity}] Out[1] = -EulerGamma/2 - Log[2]/2 + (EulerGamma - CosIntegral[2] + Log[2])/2 Out[2]= 4.1 for Microsoft Windows (November 2, 2000) which is wrong while Mathematica 4.2 works excellent In[1] := Integrate[Cos[z]^2/z, {z, 1, Infinity}] Out[1] = Integrate::idiv: Integral of...does not converge on {1, Infinity). Out[2]= 4.2 for Microsoft Windows (February 28, 2002) Best wishes, Vladimir Bondarenko Mathematical Director Symbolic Testing Group Web : http://www.CAS-testing.org/ http://maple.bug-list.org/VER2/ (under tuning) http://maple.bug-list.org/VER3/ (under tuning) http://maple.bug-list.org/VER1/ (under tuning) http://www.beautyriot.com/ (teamwork) http://www.ohaha.com/ (teamwork) Voice: (380)-652-447325 Mon-Fri 9 a.m. - 6 p.m. Mail : 76 Zalesskaya Str, Simferopol, Crimea, Ukraine ==== John, You could do something like this. points = With[{del = 2Pi/24}, Table[{Cos[t], Sin[t], t/(2Pi)}, {t, 0, 2Pi - del, del}]]; Show[Graphics[ {AbsolutePointSize[7], {Hue[Last[#]], Point[#[[{1, 2}]]]} & /@ points, Line[Drop[#, -1] & /@ points]}], AspectRatio -> Automatic, Background -> GrayLevel[0.4], ImageSize -> 400]; When I made the point list I made certain the z values were between 0 and 1. Otherwise you will have to define a color function that will associate a proper color with each value of z. David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Approved: Steven M. Christensen , Moderator ==== Many thanks to all who replied. The original problem was as follows: In a presentation I wish to use Plot to generate a sequence of frames and then animate them. The problem is that the audience sees the animation twice. Once when the frames are being generated and then again after I have closed the group and double clicked on the top graphic. However, the first showing during generation is enough (but uncontrolled). Is it possible to tidy up the generation of the graphic so that it becomes the animation? There were two main solutions which I now apply to my real problem and not the simple, generic, problem in the original post. The real problem was how to build up a probably density function (PDF) in real time. In the examples below I redraw the PDF every time I add 10 points. Initialise <500] ], {500}]; NotebookWrite[EvaluationNotebook[],Cell[CellGroupData[graphs,Closed]]]; SelectionMove[EvaluationNotebook[], All, GeneratedCell]; FrontEndExecute[{FrontEndToken[EvaluationNotebook[], SelectionAnimate]}] ] This solution works but it generates 500 frames and sometimes exceeds the memory. The paradigm here is generate all frames, then animate all frames. We really need a loop that does: generate next frame, delete last frame, show next frame Is it possible to do this? Bobby Treat also offered a solution involving SelectionMove. The second solution was from Jerry Blimbaum and uses JAVA InstallJava[]; UseFrontEndForRendering = False; createWindow[] := Module[{frame}, frame = JavaNew[com.wolfram.jlink.MathFrame, Probability Density Function]; drawArea = JavaNew[com.wolfram.jlink.MathCanvas]; drawArea@setUsesFE[UseFrontEndForRendering]; drawArea@setSize[800, 600]; JavaBlock[frame@setLayout[JavaNew[java.awt.BorderLayout]]; frame@add[drawArea, ReturnAsJavaObject[BorderLayout`CENTER]]; frame@pack[]; frame@setSize[800, 600]; frame@setLocation[100, 100]; JavaShow[frame]]; frame ] ClearAll[drawString]; drawString[] :=( data=Flatten[Join[data,RandomArray[wb,10]]]; freq=BinCounts[data,{0,50,1}]; BarChart[Transpose[{freq,midpts}],ImageSize ->500, DisplayFunction -> Identity]) LoadJavaClass[java.lang.Thread]; AnimationPlot[n_] := JavaBlock[ Block[{frm}, frm = createWindow[]; Do[(obj = drawString[]; drawArea@setMathCommand[obj]; drawArea@repaintNow[]; Thread@sleep[];) ,{n} ]]] data={}; AnimationPlot[500]; This solution works and does not use significant memory. However, we have not managed to control the speed of this animation. The JAVA code sleep does not work nor does the use of a Mathematica Pause which always rounds up to an integer (is this a bug?). Hugh Goyder ==== John, It is easier to do without ListPlot: Make some data, dat= Table[Random[],{10},{3}]; Show[Graphics[{PointSize[.05],{Hue[2/3#3],Point[{#1,#2}]}&@@@dat} ], Frame->True ]; If the points need to be joined then something like Show[Graphics[{ Line /@Partition[dat[[All,{1,2}]],2,1], PointSize[.05],{Hue[2/3#3],Point[{#1,#2}]}&@@@dat } ], Frame->True ]; -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > I have a list of points l1={xi,yi,zi} how can I make a 2D list plot of > > ==== A few years ago we made a package that do just this. See http://cern.ch/jowett/Mathematica/Graphics/ColorListPlot.html The type of plot you want is covered in the section Examples then Three-dimensional data. The Introduction gives links for downloading the package. JMJ > I have a list of points l1={xi,yi,zi} how can I make a 2D list plot of > > ==== Your problem can be solved in numerous ways of course, try something like .: Module[ { data=Flatten[Table[{x,y,Random[]},{x,10},{y,10}],1] }, Show[ Graphics[ { AbsolutePointSize[10], data/.{x_,y_,z_}[Rule]{Hue[z],Point[{x,y}]} } ] ,AspectRatio[Rule]Automatic ] ] Note that you can replace the main part, ie the transforming rule from point lists to colored point directives with a function for example .: {Hue[#3],Point[{#1,#2}]}&@@@data bye, Borut | I have a list of points l1={xi,yi,zi} how can I make a 2D list plot of Reply-To: kuska@informatik.uni-leipzig.de ==== data = Table[{Random[], Random[], Random[]}, {20}]; Show[Graphics[ {Hue[Last[#]], Point[Take[#, 2]]} & /@ data, Axes -> True ] ] Jens > > I have a list of points l1={xi,yi,zi} how can I make a 2D list plot of ==== >>I need a loop that goes generate next frame, delete old frame, show new frame so that the number of frames does not become excessive. I'm pinging the group for that. I'm just following along in this, other than the trick of taking advantage of the half-period. I'll be very interested in a solution myself, as I frequently run out of memory in animations of fairly modest size -- despite having 1024MB of RAM. Bobby -----Original Message----- need a loop that goes generate next frame, delete old frame, show new frame so that the number of frames does not become excessive. Any ideas? Hugh Goyder -----Original Message----- graphs = Rest@Join[half, Rest@Reverse@half]; NotebookWrite[EvaluationNotebook[], Cell[CellGroupData[graphs, Closed]]]; SelectionMove[EvaluationNotebook[], All, GeneratedCell]; FrontEndExecute[{FrontEndToken[EvaluationNotebook[], SelectionAnimate]}]] Bobby -----Original Message----- >showing during generation is enough (but uncontrolled). > >Hugh Goyder This creates a graphics cell from a graphics expression. GraphicCell[graphics_] := Cell[GraphicsData[PostScript, DisplayString[graphics]],Graphics] cellgroup. Block[{$DisplayFunction=Identity, graphs}, graphs = Table[GraphicCell[ Plot[Sin[t]*Sin[x], {x, 0, Pi}, PlotRange -> {{0, Pi}, {-1, 1}}, ImageSize -> 400]], {t,0,15,.1}]; NotebookWrite[EvaluationNotebook[],Cell[CellGroupData[graphs,Closed]]]; SelectionMove[EvaluationNotebook[], All, GeneratedCell]; FrontEndExecute[{FrontEndToken[EvaluationNotebook[], SelectionAnimate]}] ] -------------------------------------------------------------- Omega Consulting The final answer to your Mathematica needs Spend less time searching and more time finding. http://www.wz.com/internet/Mathematica.html ==== Actually, we don't know whether SetAccuracy succeeds, because we don't know how inexact those numbers really are. If they are known to more digits than shown in the original post, they should be entered with as much precision as they deserve. If not, there's no trick or algorithm that will give the result precision, because the value of f really is in the noise. That is, we have no idea what the value of f should be. Mathematica's failing is in returning a value without pointing out that it has no precision. Bobby -----Original Message----- f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/(2*b), Infinity]; a=77617.1; b=33096.1; a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; f -1564032114908486835197494923989618867972540153873951942813115514949 3891236234 525007719168693704591197760187988046304361497869199129319625743010292363 1246 75 /10867106143970760551000357827554793888198143135975649579607989867743572 8240 16 0653953612982932181371232436367739737604096 2) Rewriting as fractions a=776171/10; b=330961/10; f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) -(5954133808997234115690303589909929091649391296257/ 41370125000000) 3) Using Rationalize Clear[a,b,f] f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/(2*b), 0]; a=77617.1; b=33096.1; a=Rationalize[a,0];b=Rationalize[b,0]; f -(5954133808997234115690303589909929091649391296257/ 41370125000000) I use Rationalize[. , 0] besause of results like Rationalize[3.1415959] 3.1416 Rationalize[3.1415959,0] 31415959/10000000 -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > Well, first of of all, your using SetAccuracy and SetPrecision does > nothing at all here, since they do not change the value of a or b. You > should use a = SetAccuracy[a, Infinity] etc. But even then you won't > get the same answer as when you use exact numbers because of the way > you evaluate f. Here is the order of evaluation that will give you the > same answer, and should explain what is going on: > > f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* > b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; > > a = 77617.; > > > b = 33096.; > > a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; > > f > > 54767 > -(-----) > 66192 > > Andrzej Kozlowski > Toyama International University > JAPAN > > > > > > Could someone explain what is going on here, please? > > > > In[1]:= > > a = 77617.; b = 33096.; > > > > In[2]:= > > SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; > > SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; > > > > In[4]:= > > f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + > > a/(2*b) > > > > In[5]:= > > SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; > > > > In[6]:= > > f > > > > Out[6]= > > -1.1805916207174113*^21 > > > > In[7]:= > > a = 77617; b = 33096; > > > > In[8]:= > > g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + > > (55/10)*b^8 + a/(2*b) > > > > In[9]:= > > g > > > > Out[9]= > > -(54767/66192) > > > > In[10]:= > > N[%] > > > > Out[10]= > > -0.8273960599468214 > > > > > > PK > > > > > > > > ==== It seems clear to me that what Allan and what you mean by succeeds here refer to quite different things and your objection is therefore beside the point. There are obviously two ways in which one can interpret the original posting. The first interpretation, which Allan and myself adopted, was that the question concerned purely the computational mechanism of Mathematica. Or, to put it in other words, it was why are the results of these two computations not the same?. In this sense success means no more than making Mathematica return the same answer using the two different routes the original poster used. You on the other hand choose to assume that the posting shows that its author does not understand not just the mechanism of significance arithmetic used by Mathematica but also computations with inexact numbers in general. I do not think this is necessarily the correct assumption. I also don't see that Mathematica is doing anything wrong. After all, one can always check the accuracy of your answer: In[1]:= a = 77617.; b = 33096.; In[2]:= f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b) In[3]:= f Out[3]= -1.1805916207174113*^21 In[4]:= Accuracy[%] Out[4]= -5 which tells you that it can't be very reliable. What more do you demand? Andrzej Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ > Actually, we don't know whether SetAccuracy succeeds, because we > don't > know how inexact those numbers really are. If they are known to more > digits than shown in the original post, they should be entered with as > much precision as they deserve. If not, there's no trick or algorithm > that will give the result precision, because the value of f really is > in the noise. That is, we have no idea what the value of f should > be. > > Mathematica's failing is in returning a value without pointing out that > it has no precision. > > Bobby > > -----Original Message----- > Sent: Monday, September 30, 2002 11:59 AM > > Andrzej, Bobby, Peter > > It looks as if using SetAccuracy succeeds here because the inexact > numbers > that occur have finite binary representations. If we change them > slightly to > avoid this then we have to use Rationalize: > > 1) Using SetAccuracy > > Clear[a,b,f] > > > f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ > (2*b), > Infinity]; > > a=77617.1; > b=33096.1; > > a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; > > f > > > -1564032114908486835197494923989618867972540153873951942813115514949 > 3891236234 > > 52500771916869370459119776018798804630436149786919912931962574301029236 > 3 > 1246 > 75 > > / > 10867106143970760551000357827554793888198143135975649579607989867743572 > 8240 > 16 > 0653953612982932181371232436367739737604096 > > 2) Rewriting as fractions > > a=776171/10; > b=330961/10; > > f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) > > -(5954133808997234115690303589909929091649391296257/ > 41370125000000) > > 3) Using Rationalize > > Clear[a,b,f] > > > f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ > (2*b), > 0]; > > a=77617.1; > b=33096.1; > > a=Rationalize[a,0];b=Rationalize[b,0]; > > f > > -(5954133808997234115690303589909929091649391296257/ > 41370125000000) > > > I use Rationalize[. , 0] besause of results like > > Rationalize[3.1415959] > > 3.1416 > > Rationalize[3.1415959,0] > > 31415959/10000000 > -- > Allan > > --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 > > >> Well, first of of all, your using SetAccuracy and SetPrecision does >> nothing at all here, since they do not change the value of a or b. You >> should use a = SetAccuracy[a, Infinity] etc. But even then you won't >> get the same answer as when you use exact numbers because of the way >> you evaluate f. Here is the order of evaluation that will give you the >> same answer, and should explain what is going on: >> >> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* >> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; >> >> a = 77617.; >> >> >> b = 33096.; >> >> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; >> >> f >> >> 54767 >> -(-----) >> 66192 >> >> Andrzej Kozlowski >> Toyama International University >> JAPAN >> >> >> >> >>> Could someone explain what is going on here, please? >>> >>> In[1]:= >>> a = 77617.; b = 33096.; >>> >>> In[2]:= >>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; >>> >>> In[4]:= >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + >>> a/(2*b) >>> >>> In[5]:= >>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; >>> >>> In[6]:= >>> f >>> >>> Out[6]= >>> -1.1805916207174113*^21 >>> >>> In[7]:= >>> a = 77617; b = 33096; >>> >>> In[8]:= >>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>> (55/10)*b^8 + a/(2*b) >>> >>> In[9]:= >>> g >>> >>> Out[9]= >>> -(54767/66192) >>> >>> In[10]:= >>> N[%] >>> >>> Out[10]= >>> -0.8273960599468214 >>> >>> >>> PK >>> >>> >>> >> >> > > > > > > > > ==== > It seems clear to me that what Allan and what you mean by succeeds > here refer to quite different things and your objection is therefore > beside the point. There are obviously two ways in which one can > interpret the original posting. The first interpretation, which Allan > and myself adopted, was that the question concerned purely the > computational mechanism of Mathematica. Or, to put it in other words, > it was why are the results of these two computations not the same?. > In this sense success means no more than making Mathematica return > the same answer using the two different routes the original poster used. > You on the other hand choose to assume that the posting shows that its > author does not understand not just the mechanism of significance > arithmetic used by Mathematica but also computations with inexact > numbers in general. I do not think this is necessarily the correct > assumption. I also don't see that Mathematica is doing anything wrong. > After all, one can always check the accuracy of your answer: > > In[1]:= > a = 77617.; b = 33096.; > > In[2]:= > f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + > 5.5*b^8 + a/(2*b) > > In[3]:= > f > > Out[3]= > -1.1805916207174113*^21 > > In[4]:= > Accuracy[%] > > Out[4]= > -5 > > which tells you that it can't be very reliable. What more do you demand? > As you are dealing here only with machine-precision numbers, your When you do calculations with arbitrary-precision numbers, as discussed in the previous section, Mathematica always keeps track of the precision of your results, and gives only those digits which are known to be correct, given the precision of your input. When you do calculations with machine-precision numbers, however, Mathematica always gives you a machine-[CapitalEth]precision result, whether or not all the digits in the result can, in fact, be determined to be correct on the basis of your input. In practice, to rely on a numerical result, you ALWAYS have to check its accuracy. How reliable is Accuracy anyway? In[1]:= a = SetAccuracy[77617., Infinity]; b = SetAccuracy[33096., Infinity]; In[3]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity] Out[3]= 1180591620717411303424 In[4]:= Accuracy[f] Out[4]= Infinity We got completely wrong result with Infinite accuracy. In conclusion, one can not rely on Accuracy. Checking numerical results in Mathematica sounds like a tough task.:-) --PK > Andrzej > > > Andrzej Kozlowski > Yokohama, Japan > http://www.mimuw.edu.pl/~akoz/ > http://platon.c.u-tokyo.ac.jp/andrzej/ > [...] > > ==== Andrzej Yes, like you I took the original question to be about how to get the result of using the naive rational versions in place of the inexact numbers. Bobby raises the question of how we might know accuracy of the result. You answer this with > a = 77617.; b = 33096.; > > In[2]:= > f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + > 5.5*b^8 + a/(2*b) > > In[3]:= > f > > Out[3]= > -1.1805916207174113*^21 > > In[4]:= > Accuracy[%] > > Out[4]= > -5 However this computation is done in machine arithmetic, which means that Mathematica keeps no check on the accuracy and precision of the computation, and Mathematica gives the default accuracy value without any real signifcance: $MachinePrecision - Log[10,Abs[f]]//Round -5 To get meaningful accuracy and precision values we need to force the computation to be in bignums (bigfloat, arbitrary precision) arithmetic. Mathematica then keeps track of the accuracy and precision that it can guarantee. Clear[f,a,b,k] k=50; f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/(2*b),k]; a=77617.;b=33096.; a=SetAccuracy[a,k]; b=SetAccuracy[b,k]; f -0.82739605995 Accuracy[f] 11 Precision[f] 11 Which tells us that the error in the computed value of f is not more than 10^-11 The above results are rounded. More accurately we get Accuracy[f, Round->False] 11.4956 Precision[f, Round->False] 11.4133 There are several related issues here: - is the precision (accuracy) of rhe input known? - how does Mathematica interpret what one gives it? - how does Mathematica go about the calculation? -- Allan --------------------- Allan Hayes Mathematica Training and Consulting Leicester UK www.haystack.demon.co.uk hay@haystack.demon.co.uk Voice: +44 (0)116 271 4198 > It seems clear to me that what Allan and what you mean by succeeds > here refer to quite different things and your objection is therefore > beside the point. There are obviously two ways in which one can > interpret the original posting. The first interpretation, which Allan > and myself adopted, was that the question concerned purely the > computational mechanism of Mathematica. Or, to put it in other words, > it was why are the results of these two computations not the same?. > In this sense success means no more than making Mathematica return > the same answer using the two different routes the original poster used. > You on the other hand choose to assume that the posting shows that its > author does not understand not just the mechanism of significance > arithmetic used by Mathematica but also computations with inexact > numbers in general. I do not think this is necessarily the correct > assumption. I also don't see that Mathematica is doing anything wrong. > After all, one can always check the accuracy of your answer: > > In[1]:= > a = 77617.; b = 33096.; > > In[2]:= > f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + > 5.5*b^8 + a/(2*b) > > In[3]:= > f > > Out[3]= > -1.1805916207174113*^21 > > In[4]:= > Accuracy[%] > > Out[4]= > -5 > > which tells you that it can't be very reliable. What more do you demand? > > Andrzej > > > Andrzej Kozlowski > Yokohama, Japan > http://www.mimuw.edu.pl/~akoz/ > http://platon.c.u-tokyo.ac.jp/andrzej/ > > > > Actually, we don't know whether SetAccuracy succeeds, because we > > don't > > know how inexact those numbers really are. If they are known to more > > digits than shown in the original post, they should be entered with as > > much precision as they deserve. If not, there's no trick or algorithm > > that will give the result precision, because the value of f really is > > in the noise. That is, we have no idea what the value of f should > > be. > > > > Mathematica's failing is in returning a value without pointing out that > > it has no precision. > > > > Bobby > > > > -----Original Message----- > > Sent: Monday, September 30, 2002 11:59 AM > > > > Andrzej, Bobby, Peter > > > > It looks as if using SetAccuracy succeeds here because the inexact > > numbers > > that occur have finite binary representations. If we change them > > slightly to > > avoid this then we have to use Rationalize: > > > > 1) Using SetAccuracy > > > > Clear[a,b,f] > > > > > > f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ > > (2*b), > > Infinity]; > > > > a=77617.1; > > b=33096.1; > > > > a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; > > > > f > > > > > > -1564032114908486835197494923989618867972540153873951942813115514949 > > 3891236234 > > > > 52500771916869370459119776018798804630436149786919912931962574301029236 > > 3 > > 1246 > > 75 > > > > / > > 10867106143970760551000357827554793888198143135975649579607989867743572 > > 8240 > > 16 > > 0653953612982932181371232436367739737604096 > > > > 2) Rewriting as fractions > > > > a=776171/10; > > b=330961/10; > > > > f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) > > > > -(5954133808997234115690303589909929091649391296257/ > > 41370125000000) > > > > 3) Using Rationalize > > > > Clear[a,b,f] > > > > > > f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ > > (2*b), > > 0]; > > > > a=77617.1; > > b=33096.1; > > > > a=Rationalize[a,0];b=Rationalize[b,0]; > > > > f > > > > -(5954133808997234115690303589909929091649391296257/ > > 41370125000000) > > > > > > I use Rationalize[. , 0] besause of results like > > > > Rationalize[3.1415959] > > > > 3.1416 > > > > Rationalize[3.1415959,0] > > > > 31415959/10000000 > > -- > > Allan > > > > --------------------- > > Allan Hayes > > Mathematica Training and Consulting > > Leicester UK > > www.haystack.demon.co.uk > > hay@haystack.demon.co.uk > > Voice: +44 (0)116 271 4198 > > > > > >> Well, first of of all, your using SetAccuracy and SetPrecision does > >> nothing at all here, since they do not change the value of a or b. You > >> should use a = SetAccuracy[a, Infinity] etc. But even then you won't > >> get the same answer as when you use exact numbers because of the way > >> you evaluate f. Here is the order of evaluation that will give you the > >> same answer, and should explain what is going on: > >> > >> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* > >> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; > >> > >> a = 77617.; > >> > >> > >> b = 33096.; > >> > >> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; > >> > >> f > >> > >> 54767 > >> -(-----) > >> 66192 > >> > >> Andrzej Kozlowski > >> Toyama International University > >> JAPAN > >> > >> > >> > >> > >>> Could someone explain what is going on here, please? > >>> > >>> In[1]:= > >>> a = 77617.; b = 33096.; > >>> > >>> In[2]:= > >>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; > >>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; > >>> > >>> In[4]:= > >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + > >>> a/(2*b) > >>> > >>> In[5]:= > >>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; > >>> > >>> In[6]:= > >>> f > >>> > >>> Out[6]= > >>> -1.1805916207174113*^21 > >>> > >>> In[7]:= > >>> a = 77617; b = 33096; > >>> > >>> In[8]:= > >>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + > >>> (55/10)*b^8 + a/(2*b) > >>> > >>> In[9]:= > >>> g > >>> > >>> Out[9]= > >>> -(54767/66192) > >>> > >>> In[10]:= > >>> N[%] > >>> > >>> Out[10]= > >>> -0.8273960599468214 > >>> > >>> > >>> PK > >>> > >>> > >>> > >> > >> > > > > > > > > > > > > > > > > > > > ==== These expressions are condensation of larger ones (about 700 lines or so each) but they illustrate random substitution failures in 4.2. Question: how can the substitution Sqrt[...]->Q always be made to work? The help file under ReplaceAll, ReplaceRepeated, etc, does not address this problem. f=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]/4*F^2); Print[(f/.Sqrt[X+Y+Z]->Q)//InputForm]; B*(A + Q) + (4*C)/(F^2*Sqrt[X + Y + Z]) (* fails *) g=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]*4*F^2); Print[(g/.Sqrt[X+Y+Z]->Q)//InputForm]; B*(A + Q) + C/(4*F^2*Sqrt[X + Y + Z]) (* fails *) h=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]+4*F^2); Print[(h/.Sqrt[X+Y+Z]->Q)//InputForm]; B*(A + Q) + C/(4*F^2 + Q) (* works *) Reply-To: kuska@informatik.uni-leipzig.de ==== Sqrt[a] is internal Power[a,Rational[1,2]] and 1/Sqrt[a] is interal Power[a,Rational[-1,2]] and so a rule Sqrt[a]->q will not match with 1/Sqrt[a]. You need f = B*(A + Sqrt[X + Y + Z]) + C/(Sqrt[X + Y + Z]/4*F^2); (f /. (X + Y + Z)^Rational[n_, 2] -> Q^n) Jens > > These expressions are condensation of larger ones > (about 700 lines or so each) but they illustrate random > substitution failures in 4.2. Question: how can the > substitution Sqrt[...]->Q always be made to work? > The help file under ReplaceAll, ReplaceRepeated, etc, > does not address this problem. > > > f=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]/4*F^2); > Print[(f/.Sqrt[X+Y+Z]->Q)//InputForm]; > > B*(A + Q) + (4*C)/(F^2*Sqrt[X + Y + Z]) (* fails *) > > g=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]*4*F^2); > Print[(g/.Sqrt[X+Y+Z]->Q)//InputForm]; > > B*(A + Q) + C/(4*F^2*Sqrt[X + Y + Z]) (* fails *) > > h=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]+4*F^2); > Print[(h/.Sqrt[X+Y+Z]->Q)//InputForm]; > > B*(A + Q) + C/(4*F^2 + Q) (* works *) ==== > > > Sqrt[a] is internal Power[a,Rational[1,2]] and > 1/Sqrt[a] is interal Power[a,Rational[-1,2]] > and so a rule Sqrt[a]->q will not match with > 1/Sqrt[a]. You need > > f = B*(A + Sqrt[X + Y + Z]) + C/(Sqrt[X + Y + Z]/4*F^2); > (f /. (X + Y + Z)^Rational[n_, 2] -> Q^n) > > Jens ought to be in the help Examples under . and . I should clarify three things. First, why use of FullForm is impractical. The source expressions I am dealing with are highly complex, with thousands of terms. The subexpressions to be replaced appear in hundreds of places, in many nested combinations. Detailed eye examination after each run would take a long time. Second, the operation subexpression->letter is used as preparation for conversion of those expressions to C. The letters will be names of temps (temporary variables) in C code. Upon replacing all subexpressions the source contracts to about 1-5% of original size. Third, there is a brute force replacement method which can be used (and was): output in InputForm, save cell as text, use a smart text editor that ignores blanks and CRs (e.g. emacs), and re-input for C conversion. This is cumbersome (the text has to be telnet'd to a Unix machine and back) and error prone but available as last resort. Perhaps a future version of Mathematica may incorporate a <- operator for this kind of reverse expansion to extract common subexpressions. The output would be the contracted expression and a temps list. ==== Dear Colleagues, I calculated: Sum[1/Prime[n], {n, 15000}] // N Result: 2.74716 Now I wonder if this sum will converge or keep on growing, albeit very slowly. Matthias Bode Sal. Oppenheim jr. & Cie. KGaA Koenigsberger Strasse 29 D-60487 Frankfurt am Main GERMANY Mobile: +49(0)172 6 74 95 77 Internet: http://www.oppenheim.de ==== > I calculated: > > Sum[1/Prime[n], {n, 15000}] // N > > Result: 2.74716 > > Now I wonder if this sum will converge or keep on growing, albeit very > slowly. The latter. It is a well known fact that the series diverges (which, by the way, shows that there are infinitely many primes :-). David -- -------------------- http://NewsReader.Com/ -------------------- Usenet Newsgroup Service ==== Hm, I wonder whose failures are you referring to when write substitution failures in 4.2? When using Mathematica's pattern matching there is one fundamental rule (very frequently restated on this list) you should adhere to: check the FullForm of the expression you are trying to match. So taking just your first case: In[1]:= FullForm[f = B*(A + Sqrt[X + Y + Z]) + C/((Sqrt[X + Y + Z]/4)*F^2)] Out[2]//FullForm= Plus[Times[4,C,Power[ F,-2],Power[Plus[X,Y,Z],Rational[-1,2]]],Times[B,Plus[A,Power[Plus[X, Y,Z],Rational[1,2]]]]] This ought to make the reason for the failure of your substitution clear. To make it work you must find a way to much the right pattern and not forget that the matching is purely syntactic. In[2]:= f /. (X + Y + Z)^(Rational[x_, 2]) -> Q^x Out[2]= (4*C)/(F^2*Q) + B*(A + Q) There are of course other ways you can get this to work, e.g. In[3]:= PowerExpand[f /. X + Y + Z -> Q^2] Out[3]= (4*C)/(F^2*Q) + B*(A + Q) There is even a rather crazy method that sometimes actually works: In[4]:= ToExpression[StringReplace[ToString[Evaluate[InputForm[f]]], Sqrt[X + Y + Z] -> Q]] Out[4]= (4*C)/(F^2*Q) + B*(A + Q) Andrzej Kozlowski Toyama International University JAPAN http://sigma.tuins.ac.jp/~andrzej/ > These expressions are condensation of larger ones > (about 700 lines or so each) but they illustrate random > substitution failures in 4.2. Question: how can the > substitution Sqrt[...]->Q always be made to work? > The help file under ReplaceAll, ReplaceRepeated, etc, > does not address this problem. > > > > f=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]/4*F^2); > Print[(f/.Sqrt[X+Y+Z]->Q)//InputForm]; > > B*(A + Q) + (4*C)/(F^2*Sqrt[X + Y + Z]) (* fails *) > > g=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]*4*F^2); > Print[(g/.Sqrt[X+Y+Z]->Q)//InputForm]; > > B*(A + Q) + C/(4*F^2*Sqrt[X + Y + Z]) (* fails *) > > h=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]+4*F^2); > Print[(h/.Sqrt[X+Y+Z]->Q)//InputForm]; > > B*(A + Q) + C/(4*F^2 + Q) (* works *) > > > ==== I was trying to run one of the Mathematica Book -> Graphics Gallery -> Animations -> Rolling Square. I don't recall the exact sequence of actions I took. I believe I selected the cell with the square (the only cell in the notebook) and from the menu, Cell -> Animate Selected Graphics. This resulted in a set of animation control buttons appearing in the bottom frame of the window. I clicked on one of these buttons, but nothing happened. I looke back in the menu and saw M-y as a keyboard shorcut to run an animation. I tried that with no result. I clicked another button in graphics control set, and my X windows locked up. This included the keyboard's ability to give me another display by using Ctl+Alt+F1. I went to another system and ssh-ed in and found Mathematica had over 50% of my user resources, and was climbing. The same was true for VM. I have a gig of physical RAM. Once I killed Mathematica, my X came back to life. I've had several bad experiences with Mathematica and X. I honestly believe there isolate and fix these. Have others had such problems? STH ==== > XFree86 4.0.2. I was trying to run one of the Mathematica Book -> Graphics > Gallery -> > Animations -> Rolling Square. I don't recall the exact sequence of > actions > I took. I never had these problems. But I'm not longer able to export gif's etc. -- Hendrik van Hees Fakult.8at f.9fr Physik http://theory.gsi.de/~vanhees/ D-33615 Bielefeld ==== Dear MathGroup Members, I want to minimize a function which returns the minimizing value (arg min) of another function. For a simple example consider the following function opt which returns the arg min of x-2.5(1+Erf[x-s]). opt[s_]:=Block[{x}, x/. Last[ FindMinimum[x-2.5(1+Erf[x-s]), {x,1,3}]]] Now in a second step I want (again this is only a simple example for illustrative purposes) to minimize (opt[s]-2)^2 with respect to s. FindMininum has no problems with this. FindMinimum[(opt[s]-2)^2,{s,0.9,1.1}] {3.18689*^-23, {s -> 0.9816}}) However, NMinimize surrenders(!!!). Typing <<><><><><><><><><><><> Johannes Ludsteck Economics Department University of Regensburg Universitaetsstrasse 31 93053 Regensburg Reply-To: kuska@informatik.uni-leipzig.de ==== you guess right and if you hinder Mathematica to evaluate opt[] for symbolic arguments, with opt[s_?NumericQ] := Block[{x}, x /. Last[FindMinimum[x - 2.5(1 + Erf[x - s]), {x, 1, 3}]]] NMinimize[] works as expected. Jens > > Dear MathGroup Members, > > I want to minimize a function which returns the > minimizing value (arg min) of another function. > > For a simple example consider the following > function opt which returns the arg min of x-2.5(1+Erf[x-s]). > > opt[s_]:=Block[{x}, x/. Last[ > FindMinimum[x-2.5(1+Erf[x-s]), {x,1,3}]]] > > Now in a second step I want (again this is only > a simple example for illustrative purposes) to minimize > (opt[s]-2)^2 with respect to s. > > FindMininum has no problems with this. > > FindMinimum[(opt[s]-2)^2,{s,0.9,1.1}] > {3.18689*^-23, {s -> 0.9816}}) > > However, NMinimize surrenders(!!!). Typing > > < NMinimize[(opt[s]-2)^2,{s,0.9,1.1}] > only leads to the error message > > FindMinimum::fmnum: Objective function > 0.1 - 2.5 (1. +Erf[0.1 - 1. s]) is not real at {x} = {1.}. > > There is nothing wrong with minimand. It has exactly > one minimum in the Interval[{0.9,1.1}]. > > I guess the reason is that NMinimize calls opt[s] > not with a numerical value for s. This causes the > problem, since opt again calls FindMinimum. > Why? Can someone explain the failure and tell me > how to avoid this drawback? Wolfram Research boasts > that NMinimize can handle any function... > > I hope that nobody will recommend me to use FindMinimum > here instead. I know that the example here could of > course be solved by FindMinimum, but my real world > application can not. > > Johannes Ludsteck > > <><><><><><><><><><><><> > Johannes Ludsteck > Economics Department > University of Regensburg > Universitaetsstrasse 31 > 93053 Regensburg ==== Your Mathematica 4.2 is certainly not like the one most of us have: > This reveals us another integral which Mathematica 4.1 fails to > calculate > > In[1] := Integrate[Cos[Pi/4 - z]^2/z, {z, 1, Infinity}] // N > Out[1]= -0.0173083 > > Out[2]= 4.1 for Microsoft Windows (November 2, 2000) > > but Mathematica 4.2 handles correctly > > In[1] := Integrate[Cos[Pi/4 - z]^2/z, {z, 1, Infinity}] > Out[1] = Integrate::idiv: Integral of... does not converge on > {1, Infinity). > > Out[2]= 4.2 for Microsoft Windows (February 28, 2002) Well, actually with my 4.2 we get: In[1]:= Integrate[Cos[Pi/4 - z]^2/z, {z, 1, Infinity}] Out[1]= (1/4)*(Pi - 2*SinIntegral[2]) > > > Even simpler, > > In[1] := Integrate[Cos[z]^2/z, {z, 1, Infinity}] > Out[1] = -EulerGamma/2 - Log[2]/2 + (EulerGamma - CosIntegral[2] + > Log[2])/2 > > Out[2]= 4.1 for Microsoft Windows (November 2, 2000) > > which is wrong while Mathematica 4.2 works excellent > > In[1] := Integrate[Cos[z]^2/z, {z, 1, Infinity}] > Out[1] = Integrate::idiv: Integral of...does not converge on > {1, Infinity). > > Out[2]= 4.2 for Microsoft Windows (February 28, 2002) Not so fast: In[2]:= Integrate[Cos[z]^2/z, {z, 1, Infinity}] Out[2]= -(EulerGamma/2) - Log[2]/2 + (1/2)*(EulerGamma - CosIntegral[2] + Log[2]) I'd speculate that fixing these two integrals in the beta stage some more important ones, so the original way of doing things was restored in the released version. However, it gets even more interesting if we load: << Calculus`Limit` In[4]:= Integrate[Cos[Pi/4-z]^2/z,{z,1,Infinity}] ReplaceRepeated:: rrlim :Exiting after Interval[{-1,1}]/z + Interval[{0, 1}]/z scanned 65536 times. Integrate:: idiv :Integral of Cos[Pi]/4 - z^2/z does not converge on {1, Infinity]. Out[4]= Integrate[Cos[Pi/4 - z]^2/z, {z, 1, Infinity}] In[5]:= Integrate[Cos[z]^2/z, {z, 1, Infinity}] ReplaceRepeated::rrlim:Exiting after Interval[{0,1}] scanned 65536 times. Integrate::idiv:Integral of Cos[z]^2/z does not converge on {1, Infinity} Out[5]= Integrate[Cos[z]^2/z, {z, 1, Infinity}] In[6]:= Out[6]= 4.2 for Mac OS X (June 4, 2002) > > > Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ ==== Murray, I don't like TraditionalForm for Inline cells either. I just use MenuCellDefault Inline Format TypeStandard Form. Also, when I design my own style sheets I often define a bolder font for Inline cells so they stand out better. David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ soon as I press the Control-^ key combination, an Inline cell is created beginning with the x, and then when I type the exponent 2 everything in that Inline cell is now in Times, and the x is Italic. To change both characters to Courier is not so easy: it seems to require separately the entire Inline cell and selecting Courier does not change the exponent!) So to avoid this annoyance I normally must first type the desired expression in a separate Input cell, then copy the contents of that cell to the desired point in the Text cell. Any suggestions on a more efficient method for handling this? > In receiving notebooks from many different people I have noticed that > beginners often do not know how to use Text cells .... > > Some users may hesitate to use Text cells because they want to include a > mathematical expression in the comments.... > Just use an Inline cell within the text cell.... -- Murray Eisenberg murray@math.umass.edu Mathematics & Statistics Dept. Lederle Graduate Research Tower phone 413 549-1020 (H) University of Massachusetts 413 545-2859 (W) 710 North Pleasant Street Amherst, MA 01375 Reply-To: murray@math.umass.edu ==== That approach wouldn't work for me, since I often DO have to include within text cells some expressions in traditional mathematical notation and others in Mathematica's Standard Form. (The reason is that I'm often writing exposition as to how to express mathematical ideas and procedures in terms of Mathematica.) So probably the best way -- I'm not sure yet how to do it nicely -- would be a palette that more quickly allows me to change the format of a highlighted Inline cell to one or the other. Which reminds me of a related formatting matter. Often I need to include several paragraphs within a text cell, including displayed Inline cells on their own, separate lines. (No separate text cells would NOT meet my needs here.) The thing I usually do is to select the whole cell and from the Options Inspector sucessively select Formatting Options > Text Layout Options > ParagraphSpacing and then change the setting multiply from its default value 0 to 0.5. That provides just the right amount of inter-paragraph space. One of these days I'll figure out how to program a button on a palette to do that. General observation: It's stuff like this that makes Mathematica so much harder than a traditional tool, such as TeX/LaTeX, for typesetting mathematical exposition. Nothing beats a markup language for speed of entry. At least a sufficiently generous supply of formatting buttons would be the next best thing (far superior to having to burrow down through a nested menu in the Options Inspector). For example, I always keep open the palette FormattingTools.nb that allows changing the text face, size, or font (Courier, Times, Helvetica) at a click. (Not sure where I got the FormattingTools.nb palette from; it's not part of the standard Mathematica distribution. Perhaps it was from Publicon?) > Murray, > > I don't like TraditionalForm for Inline cells either. I just use > MenuCellDefault Inline Format TypeStandard Form.... -- Murray Eisenberg murray@math.umass.edu Mathematics & Statistics Dept. Lederle Graduate Research Tower phone 413 549-1020 (H) University of Massachusetts 413 545-2859 (W) 710 North Pleasant Street Amherst, MA 01375 ==== What is the TMJ style sheet? In any case, when people design new style sheets they would be well advised to keep Alt-7 as the key for Text style, and maintain the keys for other common styles also. Mathematica assigns keys in the order that the cell definitions appear in the style sheet. That means that new cell styles should be moved lower in the style sheet, even if that would not be their natural order. It would be nice if WRI allowed us to explicitly assign the keys for the styles. David Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.net Approved: Steven M. Christensen , Moderator ==== Here's a better version of DotPlot that is improved by Borut's and Allan's idea: Needs[Utilities`FilterOptions`]; DotPlot[data_?MatrixQ, opts___Rule] := Module[{colorfn, scale, dotsize, m, scalefn}, colorfn = ColorFunction/.{opts}/.{ColorFunction->Hue}; scale = ColorFunctionScaling/.{opts}/.{ColorFunctionScaling->True}; dotsize = DotSize/.{opts}/.{DotSize->0.01}; m = If[scale, {Min[#],Max[#]}& @ Last[Transpose[data]], {0,1}]; scalefn = (# - m[[1]])/(m[[2]]-m[[1]])& ; Show[ Graphics[ {PointSize[dotsize],colorfn[scalefn[#3]],Point[{#1,#2}]}&@@@data], FilterOptions[Graphics, opts] ] ] --Selwyn ==== Here's a refinement of my previous post. Everything is wrapped up in a function named DotPlot, which takes the following options: ColorFunction (default: Hue) ColorFunctionScaling (default: True) DotSize (default: 0.01) You can also provide other options such as Axes, Frame, Background, etc. Needs[Utilities`FilterOptions`]; DotPlot[data_?MatrixQ, opts___Rule] := Module[{colorfn, scale, dotsize}, colorfn = ColorFunction/.{opts}/.{ColorFunction->Hue}; scale = ColorFunctionScaling/.{opts}/.{ColorFunctionScaling->True}; dotsize = DotSize/.{opts}/.{DotSize->.01}; With[{m = If[scale,{Min[#],Max[#]}&@Transpose[data][[3]], {0,1}]}, Show[(Graphics[ {colorfn[#[[3]]],PointSize[dotsize],Point[{#[[1]],#[[2]]}]}]&)/@ (data /. {x_,y_,z_}-> {x,y,(z-m[[1]])/(m[[2]]-m[[1]])}), FilterOptions[Graphics, opts] ] ] ] Example: vals = Table[{Random[], Random[], .5 Random[]}, {100}]; DotPlot[vals, ColorFunction -> (Hue[#, 1, 1-#] & ), ColorFunctionScaling -> False, DotSize -> 0.02, Frame -> True] --- Selwyn Hollis ==== With verision 3.0 on a machine running Windows NT the following code produces rotated text with each of the charaters rotated as well degstr[th_]:=StringJoin[ToString[th], Degrees] str[th_,offset_]:={Point[{th,1}], Text[degstr[th],{th,1},offset,{Cos[th Degree],Sin[th Degree]}]} pic[offset_,plotrange_,angles_]:= Show[ Graphics[{PointSize[.015],str[#,offset]&/@angles},PlotRange->plotrange, AspectRati->.2,Frame[Rule]True,FrameTick->None, DefaultFon->{Courier,10}]]; pic[{0,0},{{-10,95},{-.1,2.1}},Range[0,90,10]]; With version 4.2 on a machine running Mac OS 10.2.1 the text is rotated but the characters in the text are not. That is the position of the chacters relative to each other changes but their orientation remains constant. I would like to have the character orientation change as it does when this code is run under Windows NT. Does anyone know if there is a setting somewhere that controls character orientation when text is rotated? ==== This should be on faq. Mac Mathematica lack this particular feature (This is not a bug, according to WRI). You can use HERSHEY package available from MathSource (ftp.mathsource.com). Example: s=.08; g=Plot[Sin[x],{x,0,Pi},DisplayFunction->Identity]; g1=Graphics[{AbsoluteThickness[.5], hText[x-axis values,{0,0}], Hue[0,0,1],hText[(,{-14,0}]}, AspectRatio->Automatic]; s1=s/2*Divide@@(Max@#-Min@#&/@ Transpose[Sequence@@#&/@Join[g1[[1,2]],g1[[1,4]]]/. Line->(Sequence@@#&)]); g2=Graphics[{AbsoluteThickness[.5], hText[y-axis values,{0,0},1,90 Degree], Hue[0,0,1],hText[(,{0,-14},1,90 Degree]}, AspectRatio->Automatic]; s2=s/2/Divide@@(Max@#-Min@#&/@ Transpose[Sequence@@#&/@Join[g2[[1,2]],g2[[1,4]]]/. Line->(Sequence@@#&)]); Show[Graphics@Rectangle[Scaled@{s,s},Scaled@{1,1},g], Graphics@Rectangle[Scaled@{.56-s1,0},Scaled@{.56+s1,s},g1], Graphics@Rectangle[Scaled@{0,.56-s2},Scaled@{s,.56+s2},g2], DisplayFunction->$DisplayFunction]; DH ==== Andrzej, observe: a = 77617; b = 33096; f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b) f Accuracy[f] Precision[f] Accuracy[10.^21] Precision[10.^21] 1.1805916207174113*^21 -5 16 -5 16 The two numbers supposedly have the same accuracy and precision, yet the first is in doubt by about 22 ORDERS OF MAGNITUDE -- never mind the digits! Mathematica computed this beast without giving any indication it was just noise -- without REALIZING it was noise. >>What more do you demand? I'm not demanding, objecting, or criticizing. I'm pointing out the problem. On the one hand, the poster is computing something that's simply not well-behaved. Unless he knows the coefficients with VERY high precision, he can't know even the magnitude of the result -- and that's not Mathematica's fault at all. On the other hand, Mathematica doesn't notice that precision is lost in the computation, and perhaps it should. You thought it DID notice, after all -- but it didn't. Bobby -----Original Message----- the same answer using the two different routes the original poster used. You on the other hand choose to assume that the posting shows that its author does not understand not just the mechanism of significance arithmetic used by Mathematica but also computations with inexact numbers in general. I do not think this is necessarily the correct assumption. I also don't see that Mathematica is doing anything wrong. After all, one can always check the accuracy of your answer: In[1]:= a = 77617.; b = 33096.; In[2]:= f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b) In[3]:= f Out[3]= -1.1805916207174113*^21 In[4]:= Accuracy[%] Out[4]= -5 which tells you that it can't be very reliable. What more do you demand? Andrzej Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ > Actually, we don't know whether SetAccuracy succeeds, because we > don't > know how inexact those numbers really are. If they are known to more > digits than shown in the original post, they should be entered with as > much precision as they deserve. If not, there's no trick or algorithm > that will give the result precision, because the value of f really is > in the noise. That is, we have no idea what the value of f should > be. > > Mathematica's failing is in returning a value without pointing out that > it has no precision. > > Bobby > > -----Original Message----- > Sent: Monday, September 30, 2002 11:59 AM > > Andrzej, Bobby, Peter > > It looks as if using SetAccuracy succeeds here because the inexact > numbers > that occur have finite binary representations. If we change them > slightly to > avoid this then we have to use Rationalize: > > 1) Using SetAccuracy > > Clear[a,b,f] > > > f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ > (2*b), > Infinity]; > > a=77617.1; > b=33096.1; > > a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; > > f > > > -1564032114908486835197494923989618867972540153873951942813115514949 > 3891236234 > > 52500771916869370459119776018798804630436149786919912931962574301029236 > 3 > 1246 > 75 > > / > 10867106143970760551000357827554793888198143135975649579607989867743572 > 8240 > 16 > 0653953612982932181371232436367739737604096 > > 2) Rewriting as fractions > > a=776171/10; > b=330961/10; > > f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) > > -(5954133808997234115690303589909929091649391296257/ > 41370125000000) > > 3) Using Rationalize > > Clear[a,b,f] > > > f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ > (2*b), > 0]; > > a=77617.1; > b=33096.1; > > a=Rationalize[a,0];b=Rationalize[b,0]; > > f > > -(5954133808997234115690303589909929091649391296257/ > 41370125000000) > > > I use Rationalize[. , 0] besause of results like > > Rationalize[3.1415959] > > 3.1416 > > Rationalize[3.1415959,0] > > 31415959/10000000 > -- > Allan > > --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 > > >> Well, first of of all, your using SetAccuracy and SetPrecision does >> nothing at all here, since they do not change the value of a or b. You >> should use a = SetAccuracy[a, Infinity] etc. But even then you won't >> get the same answer as when you use exact numbers because of the way >> you evaluate f. Here is the order of evaluation that will give you the >> same answer, and should explain what is going on: >> >> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* >> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; >> >> a = 77617.; >> >> >> b = 33096.; >> >> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; >> >> f >> >> 54767 >> -(-----) >> 66192 >> >> Andrzej Kozlowski >> Toyama International University >> JAPAN >> >> >> >> >>> Could someone explain what is going on here, please? >>> >>> In[1]:= >>> a = 77617.; b = 33096.; >>> >>> In[2]:= >>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; >>> >>> In[4]:= >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + >>> a/(2*b) >>> >>> In[5]:= >>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; >>> >>> In[6]:= >>> f >>> >>> Out[6]= >>> -1.1805916207174113*^21 >>> >>> In[7]:= >>> a = 77617; b = 33096; >>> >>> In[8]:= >>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>> (55/10)*b^8 + a/(2*b) >>> >>> In[9]:= >>> g >>> >>> Out[9]= >>> -(54767/66192) >>> >>> In[10]:= >>> N[%] >>> >>> Out[10]= >>> -0.8273960599468214 >>> >>> >>> PK >>> >>> >>> >> >> > > > > > > > > ==== Consider the total differential of f, with respect to the inexact numbers: Clear[a, b, x, y, f] f[a_, b_, x_, y_] := x*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + y*b^8 + a/(2*b) Simplify[Dt[f[a, b, x, y]] /. {Dt[a] -> da, Dt[b] -> db, Dt[x] -> dx, Dt[y] -> dy} /. {a -> 77617, b -> 33096, x -> 333.75, y -> 5.5}] -2.0400456966858126*^32*da + 4.784331242850472*^32*db + 1.3141745343712155*^27*dx + 1.4394747892125385*^36*dy f is sensitive to inaccuracy in the various numbers to widely varying degrees. Since the correct answer is small, we need a LOT of precision in the inputs to get there. If any of the inputs are merely machine precision numbers, we have NO precision in the result. The second and third terms are nearly the same magnitude with different signs. Even worse, the first term almost perfectly fills the gap: a = 77617; b = 33096; (33375/100)*b^6 438605750846393161930703831040 a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) -7917111779274712207494296632228773890 (55/10)*b^8 7917111340668961361101134701524942848 % + %% + %%% -2 Bobby Treat -----Original Message----- b = SetAccuracy[33096., Infinity]; In[4]:= f Out[4]= -(54767/66192) In[5]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity] Out[5]= 1180591620717411303424 Similarily: In[1]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 50]; a = SetAccuracy[77617., 100]; b = SetAccuracy[33096., 100]; In[4]:= f Out[4]= -0.8273960599468212641107299556`11.4133 In[5]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), 100]; Out[5]= 1.180591620717411303424`121.0721*^21 -PK ==== To Technical Support and the Mathematica User community, I'm writing to report what I consider to be a bug. First, I want to show a simplified example of the problem. Consider the following expression: expr=0.22 + x[0] + (3*(-0.16+ x[1]))/4 + (9*(0.546 + x[2]))/16; When simplified I expected to get some real number plus x[0]+3x[1]/4+9x[2]/16, but instead I get the following: Simplify[expr] 0.407125 + x[0] + 0.75 x[1] + 0.5625 x[2] As you can see, for some reason Mathematica converted the fractions 3/4 and 9/16 to real machine numbers. I consider this to be a bug. Now, for an example more representative of the situation that I've been coming across. expr12 = 1.001`17 + Sum[(x[i] - 1.001`17)/2^i, {i, 12}]; expr13 = 1.001`17 + Sum[(x[i] - 1.001`17)/2^i, {i, 13}]; expr55 = 1.001`17 + Sum[(x[i] - 1.001`17)/2^i, {i, 55}]; As you can see, I have replaced the real numbers by extended precision numbers. The simplified example above demonstrates that the problem exists when using machine numbers. Now, we'll see what happens when we use arbitrary precision numbers. First, let's simplify the expression with 12 terms. Simplify[expr12] (1.0010000000000 + 2048 x[1] + 1024 x[2] + 512 x[3] + 256 x[4] + 128 x[5] + 64 x[6] + 32 x[7] + 16 x[8] + 8 x[9] + 4 x[10] + 2 x[11] + x[12]) / 4096 As you can see, a sum with 12 terms upon simplification has coefficients which are still integers as they should be. However, increasing the number of terms to 13 yields Simplify[expr13] 0.0001221923828125 + 0.500000000000 x[1] + 0.250000000000 x[2] + 0.1250000000000 x[3] + 0.0625000000000 x[4] + 0.0312500000000 x[5] + 0.01562500000000 x[6] + 0.00781250000000 x[7] + 0.00390625000000 x[8] + 0.001953125000000 x[9] + 0.000976562500000 x[10] + 0.000488281250000 x[11] + 0.000244140625000 x[12] + 0.0001220703125000 x[13] Now, all of the coefficients are converted to real numbers, replicating the bug from the simplified example. Finally, let's see what happens when we have 55 terms. Rather than putting the resulting expression here, I will just leave it at the end of the post. The result though is somewhat surprising. Each of the coefficients of the x[i] are again real numbers, but now their precision is only 0! The proper result of course is the sum of some real number (with a precision close to 0 due to numerical cancellation) and an expression consisting of rational numbers multiplied by x[i]. The loss of precision of the coefficients of the x[i] is precisely what occurred to me. Of course, in this simple example, simply using Expand instead of Simplify produces the expected result, and is my workaround. I hope you agree with me that this is a bug, and one that Wolfram needs to correct. Carl Woll Physics Dept U of Washington Simplify[expr55] -16 -1 -1 -1 0. 10 + 0. x[1] + 0. x[2] + 0. 10 x[3] + 0. 10 x[4] + 0. 10 x[5] + -2 -2 -2 -3 -3 0. 10 x[6] + 0. 10 x[7] + 0. 10 x[8] + 0. 10 x[9] + 0. 10 x[10] + -3 -3 -4 -4 -4 0. 10 x[11] + 0. 10 x[12] + 0. 10 x[13] + 0. 10 x[14] + 0. 10 x[15] + -5 -5 -5 -6 -6 0. 10 x[16] + 0. 10 x[17] + 0. 10 x[18] + 0. 10 x[19] + 0. 10 x[20] + -6 -6 -7 -7 -7 0. 10 x[21] + 0. 10 x[22] + 0. 10 x[23] + 0. 10 x[24] + 0. 10 x[25] + -8 -8 -8 -9 -9 0. 10 x[26] + 0. 10 x[27] + 0. 10 x[28] + 0. 10 x[29] + 0. 10 x[30] + -9 -9 -10 -10 0. 10 x[31] + 0. 10 x[32] + 0. 10 x[33] + 0. 10 x[34] + -10 -11 -11 -11 0. 10 x[35] + 0. 10 x[36] + 0. 10 x[37] + 0. 10 x[38] + -12 -12 -12 -12 0. 10 x[39] + 0. 10 x[40] + 0. 10 x[41] + 0. 10 x[42] + -13 -13 -13 -14 0. 10 x[43] + 0. 10 x[44] + 0. 10 x[45] + 0. 10 x[46] + -14 -14 -15 -15 0. 10 x[47] + 0. 10 x[48] + 0. 10 x[49] + 0. 10 x[50] + -15 -15 -16 -16 - 16 0. 10 x[51] + 0. 10 x[52] + 0. 10 x[53] + 0. 10 x[54] + 0. 10 x[55] ==== Look at the FullForm of your expressions: FullForm[Sqrt[X + Y + Z]] Power[Plus[X, Y, Z], Rational[1, 2]] FullForm[f] Plus[Times[Rational[1, 4], C, Power[F, -2], Power[Plus[X, Y, Z], Rational[-1, 2]]], Times[B, Plus[A, Power[Plus[X, Y, Z], Rational[1, 2]]]]] FullForm[h] Plus[Times[B, Plus[A, Power[Plus[X, Y, Z], Rational[1, 2]]]], Times[C, Power[Plus[Times[4, Power[F, 2]], Power[Plus[X, Y, Z], Rational[1, 2]]], -1]]] The square root appears in two different forms! A solution is to replace both patterns: f = B*(A + Sqrt[X + Y + Z]) + C/(Sqrt[X + Y + Z]/4*F^2); f /. {Sqrt[X + Y + Z] -> Q, 1/Sqrt[X + Y + Z] -> 1/Q} (4*C)/(F^2*Q) + B*(A + Q) This is very annoying, of course, and it may not take care of every case. Looking at the FullForm should help, when it fails. Bobby -----Original Message----- Print[(f/.Sqrt[X+Y+Z]->Q)//InputForm]; B*(A + Q) + (4*C)/(F^2*Sqrt[X + Y + Z]) (* fails *) g=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]*4*F^2); Print[(g/.Sqrt[X+Y+Z]->Q)//InputForm]; B*(A + Q) + C/(4*F^2*Sqrt[X + Y + Z]) (* fails *) h=B*(A+Sqrt[X+Y+Z])+C/(Sqrt[X+Y+Z]+4*F^2); Print[(h/.Sqrt[X+Y+Z]->Q)//InputForm]; B*(A + Q) + C/(4*F^2 + Q) (* works *) ==== >>I should add that the solution is over natural numbers.. this will probably make a big difference.. Yes, that probably removes the nearly from nearly impossible. I'd be curious to see the actual problem. Bobby -----Original Message----- eqns.. > preferably ones that will minimize the computation time for the other 14). > These equation are not linearly related.. the highest degree in any one eqn > is degree 4 i believe.. and there are some cross terms in the equations but > not every equation depends on every variable.. (some are actually rather > simple eqns). Any ideas on how to get started with this using > mathematica (any ideas for algorithms).. > > Anything will be helpful.. > > I can be reached at ngupta2@seas.upenn.edu > > Many thanks, > > Nachi > > > > ==== Group, in the help browser (V4.2) there is the following about Return[]: (* Return [expr] returns the value expr, existing all procedures and loops in a function *). Besides, from 2.5.9 Loops and control structures (also in Help browser) we can conclude that While, Do, Module, With, etc. must be understood as loops and/or control structures. But let us have a look at these two function definitions: yes[a_]:=With[{Ever=True},While[Ever,If[a==7,Return[Terminate] ];Print[loop]]], In (*yes[]*), Return[] exits the function breaking all loops and control structures, and yields Terminate as expected. not[a_]:=With[{Ever=True},While[Ever,While[Ever,If[a==7,Do[Return[Terminat e];Print[foo] ]];Print[loop]]]], However, in (*not[]*), Return[] only breaks the Do[] construct and keeps looping inside While[]. To my understanding, this contradicts that Return [expr] returns the value expr, existing all procedures and loops in a function as stated executing Return[] anywhere at any deep of nested looping constructs inside the function. I would appreciate any feedback. Emilio Martin-Serrano ==== You are of course right, I forgot that Mathematica does not try to keep precision or accuracy of machine arithmetic computations. But I think the actual precision you set need not be higher than machine precision, it just has to be set explicitely (is that right?). For example: In[1]:= Clear[f,a,b,k] In[2]:= k = $MachinePrecision; In[3]:= f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ (2*b),k]; In[4]:= a=77617.;b=33096.; In[5]:= a=SetAccuracy[a,k]; In[6]:= b=SetAccuracy[b,k]; In[7]:= f Out[7]= !((-5.51716400890319`-2.8311*^19)) In[8]:= Accuracy[f] Accuracy::mnprec: Value -23 would be inconsistent with $MinPrecision; bounding by $MinPrecision instead. Out[8]= -20 Andrzej > Andrzej > > Yes, like you I took the original question to be about how to get the > result > of using the naive rational versions in place of the inexact numbers. > Bobby raises the question of how we might know accuracy of the result. > > You answer this with >> a = 77617.; b = 33096.; >> >> In[2]:= >> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >> 5.5*b^8 + a/(2*b) >> >> In[3]:= >> f >> >> Out[3]= >> -1.1805916207174113*^21 >> >> In[4]:= >> Accuracy[%] >> >> Out[4]= >> -5 > > However this computation is done in machine arithmetic, which means > that > Mathematica keeps no check on the accuracy and precision of the > computation, > and Mathematica gives the default accuracy value without any real > signifcance: > > $MachinePrecision - Log[10,Abs[f]]//Round > > -5 > > To get meaningful accuracy and precision values we need to force the > computation to be in bignums (bigfloat, arbitrary precision) > arithmetic. > Mathematica then keeps track of the accuracy and precision that it can > guarantee. > > Clear[f,a,b,k] > k=50; > > f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ > (2*b),k]; > a=77617.;b=33096.; > a=SetAccuracy[a,k]; > b=SetAccuracy[b,k]; > f > > -0.82739605995 > > Accuracy[f] > > 11 > > Precision[f] > > 11 > > Which tells us that the error in the computed value of f is not more > than > > 10^-11 > > > The above results are rounded. More accurately we get > > Accuracy[f, Round->False] > > 11.4956 > > Precision[f, Round->False] > > 11.4133 > > There are several related issues here: > - is the precision (accuracy) of rhe input known? > - how does Mathematica interpret what one gives it? > - how does Mathematica go about the calculation? > > -- > Allan > > --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 > > >> It seems clear to me that what Allan and what you mean by succeeds >> here refer to quite different things and your objection is therefore >> beside the point. There are obviously two ways in which one can >> interpret the original posting. The first interpretation, which Allan >> and myself adopted, was that the question concerned purely the >> computational mechanism of Mathematica. Or, to put it in other words, >> it was why are the results of these two computations not the same?. >> In this sense success means no more than making Mathematica return >> the same answer using the two different routes the original poster >> used. >> You on the other hand choose to assume that the posting shows that its >> author does not understand not just the mechanism of significance >> arithmetic used by Mathematica but also computations with inexact >> numbers in general. I do not think this is necessarily the correct >> assumption. I also don't see that Mathematica is doing anything wrong. >> After all, one can always check the accuracy of your answer: >> >> In[1]:= >> a = 77617.; b = 33096.; >> >> In[2]:= >> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >> 5.5*b^8 + a/(2*b) >> >> In[3]:= >> f >> >> Out[3]= >> -1.1805916207174113*^21 >> >> In[4]:= >> Accuracy[%] >> >> Out[4]= >> -5 >> >> which tells you that it can't be very reliable. What more do you >> demand? >> >> Andrzej >> >> >> Andrzej Kozlowski >> Yokohama, Japan >> http://www.mimuw.edu.pl/~akoz/ >> http://platon.c.u-tokyo.ac.jp/andrzej/ >> >> >>> Actually, we don't know whether SetAccuracy succeeds, because we >>> don't >>> know how inexact those numbers really are. If they are known to more >>> digits than shown in the original post, they should be entered with >>> as >>> much precision as they deserve. If not, there's no trick or >>> algorithm >>> that will give the result precision, because the value of f really is >>> in the noise. That is, we have no idea what the value of f should >>> be. >>> >>> Mathematica's failing is in returning a value without pointing out >>> that >>> it has no precision. >>> >>> Bobby >>> >>> -----Original Message----- >>> Sent: Monday, September 30, 2002 11:59 AM >>> >>> Andrzej, Bobby, Peter >>> >>> It looks as if using SetAccuracy succeeds here because the inexact >>> numbers >>> that occur have finite binary representations. If we change them >>> slightly to >>> avoid this then we have to use Rationalize: >>> >>> 1) Using SetAccuracy >>> >>> Clear[a,b,f] >>> >>> >>> f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>> (2*b), >>> Infinity]; >>> >>> a=77617.1; >>> b=33096.1; >>> >>> a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; >>> >>> f >>> >>> >>> -1564032114908486835197494923989618867972540153873951942813115514949 >>> 3891236234 >>> >>> 525007719168693704591197760187988046304361497869199129319625743010292 >>> 36 >>> 3 >>> 1246 >>> 75 >>> >>> / >>> 108671061439707605510003578275547938881981431359756495796079898677435 >>> 72 >>> 8240 >>> 16 >>> 0653953612982932181371232436367739737604096 >>> >>> 2) Rewriting as fractions >>> >>> a=776171/10; >>> b=330961/10; >>> >>> f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) >>> >>> -(5954133808997234115690303589909929091649391296257/ >>> 41370125000000) >>> >>> 3) Using Rationalize >>> >>> Clear[a,b,f] >>> >>> >>> f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>> (2*b), >>> 0]; >>> >>> a=77617.1; >>> b=33096.1; >>> >>> a=Rationalize[a,0];b=Rationalize[b,0]; >>> >>> f >>> >>> -(5954133808997234115690303589909929091649391296257/ >>> 41370125000000) >>> >>> >>> I use Rationalize[. , 0] besause of results like >>> >>> Rationalize[3.1415959] >>> >>> 3.1416 >>> >>> Rationalize[3.1415959,0] >>> >>> 31415959/10000000 >>> -- >>> Allan >>> >>> --------------------- >>> Allan Hayes >>> Mathematica Training and Consulting >>> Leicester UK >>> www.haystack.demon.co.uk >>> hay@haystack.demon.co.uk >>> Voice: +44 (0)116 271 4198 >>> >>> >>>> Well, first of of all, your using SetAccuracy and SetPrecision does >>>> nothing at all here, since they do not change the value of a or b. >>>> You >>>> should use a = SetAccuracy[a, Infinity] etc. But even then you won't >>>> get the same answer as when you use exact numbers because of the >>>> way >>>> you evaluate f. Here is the order of evaluation that will give you >>>> the >>>> same answer, and should explain what is going on: >>>> >>>> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* >>>> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; >>>> >>>> a = 77617.; >>>> >>>> >>>> b = 33096.; >>>> >>>> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; >>>> >>>> f >>>> >>>> 54767 >>>> -(-----) >>>> 66192 >>>> >>>> Andrzej Kozlowski >>>> Toyama International University >>>> JAPAN >>>> >>>> >>>> >>>> >>>>> Could someone explain what is going on here, please? >>>>> >>>>> In[1]:= >>>>> a = 77617.; b = 33096.; >>>>> >>>>> In[2]:= >>>>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >>>>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; >>>>> >>>>> In[4]:= >>>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + >>>>> a/(2*b) >>>>> >>>>> In[5]:= >>>>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; >>>>> >>>>> In[6]:= >>>>> f >>>>> >>>>> Out[6]= >>>>> -1.1805916207174113*^21 >>>>> >>>>> In[7]:= >>>>> a = 77617; b = 33096; >>>>> >>>>> In[8]:= >>>>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>>> (55/10)*b^8 + a/(2*b) >>>>> >>>>> In[9]:= >>>>> g >>>>> >>>>> Out[9]= >>>>> -(54767/66192) >>>>> >>>>> In[10]:= >>>>> N[%] >>>>> >>>>> Out[10]= >>>>> -0.8273960599468214 >>>>> >>>>> >>>>> PK >>>>> >>>>> >>>>> >>>> >>>> >>> >>> >>> >>> >>> >>> >>> >>> >> >> >> > > > > > > > > > > Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ ==== So? It's just as it should be, isn't it? Andrzej Andrzej Kozlowski Toyama International University JAPAN http://sigma.tuins.ac.jp/~andrzej/ > Go one step further: > > Clear[f, a, b, k] > k = $MachinePrecision; > f = SetAccuracy[333.75*b^6 + > a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + > 5.5*b^8 + a/(2*b), k]; > a = 77617.; b = 33096.; > a = SetAccuracy[a, k]; > b = SetAccuracy[b, k]; > f > Accuracy[f] > Precision[f] > > -5.517164009`0*^19 > > Accuracy::mnprec:Value -23 would be inconsistent with $MinPrecision; > bounding by $MinPrecision instead. > > -20 > > 0 > > Accuracy::mnprec:Value !(-23) would be inconsistent with > $MinPrecision; > bounding by $MinPrecision instead. > > See that? NO precision. > > Bobby > > -----Original Message----- > Sent: Tuesday, October 01, 2002 5:53 PM > Cc: drbob@bigfoot.com > > You are of course right, I forgot that Mathematica does not try to keep > > precision or accuracy of machine arithmetic computations. But I think > the actual precision you set need not be higher than machine precision, > > it just has to be set explicitely (is that right?). For example: > > In[1]:= > Clear[f,a,b,k] > > In[2]:= > k = $MachinePrecision; > > In[3]:= > f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ > (2*b),k]; > > In[4]:= > a=77617.;b=33096.; > > In[5]:= > a=SetAccuracy[a,k]; > > In[6]:= > b=SetAccuracy[b,k]; > > > In[7]:= > f > > Out[7]= > !((-5.51716400890319`-2.8311*^19)) > > In[8]:= > Accuracy[f] > > Accuracy::mnprec: Value -23 would be inconsistent with $MinPrecision; > bounding by $MinPrecision instead. > > Out[8]= > -20 > > > Andrzej > > > >> Andrzej >> >> Yes, like you I took the original question to be about how to get the > >> result >> of using the naive rational versions in place of the inexact numbers. >> Bobby raises the question of how we might know accuracy of the result. >> >> You answer this with >>> a = 77617.; b = 33096.; >>> >>> In[2]:= >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>> 5.5*b^8 + a/(2*b) >>> >>> In[3]:= >>> f >>> >>> Out[3]= >>> -1.1805916207174113*^21 >>> >>> In[4]:= >>> Accuracy[%] >>> >>> Out[4]= >>> -5 >> >> However this computation is done in machine arithmetic, which means >> that >> Mathematica keeps no check on the accuracy and precision of the >> computation, >> and Mathematica gives the default accuracy value without any real >> signifcance: >> >> $MachinePrecision - Log[10,Abs[f]]//Round >> >> -5 >> >> To get meaningful accuracy and precision values we need to force the >> computation to be in bignums (bigfloat, arbitrary precision) >> arithmetic. >> Mathematica then keeps track of the accuracy and precision that it can >> guarantee. >> >> Clear[f,a,b,k] >> k=50; >> >> f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ >> (2*b),k]; >> a=77617.;b=33096.; >> a=SetAccuracy[a,k]; >> b=SetAccuracy[b,k]; >> f >> >> -0.82739605995 >> >> Accuracy[f] >> >> 11 >> >> Precision[f] >> >> 11 >> >> Which tells us that the error in the computed value of f is not more > >> than >> >> 10^-11 >> >> >> The above results are rounded. More accurately we get >> >> Accuracy[f, Round->False] >> >> 11.4956 >> >> Precision[f, Round->False] >> >> 11.4133 >> >> There are several related issues here: >> - is the precision (accuracy) of rhe input known? >> - how does Mathematica interpret what one gives it? >> - how does Mathematica go about the calculation? >> >> -- >> Allan >> >> --------------------- >> Allan Hayes >> Mathematica Training and Consulting >> Leicester UK >> www.haystack.demon.co.uk >> hay@haystack.demon.co.uk >> Voice: +44 (0)116 271 4198 >> >> >>> It seems clear to me that what Allan and what you mean by succeeds >>> here refer to quite different things and your objection is therefore >>> beside the point. There are obviously two ways in which one can >>> interpret the original posting. The first interpretation, which Allan >>> and myself adopted, was that the question concerned purely the >>> computational mechanism of Mathematica. Or, to put it in other words, >>> it was why are the results of these two computations not the same?. >>> In this sense success means no more than making Mathematica return >>> the same answer using the two different routes the original poster >>> used. >>> You on the other hand choose to assume that the posting shows that > its >>> author does not understand not just the mechanism of significance >>> arithmetic used by Mathematica but also computations with inexact >>> numbers in general. I do not think this is necessarily the correct >>> assumption. I also don't see that Mathematica is doing anything > wrong. >>> After all, one can always check the accuracy of your answer: >>> >>> In[1]:= >>> a = 77617.; b = 33096.; >>> >>> In[2]:= >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>> 5.5*b^8 + a/(2*b) >>> >>> In[3]:= >>> f >>> >>> Out[3]= >>> -1.1805916207174113*^21 >>> >>> In[4]:= >>> Accuracy[%] >>> >>> Out[4]= >>> -5 >>> >>> which tells you that it can't be very reliable. What more do you >>> demand? >>> >>> Andrzej >>> >>> >>> Andrzej Kozlowski >>> Yokohama, Japan >>> http://www.mimuw.edu.pl/~akoz/ >>> http://platon.c.u-tokyo.ac.jp/andrzej/ >>> >>> >>>> Actually, we don't know whether SetAccuracy succeeds, because we >>>> don't >>>> know how inexact those numbers really are. If they are known to > more >>>> digits than shown in the original post, they should be entered with > >>>> as >>>> much precision as they deserve. If not, there's no trick or >>>> algorithm >>>> that will give the result precision, because the value of f really > is >>>> in the noise. That is, we have no idea what the value of f should >>>> be. >>>> >>>> Mathematica's failing is in returning a value without pointing out >>>> that >>>> it has no precision. >>>> >>>> Bobby >>>> >>>> -----Original Message----- >>>> Sent: Monday, September 30, 2002 11:59 AM >>>> >>>> Andrzej, Bobby, Peter >>>> >>>> It looks as if using SetAccuracy succeeds here because the inexact >>>> numbers >>>> that occur have finite binary representations. If we change them >>>> slightly to >>>> avoid this then we have to use Rationalize: >>>> >>>> 1) Using SetAccuracy >>>> >>>> Clear[a,b,f] >>>> >>>> >>>> f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>>> (2*b), >>>> Infinity]; >>>> >>>> a=77617.1; >>>> b=33096.1; >>>> >>>> a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; >>>> >>>> f >>>> >>>> >>>> -1564032114908486835197494923989618867972540153873951942813115514949 >>>> 3891236234 >>>> >>>> > 525007719168693704591197760187988046304361497869199129319625743010292 >>>> 36 >>>> 3 >>>> 1246 >>>> 75 >>>> >>>> / >>>> > 108671061439707605510003578275547938881981431359756495796079898677435 >>>> 72 >>>> 8240 >>>> 16 >>>> 0653953612982932181371232436367739737604096 >>>> >>>> 2) Rewriting as fractions >>>> >>>> a=776171/10; >>>> b=330961/10; >>>> >>>> f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) >>>> >>>> -(5954133808997234115690303589909929091649391296257/ >>>> 41370125000000) >>>> >>>> 3) Using Rationalize >>>> >>>> Clear[a,b,f] >>>> >>>> >>>> f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>>> (2*b), >>>> 0]; >>>> >>>> a=77617.1; >>>> b=33096.1; >>>> >>>> a=Rationalize[a,0];b=Rationalize[b,0]; >>>> >>>> f >>>> >>>> -(5954133808997234115690303589909929091649391296257/ >>>> 41370125000000) >>>> >>>> >>>> I use Rationalize[. , 0] besause of results like >>>> >>>> Rationalize[3.1415959] >>>> >>>> 3.1416 >>>> >>>> Rationalize[3.1415959,0] >>>> >>>> 31415959/10000000 >>>> -- >>>> Allan >>>> >>>> --------------------- >>>> Allan Hayes >>>> Mathematica Training and Consulting >>>> Leicester UK >>>> www.haystack.demon.co.uk >>>> hay@haystack.demon.co.uk >>>> Voice: +44 (0)116 271 4198 >>>> >>>> >>>>> Well, first of of all, your using SetAccuracy and SetPrecision does >>>>> nothing at all here, since they do not change the value of a or b. > >>>>> You >>>>> should use a = SetAccuracy[a, Infinity] etc. But even then you > won't >>>>> get the same answer as when you use exact numbers because of the >>>>> way >>>>> you evaluate f. Here is the order of evaluation that will give you > >>>>> the >>>>> same answer, and should explain what is going on: >>>>> >>>>> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* >>>>> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; >>>>> >>>>> a = 77617.; >>>>> >>>>> >>>>> b = 33096.; >>>>> >>>>> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; >>>>> >>>>> f >>>>> >>>>> 54767 >>>>> -(-----) >>>>> 66192 >>>>> >>>>> Andrzej Kozlowski >>>>> Toyama International University >>>>> JAPAN >>>>> >>>>> >>>>> >>>>> >>>>>> Could someone explain what is going on here, please? >>>>>> >>>>>> In[1]:= >>>>>> a = 77617.; b = 33096.; >>>>>> >>>>>> In[2]:= >>>>>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >>>>>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; >>>>>> >>>>>> In[4]:= >>>>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + >>>>>> a/(2*b) >>>>>> >>>>>> In[5]:= >>>>>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; >>>>>> >>>>>> In[6]:= >>>>>> f >>>>>> >>>>>> Out[6]= >>>>>> -1.1805916207174113*^21 >>>>>> >>>>>> In[7]:= >>>>>> a = 77617; b = 33096; >>>>>> >>>>>> In[8]:= >>>>>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>>>> (55/10)*b^8 + a/(2*b) >>>>>> >>>>>> In[9]:= >>>>>> g >>>>>> >>>>>> Out[9]= >>>>>> -(54767/66192) >>>>>> >>>>>> In[10]:= >>>>>> N[%] >>>>>> >>>>>> Out[10]= >>>>>> -0.8273960599468214 >>>>>> >>>>>> >>>>>> PK >>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>> >>> >>> >> >> >> >> >> >> >> >> >> >> > Andrzej Kozlowski > Yokohama, Japan > http://www.mimuw.edu.pl/~akoz/ > http://platon.c.u-tokyo.ac.jp/andrzej/ > > > > ==== I think I'd prefer that Mathematica gave me a clue -- without being asked explicitly in JUST the right way that only you know and had forgotten -- that there's a precision problem. Oddly enough, when you DO ask it nicely, you get error messages, but if you're not aware there's a problem, it lets you go on your merry way, working with noise. Bobby -----Original Message----- > > Clear[f, a, b, k] > k = $MachinePrecision; > f = SetAccuracy[333.75*b^6 + > a^2*(11*a^2*b^2 - b^6 - > 121*b^4 - 2) + > 5.5*b^8 + a/(2*b), k]; > a = 77617.; b = 33096.; > a = SetAccuracy[a, k]; > b = SetAccuracy[b, k]; > f > Accuracy[f] > Precision[f] > > -5.517164009`0*^19 > > Accuracy::mnprec:Value -23 would be inconsistent with $MinPrecision; > bounding by $MinPrecision instead. > > -20 > > 0 > > Accuracy::mnprec:Value !(-23) would be inconsistent with > $MinPrecision; > bounding by $MinPrecision instead. > > See that? NO precision. > > Bobby > > -----Original Message----- > Sent: Tuesday, October 01, 2002 5:53 PM > Cc: drbob@bigfoot.com > > You are of course right, I forgot that Mathematica does not try to keep > > precision or accuracy of machine arithmetic computations. But I think > the actual precision you set need not be higher than machine precision, > > it just has to be set explicitely (is that right?). For example: > > In[1]:= > Clear[f,a,b,k] > > In[2]:= > k = $MachinePrecision; > > In[3]:= > f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ > (2*b),k]; > > In[4]:= > a=77617.;b=33096.; > > In[5]:= > a=SetAccuracy[a,k]; > > In[6]:= > b=SetAccuracy[b,k]; > > > In[7]:= > f > > Out[7]= > !((-5.51716400890319`-2.8311*^19)) > > In[8]:= > Accuracy[f] > > Accuracy::mnprec: Value -23 would be inconsistent with $MinPrecision; > bounding by $MinPrecision instead. > > Out[8]= > -20 > > > Andrzej > > > >> Andrzej >> >> Yes, like you I took the original question to be about how to get the > >> result >> of using the naive rational versions in place of the inexact numbers. >> Bobby raises the question of how we might know accuracy of the result. >> >> You answer this with >>> a = 77617.; b = 33096.; >>> >>> In[2]:= >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>> 5.5*b^8 + a/(2*b) >>> >>> In[3]:= >>> f >>> >>> Out[3]= >>> -1.1805916207174113*^21 >>> >>> In[4]:= >>> Accuracy[%] >>> >>> Out[4]= >>> -5 >> >> However this computation is done in machine arithmetic, which means >> that >> Mathematica keeps no check on the accuracy and precision of the >> computation, >> and Mathematica gives the default accuracy value without any real >> signifcance: >> >> $MachinePrecision - Log[10,Abs[f]]//Round >> >> -5 >> >> To get meaningful accuracy and precision values we need to force the >> computation to be in bignums (bigfloat, arbitrary precision) >> arithmetic. >> Mathematica then keeps track of the accuracy and precision that it can >> guarantee. >> >> Clear[f,a,b,k] >> k=50; >> >> f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ >> (2*b),k]; >> a=77617.;b=33096.; >> a=SetAccuracy[a,k]; >> b=SetAccuracy[b,k]; >> f >> >> -0.82739605995 >> >> Accuracy[f] >> >> 11 >> >> Precision[f] >> >> 11 >> >> Which tells us that the error in the computed value of f is not more > >> than >> >> 10^-11 >> >> >> The above results are rounded. More accurately we get >> >> Accuracy[f, Round->False] >> >> 11.4956 >> >> Precision[f, Round->False] >> >> 11.4133 >> >> There are several related issues here: >> - is the precision (accuracy) of rhe input known? >> - how does Mathematica interpret what one gives it? >> - how does Mathematica go about the calculation? >> >> -- >> Allan >> >> --------------------- >> Allan Hayes >> Mathematica Training and Consulting >> Leicester UK >> www.haystack.demon.co.uk >> hay@haystack.demon.co.uk >> Voice: +44 (0)116 271 4198 >> >> >>> It seems clear to me that what Allan and what you mean by succeeds >>> here refer to quite different things and your objection is therefore >>> beside the point. There are obviously two ways in which one can >>> interpret the original posting. The first interpretation, which Allan >>> and myself adopted, was that the question concerned purely the >>> computational mechanism of Mathematica. Or, to put it in other words, >>> it was why are the results of these two computations not the same?. >>> In this sense success means no more than making Mathematica return >>> the same answer using the two different routes the original poster >>> used. >>> You on the other hand choose to assume that the posting shows that > its >>> author does not understand not just the mechanism of significance >>> arithmetic used by Mathematica but also computations with inexact >>> numbers in general. I do not think this is necessarily the correct >>> assumption. I also don't see that Mathematica is doing anything > wrong. >>> After all, one can always check the accuracy of your answer: >>> >>> In[1]:= >>> a = 77617.; b = 33096.; >>> >>> In[2]:= >>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>> 5.5*b^8 + a/(2*b) >>> >>> In[3]:= >>> f >>> >>> Out[3]= >>> -1.1805916207174113*^21 >>> >>> In[4]:= >>> Accuracy[%] >>> >>> Out[4]= >>> -5 >>> >>> which tells you that it can't be very reliable. What more do you >>> demand? >>> >>> Andrzej >>> >>> >>> Andrzej Kozlowski >>> Yokohama, Japan >>> http://www.mimuw.edu.pl/~akoz/ >>> http://platon.c.u-tokyo.ac.jp/andrzej/ >>> >>> >>>> Actually, we don't know whether SetAccuracy succeeds, because we >>>> don't >>>> know how inexact those numbers really are. If they are known to > more >>>> digits than shown in the original post, they should be entered with > >>>> as >>>> much precision as they deserve. If not, there's no trick or >>>> algorithm >>>> that will give the result precision, because the value of f really > is >>>> in the noise. That is, we have no idea what the value of f should >>>> be. >>>> >>>> Mathematica's failing is in returning a value without pointing out >>>> that >>>> it has no precision. >>>> >>>> Bobby >>>> >>>> -----Original Message----- >>>> Sent: Monday, September 30, 2002 11:59 AM >>>> >>>> Andrzej, Bobby, Peter >>>> >>>> It looks as if using SetAccuracy succeeds here because the inexact >>>> numbers >>>> that occur have finite binary representations. If we change them >>>> slightly to >>>> avoid this then we have to use Rationalize: >>>> >>>> 1) Using SetAccuracy >>>> >>>> Clear[a,b,f] >>>> >>>> >>>> f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>>> (2*b), >>>> Infinity]; >>>> >>>> a=77617.1; >>>> b=33096.1; >>>> >>>> a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; >>>> >>>> f >>>> >>>> >>>> -1564032114908486835197494923989618867972540153873951942813115514949 >>>> 3891236234 >>>> >>>> > 525007719168693704591197760187988046304361497869199129319625743010292 >>>> 36 >>>> 3 >>>> 1246 >>>> 75 >>>> >>>> / >>>> > 108671061439707605510003578275547938881981431359756495796079898677435 >>>> 72 >>>> 8240 >>>> 16 >>>> 0653953612982932181371232436367739737604096 >>>> >>>> 2) Rewriting as fractions >>>> >>>> a=776171/10; >>>> b=330961/10; >>>> >>>> f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) >>>> >>>> -(5954133808997234115690303589909929091649391296257/ >>>> 41370125000000) >>>> >>>> 3) Using Rationalize >>>> >>>> Clear[a,b,f] >>>> >>>> >>>> f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>>> (2*b), >>>> 0]; >>>> >>>> a=77617.1; >>>> b=33096.1; >>>> >>>> a=Rationalize[a,0];b=Rationalize[b,0]; >>>> >>>> f >>>> >>>> -(5954133808997234115690303589909929091649391296257/ >>>> 41370125000000) >>>> >>>> >>>> I use Rationalize[. , 0] besause of results like >>>> >>>> Rationalize[3.1415959] >>>> >>>> 3.1416 >>>> >>>> Rationalize[3.1415959,0] >>>> >>>> 31415959/10000000 >>>> -- >>>> Allan >>>> >>>> --------------------- >>>> Allan Hayes >>>> Mathematica Training and Consulting >>>> Leicester UK >>>> www.haystack.demon.co.uk >>>> hay@haystack.demon.co.uk >>>> Voice: +44 (0)116 271 4198 >>>> >>>> >>>>> Well, first of of all, your using SetAccuracy and SetPrecision does >>>>> nothing at all here, since they do not change the value of a or b. > >>>>> You >>>>> should use a = SetAccuracy[a, Infinity] etc. But even then you > won't >>>>> get the same answer as when you use exact numbers because of the >>>>> way >>>>> you evaluate f. Here is the order of evaluation that will give you > >>>>> the >>>>> same answer, and should explain what is going on: >>>>> >>>>> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* >>>>> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; >>>>> >>>>> a = 77617.; >>>>> >>>>> >>>>> b = 33096.; >>>>> >>>>> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; >>>>> >>>>> f >>>>> >>>>> 54767 >>>>> -(-----) >>>>> 66192 >>>>> >>>>> Andrzej Kozlowski >>>>> Toyama International University >>>>> JAPAN >>>>> >>>>> >>>>> >>>>> >>>>>> Could someone explain what is going on here, please? >>>>>> >>>>>> In[1]:= >>>>>> a = 77617.; b = 33096.; >>>>>> >>>>>> In[2]:= >>>>>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >>>>>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; >>>>>> >>>>>> In[4]:= >>>>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + >>>>>> a/(2*b) >>>>>> >>>>>> In[5]:= >>>>>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; >>>>>> >>>>>> In[6]:= >>>>>> f >>>>>> >>>>>> Out[6]= >>>>>> -1.1805916207174113*^21 >>>>>> >>>>>> In[7]:= >>>>>> a = 77617; b = 33096; >>>>>> >>>>>> In[8]:= >>>>>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>>>> (55/10)*b^8 + a/(2*b) >>>>>> >>>>>> In[9]:= >>>>>> g >>>>>> >>>>>> Out[9]= >>>>>> -(54767/66192) >>>>>> >>>>>> In[10]:= >>>>>> N[%] >>>>>> >>>>>> Out[10]= >>>>>> -0.8273960599468214 >>>>>> >>>>>> >>>>>> PK >>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>>> >>> >>> >>> >> >> >> >> >> >> >> >> >> >> > Andrzej Kozlowski > Yokohama, Japan > http://www.mimuw.edu.pl/~akoz/ > http://platon.c.u-tokyo.ac.jp/andrzej/ > > > > ==== I've been looking over the file and directory manipulation functions in Mathematica 4.1, and I'm not finding good examples of how to test for the existence of a file or directory before I attempt to create, access, or delete it. Are there any good examples of this somewhere? STH ==== FileNames[] does what you want. Steve Luttrell > I've been looking over the file and directory manipulation functions in Mathematica > 4.1, and I'm not finding good examples of how to test for the existence of > a file or directory before I attempt to create, access, or delete it. Are > there any good examples of this somewhere? > > STH > Reply-To: kuska@informatik.uni-leipzig.de ==== FileNames[] returns an empty list if there are no files. So FileNames[blub.txt] will return {} if the file does not exist. Jens > > I've been looking over the file and directory manipulation functions in Mathematica > 4.1, and I'm not finding good examples of how to test for the existence of > a file or directory before I attempt to create, access, or delete it. Are > there any good examples of this somewhere? > > STH ==== If you don't wont error messages all you need to do is to set MinPrecision low enough: In[1]:= $MinPrecision = -3; In[2]:= k = $MachinePrecision; In[3]:= f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 + a/(2*b), k]; In[4]:= a = 77617.; b = 33096.; In[5]:= a = SetAccuracy[a, k]; In[6]:= b = SetAccuracy[b, k]; In[7]:= f Out[7]= -5.517164009`-2.8311*^19 In[8]:= Accuracy[f] Out[8]= -23 In[9]:= Precision[f] Out[9]= -3 In any case the answer you get is meaningless. > I think I'd prefer that Mathematica gave me a clue -- without being > asked explicitly in JUST the right way that only you know and had > forgotten -- that there's a precision problem. > > Oddly enough, when you DO ask it nicely, you get error messages, but if > you're not aware there's a problem, it lets you go on your merry way, > working with noise. > > Bobby > > -----Original Message----- > Sent: Tuesday, October 01, 2002 8:05 PM > Cc: 'Allan Hayes'; mathgroup@smc.vnet.net > > So? It's just as it should be, isn't it? > > Andrzej > > Andrzej Kozlowski > Toyama International University > JAPAN > http://sigma.tuins.ac.jp/~andrzej/ > > > > >> Go one step further: >> >> Clear[f, a, b, k] >> k = $MachinePrecision; >> f = SetAccuracy[333.75*b^6 + >> a^2*(11*a^2*b^2 - b^6 - >> 121*b^4 - 2) + >> 5.5*b^8 + a/(2*b), k]; >> a = 77617.; b = 33096.; >> a = SetAccuracy[a, k]; >> b = SetAccuracy[b, k]; >> f >> Accuracy[f] >> Precision[f] >> >> -5.517164009`0*^19 >> >> Accuracy::mnprec:Value -23 would be inconsistent with $MinPrecision; >> bounding by $MinPrecision instead. >> >> -20 >> >> 0 >> >> Accuracy::mnprec:Value !(-23) would be inconsistent with >> $MinPrecision; >> bounding by $MinPrecision instead. >> >> See that? NO precision. >> >> Bobby >> >> -----Original Message----- >> Sent: Tuesday, October 01, 2002 5:53 PM >> Cc: drbob@bigfoot.com >> >> You are of course right, I forgot that Mathematica does not try to > keep >> >> precision or accuracy of machine arithmetic computations. But I think >> the actual precision you set need not be higher than machine > precision, >> >> it just has to be set explicitely (is that right?). For example: >> >> In[1]:= >> Clear[f,a,b,k] >> >> In[2]:= >> k = $MachinePrecision; >> >> In[3]:= >> f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ >> (2*b),k]; >> >> In[4]:= >> a=77617.;b=33096.; >> >> In[5]:= >> a=SetAccuracy[a,k]; >> >> In[6]:= >> b=SetAccuracy[b,k]; >> >> >> In[7]:= >> f >> >> Out[7]= >> !((-5.51716400890319`-2.8311*^19)) >> >> In[8]:= >> Accuracy[f] >> >> Accuracy::mnprec: Value -23 would be inconsistent with $MinPrecision; > >> bounding by $MinPrecision instead. >> >> Out[8]= >> -20 >> >> >> Andrzej >> >> >> >>> Andrzej >>> >>> Yes, like you I took the original question to be about how to get the >> >>> result >>> of using the naive rational versions in place of the inexact numbers. >>> Bobby raises the question of how we might know accuracy of the > result. >>> >>> You answer this with >>>> a = 77617.; b = 33096.; >>>> >>>> In[2]:= >>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>> 5.5*b^8 + a/(2*b) >>>> >>>> In[3]:= >>>> f >>>> >>>> Out[3]= >>>> -1.1805916207174113*^21 >>>> >>>> In[4]:= >>>> Accuracy[%] >>>> >>>> Out[4]= >>>> -5 >>> >>> However this computation is done in machine arithmetic, which means >>> that >>> Mathematica keeps no check on the accuracy and precision of the >>> computation, >>> and Mathematica gives the default accuracy value without any real >>> signifcance: >>> >>> $MachinePrecision - Log[10,Abs[f]]//Round >>> >>> -5 >>> >>> To get meaningful accuracy and precision values we need to force the >>> computation to be in bignums (bigfloat, arbitrary precision) >>> arithmetic. >>> Mathematica then keeps track of the accuracy and precision that it > can >>> guarantee. >>> >>> Clear[f,a,b,k] >>> k=50; >>> >>> f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ >>> (2*b),k]; >>> a=77617.;b=33096.; >>> a=SetAccuracy[a,k]; >>> b=SetAccuracy[b,k]; >>> f >>> >>> -0.82739605995 >>> >>> Accuracy[f] >>> >>> 11 >>> >>> Precision[f] >>> >>> 11 >>> >>> Which tells us that the error in the computed value of f is not more >> >>> than >>> >>> 10^-11 >>> >>> >>> The above results are rounded. More accurately we get >>> >>> Accuracy[f, Round->False] >>> >>> 11.4956 >>> >>> Precision[f, Round->False] >>> >>> 11.4133 >>> >>> There are several related issues here: >>> - is the precision (accuracy) of rhe input known? >>> - how does Mathematica interpret what one gives it? >>> - how does Mathematica go about the calculation? >>> >>> -- >>> Allan >>> >>> --------------------- >>> Allan Hayes >>> Mathematica Training and Consulting >>> Leicester UK >>> www.haystack.demon.co.uk >>> hay@haystack.demon.co.uk >>> Voice: +44 (0)116 271 4198 >>> >>> >>>> It seems clear to me that what Allan and what you mean by succeeds >>>> here refer to quite different things and your objection is therefore >>>> beside the point. There are obviously two ways in which one can >>>> interpret the original posting. The first interpretation, which > Allan >>>> and myself adopted, was that the question concerned purely the >>>> computational mechanism of Mathematica. Or, to put it in other > words, >>>> it was why are the results of these two computations not the > same?. >>>> In this sense success means no more than making Mathematica return >>>> the same answer using the two different routes the original poster >>>> used. >>>> You on the other hand choose to assume that the posting shows that >> its >>>> author does not understand not just the mechanism of significance >>>> arithmetic used by Mathematica but also computations with inexact >>>> numbers in general. I do not think this is necessarily the correct >>>> assumption. I also don't see that Mathematica is doing anything >> wrong. >>>> After all, one can always check the accuracy of your answer: >>>> >>>> In[1]:= >>>> a = 77617.; b = 33096.; >>>> >>>> In[2]:= >>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>> 5.5*b^8 + a/(2*b) >>>> >>>> In[3]:= >>>> f >>>> >>>> Out[3]= >>>> -1.1805916207174113*^21 >>>> >>>> In[4]:= >>>> Accuracy[%] >>>> >>>> Out[4]= >>>> -5 >>>> >>>> which tells you that it can't be very reliable. What more do you >>>> demand? >>>> >>>> Andrzej >>>> >>>> >>>> Andrzej Kozlowski >>>> Yokohama, Japan >>>> http://www.mimuw.edu.pl/~akoz/ >>>> http://platon.c.u-tokyo.ac.jp/andrzej/ >>>> >>>> >>>>> Actually, we don't know whether SetAccuracy succeeds, because we >>>>> don't >>>>> know how inexact those numbers really are. If they are known to >> more >>>>> digits than shown in the original post, they should be entered with >> >>>>> as >>>>> much precision as they deserve. If not, there's no trick or >>>>> algorithm >>>>> that will give the result precision, because the value of f really >> is >>>>> in the noise. That is, we have no idea what the value of f > should >>>>> be. >>>>> >>>>> Mathematica's failing is in returning a value without pointing out >>>>> that >>>>> it has no precision. >>>>> >>>>> Bobby >>>>> >>>>> -----Original Message----- >>>>> Sent: Monday, September 30, 2002 11:59 AM >>>>> >>>>> Andrzej, Bobby, Peter >>>>> >>>>> It looks as if using SetAccuracy succeeds here because the inexact >>>>> numbers >>>>> that occur have finite binary representations. If we change them >>>>> slightly to >>>>> avoid this then we have to use Rationalize: >>>>> >>>>> 1) Using SetAccuracy >>>>> >>>>> Clear[a,b,f] >>>>> >>>>> >>>>> f=SetAccuracy[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>>>> (2*b), >>>>> Infinity]; >>>>> >>>>> a=77617.1; >>>>> b=33096.1; >>>>> >>>>> a=SetAccuracy[a,Infinity];b=SetAccuracy[b,Infinity]; >>>>> >>>>> f >>>>> >>>>> >>>>> > -1564032114908486835197494923989618867972540153873951942813115514949 >>>>> 3891236234 >>>>> >>>>> >> 525007719168693704591197760187988046304361497869199129319625743010292 >>>>> 36 >>>>> 3 >>>>> 1246 >>>>> 75 >>>>> >>>>> / >>>>> >> 108671061439707605510003578275547938881981431359756495796079898677435 >>>>> 72 >>>>> 8240 >>>>> 16 >>>>> 0653953612982932181371232436367739737604096 >>>>> >>>>> 2) Rewriting as fractions >>>>> >>>>> a=776171/10; >>>>> b=330961/10; >>>>> >>>>> > f=33374/100*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+54/10*b^8+a/(2*b) >>>>> >>>>> -(5954133808997234115690303589909929091649391296257/ >>>>> 41370125000000) >>>>> >>>>> 3) Using Rationalize >>>>> >>>>> Clear[a,b,f] >>>>> >>>>> >>>>> f=Rationalize[333.74*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.4*b^8+a/ >>>>> (2*b), >>>>> 0]; >>>>> >>>>> a=77617.1; >>>>> b=33096.1; >>>>> >>>>> a=Rationalize[a,0];b=Rationalize[b,0]; >>>>> >>>>> f >>>>> >>>>> -(5954133808997234115690303589909929091649391296257/ >>>>> 41370125000000) >>>>> >>>>> >>>>> I use Rationalize[. , 0] besause of results like >>>>> >>>>> Rationalize[3.1415959] >>>>> >>>>> 3.1416 >>>>> >>>>> Rationalize[3.1415959,0] >>>>> >>>>> 31415959/10000000 >>>>> -- >>>>> Allan >>>>> >>>>> --------------------- >>>>> Allan Hayes >>>>> Mathematica Training and Consulting >>>>> Leicester UK >>>>> www.haystack.demon.co.uk >>>>> hay@haystack.demon.co.uk >>>>> Voice: +44 (0)116 271 4198 >>>>> >>>>> >>>>>> Well, first of of all, your using SetAccuracy and SetPrecision > does >>>>>> nothing at all here, since they do not change the value of a or b. >> >>>>>> You >>>>>> should use a = SetAccuracy[a, Infinity] etc. But even then you >> won't >>>>>> get the same answer as when you use exact numbers because of the >>>>>> way >>>>>> you evaluate f. Here is the order of evaluation that will give you >> >>>>>> the >>>>>> same answer, and should explain what is going on: >>>>>> >>>>>> f = SetAccuracy[333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121* >>>>>> b^4 - 2) + 5.5*b^8 + a/(2*b), Infinity]; >>>>>> >>>>>> a = 77617.; >>>>>> >>>>>> >>>>>> b = 33096.; >>>>>> >>>>>> a = SetAccuracy[a, Infinity]; b = SetAccuracy[b, Infinity]; >>>>>> >>>>>> f >>>>>> >>>>>> 54767 >>>>>> -(-----) >>>>>> 66192 >>>>>> >>>>>> Andrzej Kozlowski >>>>>> Toyama International University >>>>>> JAPAN >>>>>> >>>>>> >>>>>> >>>>>> >>>>>>> Could someone explain what is going on here, please? >>>>>>> >>>>>>> In[1]:= >>>>>>> a = 77617.; b = 33096.; >>>>>>> >>>>>>> In[2]:= >>>>>>> SetAccuracy[a, Infinity]; SetAccuracy[b, Infinity]; >>>>>>> SetPrecision[a, Infinity]; SetPrecision[b, Infinity]; >>>>>>> >>>>>>> In[4]:= >>>>>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + 5.5*b^8 > + >>>>>>> a/(2*b) >>>>>>> >>>>>>> In[5]:= >>>>>>> SetAccuracy[f, Infinity]; SetPrecision[f, Infinity]; >>>>>>> >>>>>>> In[6]:= >>>>>>> f >>>>>>> >>>>>>> Out[6]= >>>>>>> -1.1805916207174113*^21 >>>>>>> >>>>>>> In[7]:= >>>>>>> a = 77617; b = 33096; >>>>>>> >>>>>>> In[8]:= >>>>>>> g := (33375/100)*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>>>>> (55/10)*b^8 + a/(2*b) >>>>>>> >>>>>>> In[9]:= >>>>>>> g >>>>>>> >>>>>>> Out[9]= >>>>>>> -(54767/66192) >>>>>>> >>>>>>> In[10]:= >>>>>>> N[%] >>>>>>> >>>>>>> Out[10]= >>>>>>> -0.8273960599468214 >>>>>>> >>>>>>> >>>>>>> PK >>>>>>> >>>>>>> >>>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> >>>> >>>> >>>> >>> >>> >>> >>> >>> >>> >>> >>> >>> >>> >> Andrzej Kozlowski >> Yokohama, Japan >> http://www.mimuw.edu.pl/~akoz/ >> http://platon.c.u-tokyo.ac.jp/andrzej/ >> >> >> >> > > > > > Andrzej Kozlowski Yokohama, Japan http://www.mimuw.edu.pl/~akoz/ http://platon.c.u-tokyo.ac.jp/andrzej/ ==== >>In any case the answer you get is meaningless. Precisely. Bobby -----Original Message----- 121*b^4 - 2) + 5.5*b^8 + a/(2*b), k]; In[4]:= a = 77617.; b = 33096.; In[5]:= a = SetAccuracy[a, k]; In[6]:= b = SetAccuracy[b, k]; In[7]:= f Out[7]= -5.517164009`-2.8311*^19 In[8]:= Accuracy[f] Out[8]= -23 In[9]:= Precision[f] Out[9]= -3 In any case the answer you get is meaningless. > I think I'd prefer that Mathematica gave me a clue -- without being > asked explicitly in JUST the right way that only you know and had > forgotten -- that there's a precision problem. > > Oddly enough, when you DO ask it nicely, you get error messages, but if > you're not aware there's a problem, it lets you go on your merry way, > working with noise. > > Bobby > > -----Original Message----- > Sent: Tuesday, October 01, 2002 8:05 PM > Cc: 'Allan Hayes'; mathgroup@smc.vnet.net > > So? It's just as it should be, isn't it? > > Andrzej > > Andrzej Kozlowski > Toyama International University > JAPAN > http://sigma.tuins.ac.jp/~andrzej/ > > > > >> Go one step further: >> >> Clear[f, a, b, k] >> k = $MachinePrecision; >> f = SetAccuracy[333.75*b^6 + >> a^2*(11*a^2*b^2 - b^6 - >> 121*b^4 - 2) + >> 5.5*b^8 + a/(2*b), k]; >> a = 77617.; b = 33096.; >> a = SetAccuracy[a, k]; >> b = SetAccuracy[b, k]; >> f >> Accuracy[f] >> Precision[f] >> >> -5.517164009`0*^19 >> >> Accuracy::mnprec:Value -23 would be inconsistent with $MinPrecision; >> bounding by $MinPrecision instead. >> >> -20 >> >> 0 >> >> Accuracy::mnprec:Value !(-23) would be inconsistent with >> $MinPrecision; >> bounding by $MinPrecision instead. >> >> See that? NO precision. >> >> Bobby >> >> -----Original Message----- >> Sent: Tuesday, October 01, 2002 5:53 PM >> Cc: drbob@bigfoot.com >> >> You are of course right, I forgot that Mathematica does not try to > keep >> >> precision or accuracy of machine arithmetic computations. But I think >> the actual precision you set need not be higher than machine > precision, >> >> it just has to be set explicitely (is that right?). For example: >> >> In[1]:= >> Clear[f,a,b,k] >> >> In[2]:= >> k = $MachinePrecision; >> >> In[3]:= >> f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ >> (2*b),k]; >> >> In[4]:= >> a=77617.;b=33096.; >> >> In[5]:= >> a=SetAccuracy[a,k]; >> >> In[6]:= >> b=SetAccuracy[b,k]; >> >> >> In[7]:= >> f >> >> Out[7]= >> !((-5.51716400890319`-2.8311*^19)) >> >> In[8]:= >> Accuracy[f] >> >> Accuracy::mnprec: Value -23 would be inconsistent with $MinPrecision; > >> bounding by $MinPrecision instead. >> >> Out[8]= >> -20 >> >> >> Andrzej >> >> >> >>> Andrzej >>> >>> Yes, like you I took the original question to be about how to get the >> >>> result >>> of using the naive rational versions in place of the inexact numbers. >>> Bobby raises the question of how we might know accuracy of the > result. >>> >>> You answer this with >>>> a = 77617.; b = 33096.; >>>> >>>> In[2]:= >>>> f := 333.75*b^6 + a^2*(11*a^2*b^2 - b^6 - 121*b^4 - 2) + >>>> 5.5*b^8 + a/(2*b) >>>> >>>> In[3]:= >>>> f >>>> >>>> Out[3]= >>>> -1.1805916207174113*^21 >>>> >>>> In[4]:= >>>> Accuracy[%] >>>> >>>> Out[4]= >>>> -5 >>> >>> However this computation is done in machine arithmetic, which means >>> that >>> Mathematica keeps no check on the accuracy and precision of the >>> computation, >>> and Mathematica gives the default accuracy value without any real >>> signifcance: >>> >>> $MachinePrecision - Log[10,Abs[f]]//Round >>> >>> -5 >>> >>> To get meaningful accuracy and precision values we need to force the >>> computation to be in bignums (bigfloat, arbitrary precision) >>> arithmetic. >>> Mathematica then keeps track of the accuracy and precision that it > can >>> guarantee. >>> >>> Clear[f,a,b,k] >>> k=50; >>> >>> f=SetAccuracy[333.75*b^6+a^2*(11*a^2*b^2-b^6-121*b^4-2)+5.5*b^8+a/ >>> (2*b),k]; >>> a=77617.;b=33096.; >>> a=SetAccuracy[a,k]; >>> b=SetAccuracy[b,k]; >>> f >>> >>> -0.82739605995 >>> >>> Accuracy[f] >>> >>> 11 >>> >>> Precision[f] >>> >>> 11 >>> >>> Which tells us that the error in the computed value of f is not more >> >>> than >>> >>> 10^-11 >>> >>> >>> The above results are rounded. More accurately we get >>> >>> Accuracy[f, Round->False] >>> >>> 11.4956 >>> >>> Precision[f, Round->False] >>> >>> 11.4133 >>> >>> There are several related issues here: >>> - is the precision (accuracy) of rhe input known? >>> - how does Mathematica interpret what one gives it? >>> - how does Mathematica go about the calculation? >>> >>> -- >>> Allan >>> >>> --------------------- >>> Allan Hayes >>> Mathematica Training and Consulting >>> Leicester UK >>> www.haystack.demon.co.uk >>> hay@haystack.demon.co.uk >>> Voice: +44 (0)116 271 4198 >>> >>> message >>>> It seems clear to me that what Allan and what you mean by succeeds >>>> here refer to quite different things and your objection is therefore >>>> beside the point. There are obviously two ways in which one can >>>> interpret the original posting. The first interpretation, which > Allan >>>> and myself adopted, was that the question concerned purely the >>>> computational mechanism of Mathematica. Or, to put it in other > words, >>>> it was why are the results of these two computations not the > same?. >>>> In this sense success means no more than making Mathematica return >>>> the same answer using the two different routes the original poster >>>> used. >>>> You on the other hand choose to assume that the posting shows that >> its >>>> author does not understand not just the mechanism of significance >>>> arithmetic used by Mathematica but also computations with inexact