A7 == Thanks for your extensive answer but I still have some doubts about convergence of the following integral (m,nintegrers>=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>=0W[m,m+1]=1/2W[m,m+3]=-1/2Numerically we haveNIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, 0, Infinity}]NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy....0.597973NIntegrate[BesselJ[0, x]*BesselJ[1, x], {x, 0, Infinity}, Method ->Oscillatory]NIntegrate::ploss : ....0.5So I define also the corresponding numeric definitionNW[m_, n_] := NIntegrate[BesselJ[m, x]*BesselJ[n, x], {x, 0, Infinity},Method -> Oscillatory]THEORYThe 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 theoryWS[m_,n_,p_]:=Integrate[BesselJ[m, x]*BesselJ[n, x] x^-p, {x, 0, Infinity}]= A/BwhereA=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.ASYMPTOTICSThe Watson theory is in conßict with Mathematica and your notes accordingwhichthe asyntotic trend 1/x of the integrand in W[m,n] is enough forconvergernce. I divide the integral in two partsWasy[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 secondintegralso that I can writeWasy[m_,n_,a_]= int1[m,n,a]+int2[m,n,a]whereint1[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) issingularand 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]) RESULTSm=1;n=0;a=20.;WS[m,n,0]=divergentW[m,n]=1/2NW[m ,n]=0.5Wasy[m,n,a]=.49816m=2;n=0;a=20.;WS[m,n,0]=divergent W[m,n]=0.427599NW[m,n]=-2.43818Wasy[m,n,a]=-1.48052m=3;n =1;a=20.;WS[m,n,0]=divergentW[m,n]=0.639806NW[m,n]=- 2.31957Wasy[m,n,a]=-1.26822m=4;n=0;a=20.;WS[m,n,0]= divergentW[m,n]=-.852012NW[m,n]=1.45786Wasy[m,n,a]= 1.06835The 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.RobertRoberto BrambillaCESIVia Rubattino 5420134 Milanotel +39.02.2125.5875fax +39.02.2125.5492rlbrambilla@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 doneVB>>VB>> and observe that the output you had hadVB>>VB>> W[0,0]=-(2 EulerGamma + Log[4] + 4 PolyGamma[0, 1/2])/(2 Pi)VB>>VB>> = 0.84564VB>>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 severalposting 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 everywhereover the integration region and decays at z -> Infinityas 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])/zwhich 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 thatwe 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 BondarenkoMathematical DirectorSymbolic Testing GroupWeb : 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 aproper color with each value of z.David Parkdjmp@earthlink.nethttp://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 andthen animate them. The problem is that the audience sees the animationtwice. Once when the frames are being generated and then again after I haveclosed the group and double clicked on the top graphic. However, the firstshowing during generation is enough (but uncontrolled). Is it possible to tidy up the generation of the graphic so that it becomesthe animation?There were two main solutions which I now apply to my real problem and notthe simple, generic, problem in the original post. The real problem was how to build up a probably density function (PDF) inreal time. In the examples below I redraw the PDF every time I add 10points.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 thememory. The paradigm here is generate all frames, then animate all frames. Wereally need a loop that does: generate next frame, delete last frame, show next frameIs it possible to do this?Bobby Treat also offered a solution involving SelectionMove.The second solution was from Jerry Blimbaum and uses JAVAInstallJava[];UseFrontEndForRendering = False;createWindow[] := Module[{frame}, frame = JavaNew[com.wolfram.jlink.MathFrame, Probability DensityFunction]; 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 havenot managed to control the speed of this animation. The JAVA code sleep doesnot work nor does the use of a Mathematica Pause which always rounds up toan integer (is this a bug?).Hugh Goyder ==== >>I need a loop that goes generate next frame, delete old frame, shownew 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, otherthan the trick of taking advantage of the half-period.I'll be very interested in a solution myself, as I frequently run out ofmemory in animations of fairly modest size -- despite having 1024MB ofRAM.Bobby-----Original Message-----need a loop that goes generate next frame, delete old frame, show new frameso 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 GoyderThis 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 ConsultingThe final answer to your Mathematica needsSpend 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'tknow how inexact those numbers really are. If they are known to moredigits than shown in the original post, they should be entered with asmuch precision as they deserve. If not, there's no trick or algorithmthat will give the result precision, because the value of f really isin 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 thatit 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 - 15640321149084868351974949239896188679725401538739519428131155 149493891236234 52500771916869370459119776018798804630436149786919912931962574 3010292363124675/ 10867106143970760551000357827554793888198143135975649579607989 867743572824016 06539536129829321813712324363677397376040962) 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 HayesMathematica Training and ConsultingLeicester UKwww.haystack.demon.co.ukhay@haystack.demon.co.ukVoice : +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]:=fOut[3]=-1.1805916207174113*^21In[4]:= Accuracy[%]Out[4]=-5which tells you that it can't be very reliable. What more do you demand?AndrzejAndrzej KozlowskiYokohama, Japanhttp://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> - 15640321149084868351974949239896188679725401538739519428131155 14949> 3891236234 52500771916869370459119776018798804630436149786919912931962574 301029236 > 3> 1246> 75 / > 10867106143970760551000357827554793888198143135975649579607989 867743572> 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 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, yourWhen you do calculations with arbitrary-precision numbers, asdiscussed in the previous section, Mathematica always keeps track ofthe precision of your results, and gives only those digits which areknown to be correct, given the precision of your input. When you docalculations with machine-precision numbers, however, Mathematicaalways gives you a machine-[CapitalEth]precision result, whether or not all thedigits in the result can, in fact, be determined to be correct on thebasis of your input. In practice, to rely on a numerical result, you ALWAYS have to checkits 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]=1180591620717411303424In[4]:=Accuracy[f] Out[4]=InfinityWe got completely wrong result with Infinite accuracy. In conclusion,one can not rely on Accuracy. Checking numerical results inMathematica 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/> [...]> ==== AndrzejYes, like you I took the original question to be about how to get the resultof 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]=> -5However this computation is done in machine arithmetic, which means thatMathematica keeps no check on the accuracy and precision of the computation,and Mathematica gives the default accuracy value without any realsignifcance: $MachinePrecision - Log[10,Abs[f]]//Round -5To get meaningful accuracy and precision values we need to force thecomputation to be in bignums (bigßoat, arbitrary precision) arithmetic.Mathematica then keeps track of the accuracy and precision that it canguarantee. 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] 11Which tells us that the error in the computed value of f is not more than 10^-11The above results are rounded. More accurately we get Accuracy[f, Round->False] 11.4956 Precision[f, Round->False] 11.4133There 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 HayesMathematica Training and ConsultingLeicester UKwww.haystack.demon.co.ukhay@haystack.demon.co.ukVoice : +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> - 15640321149084868351974949239896188679725401538739519428131155 14949> 3891236234> 52500771916869370459119776018798804630436149786919912931962574 301029236> 3> 1246> 75>/> 10867106143970760551000357827554793888198143135975649579607989 867743572> 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 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 randomsubstitution failures in 4.2. Question: how can thesubstitution 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]] and1/Sqrt[a] is interal Power[a,Rational[-1,2]] and so a rule Sqrt[a]->q will not match with1/Sqrt[a]. You needf = 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) Jensought 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 preparationfor conversion of those expressions to C. The letters will be namesof 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 <- operatorfor 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}] // NResult: 2.74716Now I wonder if this sum will converge or keep on growing, albeit veryslowly.Matthias BodeSal. Oppenheim jr. & Cie. KGaAKoenigsberger Strasse 29D-60487 Frankfurt am MainGERMANYMobile: +49(0)172 6 74 95 77Internet: 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, bythe 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^xOut[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 KozlowskiToyama International UniversityJAPANhttp://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 theminimizing value (arg min) of another function.For a simple example consider the followingfunction 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 onlya 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 LudsteckEconomics DepartmentUniversity of RegensburgUniversitaetsstrasse 3193053 RegensburgReply-To: kuska@informatik.uni-leipzig.de ==== you guess right andif you hinder Mathematica toevaluate opt[] for symbolicarguments, withopt[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 KozlowskiYokohama, Japanhttp://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 useMenuCellDefault Inline Format TypeStandard Form. Also, when I design myown style sheets I often define a bolder font for Inline cells so they standout better.David Parkdjmp@earthlink.nethttp://home.earthlink.net/~djmp/ soon as I press the Control-^ key combination, an Inline cell is createdbeginning with the x, and then when I type the exponent 2 everything inthat Inline cell is now in Times, and the x is Italic. To change bothcharacters to Courier is not so easy: it seems to require separatelythe entire Inline cell and selecting Courier does not change the exponent!)So to avoid this annoyance I normally must first type the desiredexpression in a separate Input cell, then copy the contents of that cellto 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.eduMathematics & Statistics Dept.Lederle Graduate Research Tower phone 413 549-1020 (H)University of Massachusetts 413 545-2859 (W)710 North Pleasant StreetAmherst, MA 01375Reply-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.eduMathematics & Statistics Dept.Lederle Graduate Research Tower phone 413 549-1020 (H)University of Massachusetts 413 545-2859 (W)710 North Pleasant StreetAmherst, MA 01375 ==== What is the TMJ style sheet? In any case, when people design new stylesheets 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 keysin the order that the cell definitions appear in the style sheet. That meansthat new cell styles should be moved lower in the style sheet, even if thatwould not be their natural order. It would be nice if WRI allowed us toexplicitly assign the keys for the styles.David Parkdjmp@earthlink.nethttp://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.netApproved: 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 codeproduces rotated text with each of the charaters rotated as welldegstr[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 butthe characters in the text are not. That is the position of the chactersrelative to each other changes but their orientation remains constant.I would like to have the character orientation change as it does when thiscode is run under Windows NT. Does anyone know if there is a settingsomewhere that controls character orientation when text is rotated? ==== 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)fAccuracy[f]Precision[f]Accuracy[10.^21] Precision[10.^21]1.1805916207174113*^21-516-516The two numbers supposedly have the same accuracy and precision, yet thefirst is in doubt by about 22 ORDERS OF MAGNITUDE -- never mind thedigits! Mathematica computed this beast without giving any indicationit was just noise -- without REALIZING it was noise.>>What more do you demand?I'm not demanding, objecting, or criticizing. I'm pointing out theproblem. On the one hand, the poster is computing something that'ssimply not well-behaved. Unless he knows the coefficients with VERYhigh precision, he can't know even the magnitude of the result -- andthat's not Mathematica's fault at all. On the other hand, Mathematicadoesn't notice that precision is lost in the computation, and perhaps itshould. 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]:=fOut[3]=-1.1805916207174113*^21In[4]:= Accuracy[%]Out[4]=-5which tells you that it can't be very reliable. What more do you demand?AndrzejAndrzej KozlowskiYokohama, Japanhttp://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 outthat> 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> - 15640321149084868351974949239896188679725401538739519428131155 14949> 3891236234 52500771916869370459119776018798804630436149786919912931962574 301029236 > 3> 1246> 75 / > 10867106143970760551000357827554793888198143135975649579607989 867743572> 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 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 youthe> 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 inexactnumbers: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*dyf is sensitive to inaccuracy in the various numbers to widely varyingdegrees. Since the correct answer is small, we need a LOT ofprecision in the inputs to get there. If any of the inputs are merelymachine precision numbers, we have NO precision in the result.The second and third terms are nearly the same magnitude with differentsigns. Even worse, the first term almost perfectly fills the gap:a = 77617; b = 33096;(33375/100)*b^6438605750846393161930703831040a^2*(11*a ^2*b^2 - b^6 - 121*b^4 - 2)-7917111779274712207494296632228773890(55/10)*b^ 87917111340668961361101134701524942848% + %% + %%%-2Bobby Treat-----Original Message-----b = SetAccuracy[33096., Infinity]; In[4]:=fOut[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]=1180591620717411303424Similarily: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]:=fOut[4]=-0.8273960599468212641107299556`11.4133In [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 asimplified 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 plusx[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 and9/16 to real machine numbers. I consider this to be a bug.Now, for an example more representative of the situation that I've beencoming 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 precisionnumbers. The simplified example above demonstrates that the problem existswhen using machine numbers. Now, we'll see what happens when we usearbitrary precision numbers. First, let's simplify the expression with 12terms.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]) / 4096As you can see, a sum with 12 terms upon simplification has coefficientswhich are still integers as they should be. However, increasing the numberof terms to 13 yieldsSimplify[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 thebug from the simplified example. Finally, let's see what happens when wehave 55 terms. Rather than putting the resulting expression here, I willjust leave it at the end of the post. The result though is somewhatsurprising. Each of the coefficients of the x[i] are again real numbers, butnow their precision is only 0! The proper result of course is the sum ofsome real number (with a precision close to 0 due to numerical cancellation)and an expression consisting of rational numbers multiplied by x[i]. Theloss of precision of the coefficients of the x[i] is precisely what occurredto me. Of course, in this simple example, simply using Expand instead ofSimplify produces the expected result, and is my workaround. I hope youagree with me that this is a bug, and one that Wolfram needs to correct.Carl WollPhysics DeptU of WashingtonSimplify[expr55] -16 -1 -1 -10. 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. 10x[15] + -5 -5 -5 -6 -6 0. 10 x[16] + 0. 10 x[17] + 0. 10 x[18] + 0. 10 x[19] + 0. 10x[20] + -6 -6 -7 -7 -7 0. 10 x[21] + 0. 10 x[22] + 0. 10 x[23] + 0. 10 x[24] + 0. 10x[25] + -8 -8 -8 -9 -9 0. 10 x[26] + 0. 10 x[27] + 0. 10 x[28] + 0. 10 x[29] + 0. 10x[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. 10x[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 toreplace 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 everycase. 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 willprobably 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 other14).> These equation are not linearly related.. the highest degree in anyone eqn> is degree 4 i believe.. and there are some cross terms in theequations but> not every equation depends on every variable.. (some are actuallyrather> 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[Terminate];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]:=fOut[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]=-20Andrzej> 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 (bigßoat, 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 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> - 15640321149084868351974949239896188679725401538739519428131155 14949> 3891236234 52500771916869370459119776018798804630436149786919912931962574 3010292 > 36> 3> 1246> 75 /> 10867106143970760551000357827554793888198143135975649579607989 8677435 > 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 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> PKAndrzej KozlowskiYokohama, Japanhttp://www.mimuw.edu.pl/~akoz/http:// platon.c.u-tokyo.ac.jp/andrzej/ ==== So? It's just as it should be, isn't it?AndrzejAndrzej KozlowskiToyama International UniversityJAPANhttp://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 (bigßoat, 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 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> - 15640321149084868351974949239896188679725401538739519428131155 14949> 3891236234> 52500771916869370459119776018798804630436149786919912931962574 3010292> 36> 3> 1246> 75 / 10867106143970760551000357827554793888198143135975649579607989 8677435> 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 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/ ==== By taking sums and differences of the counterweights 1, 3, ... 3^(k-1)you can precisely express all integers up to the integer whose base 3notation is composed of k ones.By doubling those counterweights, you can precisely express all the EVENintegers up to TWICE that limit, or -1 + 3^k. If you know that unknownsare limited to the integers from 1 to 3^k, this allows you to preciselyweight the even numbers and bracket the odd numbers. Hencecounterweights 2, 6, ... 2*3^(k-1) are sufficient for unknowns up to3^k.Notice that this is an efficient coding scheme: Each counterweight ismultiplied by -1, 0, or 1, so there can be at most 3^k distinctcounterweight sums. More than half are negative or zero, however; thepositive sums and differences number at most (3^k - 1)/2. Since we needto enumerate (and CAN enumerate with even weights) only the even numbersand because we get the maximum unknown 3^k by elimination, we cover upto twice as many as we can enumerate, plus one. Twice (3^k - 1)/2 plusone is 3^k, so we're getting the maximum coverage possible for a givennumber of counterweights.Bobby Treat-----Original Message-----For these 40 positive integers to match 1 to 40 correspondingly, thebiggestshould be 40 (also, the smallest should be 1). So we have:a+b+c+d=40 or 1+b+c+d=40 (we state that a I have a very interesting math problem:If I have a scales,and I> have 40 things that their mass range from 1~40 which each is anature> number,and now I can only make 4 counterweights to measure out each> mass of those things.Question:What mass should the counterweights> be???> The answer is that 1,3,9,27 and I wnat to use mathematica tosolve> this problem.> In fact,I think that this physical problem has various> answer,ex.2,4,10,28> this way also work,because if I have a thing which weight 3 , and I> can measure out by comparing 2<3<4 . But,If I want to solve thismath> problem:> {x|x=k1*a+k2*b+k3*c+k4*d}={1,2,3,4,,,,,,40} where a,b,c,d is nature> numbers.> and {k1,k2,k3,k4}={1,0,-1}> How to solve it ??> mathematica solving method. appreciate any idea sharing> sincerely> bryan Just use brute force. Needs[DiscreteMath`Combinatorica`]; var = {a, b, c, d}; n = Length[var]; s = Outer[Times, var, {-1, 0, 1} ]; f = Flatten[Outer[Plus, Sequence@@s]]; Since the length of f is just 3^n then the range of numbers> to be covered should be {-(3^n-1)/2, (3^n-1)/2}.> Consequently, the largest of the weights can not exceed> (3^n-1)/2 - (1+2+...+(n-1)) or ((3^n-1) - n(n-1))/2 34 Thread[var->#]& /@ (First /@ Select[{var,f} /. Thread[var->#]& /@ KSubsets[Range[((3^n-1) - n(n-1))/2], n], Sort[#[[2]]] == Range[-(3^n-1)/2,(3^n-1)/2]&]) {{a -> 1, b -> 3, c -> 9, d -> 27}}> Bob Hanlon> Chantilly, VA USA> ==== Plot[Evaluate[Sin[x]/Sin[x]], {x, -2.5, -1.5}]Gives the expected graph.Bobby Treat-----Original Message-----substituting a ßoating-point value for x results in composite divisionbeing performed numerically, with consequent inaccuracy. This can beverified by the following:In[1]:= FullForm[HoldForm[Sin[x]/Sin[x]]Out[1]//FullForm= HoldForm[Times[Sin[x], Power[Sin[x], -1]]]By default, Mathematica scales the y axis such that the smalldiscrepanciesare visible. The PlotRange option may be given to explicitly specify therange for the y axis (for example PlotRange->{0, 2}).Ian McInnes.> I would expect that Plot[Sin[x]/Sin[x], {x, -2.5, -1.5}] would produce a horizontal line at y=1. However, on my Windows XP computer it produces a graph where the yvalueis> less than 1 at several points. Most notibly between x=-1.6 and x=-1.8 Is this just an isolated case? Or does it happen to others? If so -why? Ken.> ==== >>But my actual output reproduces the input form, i.e., it is simplyf[...listA,newEntry,...].That means your arguments don't fit a pattern for which the function isdefined. That is, it has nothing to do with what's on the right side of:=. It's about making sure there are the right number and type ofarguments passed to the function. You haven't told us what the patternis or what the arguments are, so we can't help you.Bobby-----Original Message----- I define a function of the form: f[...,listA_, newEntry_,...]:= Which[...]Here Which has the form of a set of logically exclusive and exhausivetests, each test with its ownaction that modifies the concrete list subsituted for listA_. For instance, if there were only two tests the rhs above wouldbeone that puts the newEntryat the start of listA and the other that puts it at the end: Which[ test1, listA = Insert[listA, newEntry, 1], test2, listA = Insert[listA, newEntry, -1] ]Then I try to use the function, substituting values for the arguments:f[..., listA, newEntry, ...]. Theoutput should be the appropriately modified list. But my actual outputreproduces the input form,i.e., it is simply f[...listA,newEntry,...]. Yet, the program does do something, as the use of evaluate inplaceshows: Suppose that test1 is satisfied by the particular values of thearguments. Then whenI apply evaluate in place to the lhs of the action that is supposed tooccur when test1 is true,in fact the value of Insert[listA,newEntry,1] is correct, but the rhsgivesonly the initial value oflistA and not the modified value that is on the rhs. Concretely, suppose that in the function argument listA is {121}andnewEntry is 200.Then after evaluation of the function (input to the kernel), that linereads, according toevaluate in place: true, {121} = {200, 121}So it seems that only part of the correct action was taken: 200 was putatthe start of listA. Butthe assignment of this value -- the extended list -- to be the new valueoflistA did not take place. I then tried using a new name for the modified list, substituting thisprogram line for test1: test1, newListA = Insert[listA,newEntry,1]What I get here, applying evaluate in place after evaluating thefunctionwith concrete values as aboveis: true, newListA = {200, 121} Trace doesn't help because all I get back is the function namewithits concrete arguments. Any thoughts? Tom ==== The problem is that I do not think there is any bounded 3d body described by your conditions. Anyway, this is how you can use Mathematica (in general) to solve this sort of problem.You need two packages:In[1]:=< 3x && y < 4 - x^2 && z < x^2 + 4, {x}, {y}, {z}]However, you will just get some error messages and a picture that looks two dimensional. If you wanted the volume, evaluate:In[3]:=Integrate[Boole[y>3x&&y<4-x^2&&z3x&&y<4-x^2&&z3x&&y<4-x^2&&z equation :> z=x^2 +4 (as bottom surface)> and on the xy plane which bounded by a parabola y=4-x^2 and y=3x line. How would I use Mathematica to plot out this 3D object or find out its > volumn with only the equation given? Shz Shz> Andrzej Kozlowski 3/September/2002 04:33pm Mathematica could do this sort of thing if there were a three> dimensional object described by your equations (as boundaries) but> there isn't one. More precisely, the pair of equations {z=x^2 +4,> y=4-x^2} describes a parabola in three space which you can plot with: g1=ParametricPlot3D[{x, 4 - x^2, x^2 + 4}, {x, -5, 5}] The equation y=3x describes the plane: g2 = ParametricPlot3D[{x, 3x, z}, {x, -5, 5}, {z, -25, 25}]> You can see the two together in> < Show[{g1,g2}] There are clearly two points of intersection. They can be found with:> Solve[{z == x^2 + 4, y == 4 - x^2, y == 3*x}, {x, y, z}]> {{z -> 5, y -> 3, x -> 1}, {z -> 20, y -> -12, x -> -4}} So where is the 3D object whose volume you want to find? Andrzej Kozlowski> Toyama International University> JAPAN Can I use Mathematica to find out the volumn of this 3 dimensional> object from> the equations : z=x^2 +4, y=4-x^2, y=3x> Shz Shz ==== How can I run mathematica program in background? Is there any waysimilar to fortran or c where one can runs his program in thebackground? Please suggest. RajReply-To: kuska@informatik.uni-leipzig.de ==== math < someresults.m &???BTW this is a problem of your operating systemand not of Mathematica. Jens > How can I run mathematica program in background? Is there any way> similar to fortran or c where one can runs his program in the> background? Please suggest. Raj ==== All works fine but now I wanted to generate an animated gif with help of Now the same file which worked earlier produces errors:In[4]:=Export[schwing.gif, schwing, GIF]LinkConnect::linkc: !([LeftSkeleton] 1 [RightSkeleton]) is dead; attempt to connect failed.Out[4]=schwing.gifAlso a dialogue box appears withsuch file or directorysuch file or directoryThe files are at the place and are executable etc. But if one tries to execute them always this No such file or directory error message appears.Does anyone know what to do?-- Hendrik van Hees Fakult.8at f.9fr Physik http://theory.gsi.de/~vanhees/ D-33615 Bielefeld ==== First define the functionsz[x_] := x^2 + 4;y1[x_] := 4 - x^2;y2[x_] := 3x;Figure out where the y limits coincide:Solve[y1[x] == y2[x], x]{{x -> -4}, {x -> 1}}Which is bigger in between?y1[0] - y2[0]4y1 is. Confirm that with a plot. (y2 is linear):Plot[{y1[x], y2[x]}, {x, -5, 5}]Next figure out where z is zero:Solve[z[x] == 0, x]{{x -> -2*I}, {x -> 2*I}}It's never zero for real x, but is it positive, or negative?z[0]4Look at the plot, just for fun:Plot[{z[x]}, {x, -4, 1}]z is positive for all x (but only the interval [-4,1] MATTERS). Thevolume you want, therefore, isIntegrate[z[x]*(y1[x] - y2[x]), {x, -4, 1}]625/4or:g[x_] = Integrate[z[x]*(y1[x] - y2[x]), x]g[1] - g[-4]16*x - 6*x^2 - (3*x^4)/4 - x^5/5625/4Bobby Treat-----Original Message-----Sender: steve@smc.vnet.netApproved: Steven M. Christensen , Moderator ==== >-----Original Message----->Sent: Wednesday, September 04, 2002 8:56 AM>previous output?>>What is a efficient way to use output expression to define new >function. PeterPeter,I'd say none,... Ôcause look at In[1]:= y = x^2Out[1]= x^2In[2]:= f1[x_] = Out[1]Out[2]= x^2In[3]:= f2[x_] := Out[1]In[4]:= x = 5Out[4]= 5In[5]:= f3[x_] = Out[1]Out[5]= 25In[6]:= Block[{x}, f4[x_] = Out[1]]Out[6]= 25In[7]:= Block[{x}, f5[x_] := Out[1]]In[9]:= f6[x_] = yOut[9]= 25In[10]:= Block[{x}, f7[x_] = y]Out[10]= 25In[11]:= f8[x_] := yIn[12]:= Through[{f1, f2, f3, f4, f5, f6, f7, f8}[#]] & /@ {0, 1}Out[12]= {{0, 25, 25, 0, 25, 25, 0, 25}, {1, 25, 25, 1, 25, 25, 1, 25}}In[13]:= Block[{x}, f9[x_] = Out[1]]Out[13]= %1In[14]:= Through[{f1, f2, f3, f4, f5, f6, f7, f8, f9}[#]] & /@ {0, 1}Out[14]= {{0, %1, 25, 0, %1, 25, 0, 25, %1}, {1, %1, 25, 1, %1, 25, 1, 25, %1}}f1, f4 and f7 come out right, but f1 is based on the condition that thepattern variables of the expression have not been redefined, and f4 is onlycorrect if the history is remembered well. So what rests is f7, but thatwasn't you question.--Hartmut Wolf ==== If we are in the business of giving good links to Doppler demonstrations,here is another one (freeware):http://www.phy.ntnu.edu.tw/java/Doppler/ Doppler.htmlFor anyone interested in good Java implementations of physics phenomena,NTNU Virtual Physics Laboratory athttp://www.phy.ntnu.edu.tw/java/index.html is a treasure chest.This is not Mathematica stuff (as far as I know), but might be of interestto several Mathgroup readers.Ingolf DahlChalmers UniversitySwedenf9aid@fy.chalmers.se-----Original Message-----Sender: steve@smc.vnet.netApproved: Steven M. Christensen , Moderator ==== My main application for this is in drawing a piece of a plane at a pointnormal to a vector. For me, efficiency is not an important criterion. But Iwant orthogonal vectors so the plane will be composed of squares and notdiamonds. I want unit vectors so that I can easily specify the size of thepiece I am drawing.I posted a summary last evening and it doesn't seem to have appeared today.In any case, I have realized that none of the routines really require theLinearAlgebra`Orthogonalization package, and run faster without it. I havemodified the original routines to use a pure function instead of Normalizeor GramSchmidt.Here is a summary of the responses so far. The timings were done with thefollowing exact vector. They are faster with a machine precision vector.vtest = {1, 2, 1};__________________________________________________________ Hugh Goyder...The original routine wasOrthogonalUnitVectors[v : {_, _, _}] := GramSchmidt[NullSpace[{v}]]but I changed it toOrthogonalUnitVectors[v : {_, _, _}] := #/Sqrt[#.#] & /@ NullSpace[{v}]OrthogonalUnitVectors[vtest]{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {-(1/Sqrt[3]), 1/Sqrt[3], -(1/Sqrt[3])}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{2.53 Second, Null}This appears to be the shortest and fastest routine. I would give it theprize for elegance. It is certainly easy to remember.____________________________________________________ ____Ingolf Dahl...OrthogonalUnitVectors[v:{_, _, _}] := (#1/Sqrt[#1 . #1] & ) /@ NullSpace[{v, {0, 0, 0}, {0, 0, 0}}]OrthogonalUnitVectors[vtest]{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {-(2/Sqrt[5]), 1/Sqrt[5], 0}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{2.59 Second, Null}This is really the same as Hugh Goyder's routine. It isn't necessary to addthe zero rows to the matrix. NullSpace automatically does it.________________________________________________________ Selwyn Hollis...OrthogonalUnitVectors[v : {_, _, _}] := Module[{u, w}, u = Which[(w = {0, v[[3]], -v[[2]]}).w != 0, w, (w = {v[[3]], 0, -v[[1]]}).w != 0, w, (w = {v[[2]], -v[[1]], 0}).w != 0, w, True, Return[$Failed]]; #/Sqrt[#.#] & /@ {u, Cross[u, v]} ]OrthogonalUnitVectors[vtest]{{0, 1/Sqrt[5], -(2/Sqrt[5])}, {Sqrt[5/6], -Sqrt[2/15], -(1/Sqrt[30])}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{11.92 Second, Null}I wish I understood the reason for using rows of the v Cross matrix as testvectors. It does have the advantage that a dot product is used to give ascalar test for zero. Why not...OrthogonalUnitVectors[v : {_, _, _}] := Module[{u, w}, u = Which[(w = {1, 0, 0})[Cross]v != {0, 0, 0}, w, (w = {0, 1, 0})[Cross]v != {0, 0, 0}, w, True, Return[$Failed]]; #/Sqrt[#.#] & /@ {u, Cross[u, v]} ]OrthogonalUnitVectors[vtest]{{1, 0, 0}, {0, -(1/Sqrt[5]), 2/Sqrt[5]}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{10.33 Second, Null}________________________________________________________ ____John Browne...OrthogonalUnitVectors[v : {_, _, _}] := Module[{r, v1, v2}, r = {Random[], Random[], Random[]}; v1 = Cross[v, r]; v2 = Cross[v1, v]; {v1/Sqrt[Dot[v1, v1]], v2/Sqrt[Dot[v2, v2]]}]OrthogonalUnitVectors[vtest]{{0.26601, -0.534209, 0.802408}, {-0.873254, 0.218984, 0.435286}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{9.17 Second, Null}This has the problem that the r vector might be parallel to v. But in 100000test cases it never happened.____________________________________________________ ________Daniel Lichblau...(Daniel left out the normalization, which I added.)OrthogonalUnitVectors[v:{_, _, _}] := With[{vecs = NullSpace[{v}]}, (#1/Sqrt[#1 . #1] & ) /@ {vecs[[1]], vecs[[2]] - vecs[[2]] . vecs[[1]]*vecs[[1]]}]OrthogonalUnitVectors[vtest]{{-(1/Sqrt[ 2]), 0, 1/Sqrt[2]}, {0, 1/Sqrt[5], -(2/Sqrt[5])}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{3.35 Second, Null}________________________________________________________ ______Garry HelzerOrthogonalVectors[v : {a1_, a2_, a3_}] := With[{w = First[Sort[{{a2, -a1, 0}, {a3, 0, -a1}, {0, a3, -a2}}, OrderedQ[{Plus @@ Abs[#2], Plus @@ Abs[#1]}] &]]}, {w, Cross[v, w]} ]OrthogonalUnitVectors[vtest]{{-(1/Sqrt[2]), 0, 1/Sqrt[2]}, {-(2/Sqrt[5]), 1/Sqrt[5], 0}}Do[OrthogonalUnitVectors[vtest], {10000}] // Timing{2.64 Second, Null}________________________________________________________ __________Tom BurtonA nice feature of the following simple scheme is that you know that the twopoints of discontinuity are where r and -r cross the unit sphere.triad[u_, r_:{1, 0, 0}] := Module[{v = Cross[u, r]}, #/Sqrt[# . #] & /@{u, v, Cross[u, v]}]triad[vtest]{{1/Sqrt[6], Sqrt[2/3], 1/Sqrt[6]}, {0, 1/Sqrt[5], -(2/Sqrt[5])}, {-Sqrt[5/6], Sqrt[2/15], 1/Sqrt[30]}}Do[triad[vtest], {10000}] // Timing{17.02 Second, Null}David Parkdjmp@earthlink.nethttp://home.earthlink.net/~djmp/ djmp@earthlink.nethttp://home.earthlink.net/~djmp/ ==== CeZaR,Sort of confused about what you really wanted. If you have forgotten toassign x = ? then you can do the following,Clear[x, eq, m, l]eq = m x^2 + 2(m + 1)x + m + 2l = {};x = 1;For[i = -30, i < 30, i++, AppendTo[l, eq /. m -> i]]ListPlot[l]But if you really intend to make a list of quadratics and then plot theresults for each of the quadratics then you can try the followingClear[x, eq, m, l]eq = m x^2 + 2(m + 1)x + m + 2;l = {};For[i = -30, i < 30, i++, AppendTo[l, eq /. m -> i]];test = Table[l[[i]], {i, 1, Length[l], 1}, {x, -50, 50}];MultipleListPlot[test[[5]], test[[10]], test[[15]], test[[20]],test[[25]], test[[30]], test[[35]], test[[40]], test[[45]], test[[50]], test[[55]], test[[60]], PlotJoined -> True, SymbolShape -> None, ImageSize -> 500 ]Essentially your For loop generates a list of quadratics if x is left out.Yas I want to do these operations: eq = m x^2 + 2(m + 1)x + m + 2> l = {}> For[i = -30, i < 30, i++, AppendTo[l, eq /. m -> i]]> Plot[l, {x, -50, 50}] but i get this error:> l is not a machine-size real number at x = -49.99999583333334 Can someone exaplain to me what i am doing wrong here??> CeZaR ==== Plot[l // Evaluate, {x, -30, 30}, Frame -> True];CeZaR ==== I want to get some F and R such that: F[n,p] + R[n,p] = Sum[Binomial[n,k] p^(n-k) (1-p)^k, {k, 0, Floor[n/2] - 1}],> when F[n,p] is an approximation to the sum and the R is the remaining error. Constantine. >>I'm looking for a way of finding the approximation for partitial binomial>>sum.>>I'll be pleasant for any hint..> [...]> Office: Taub 411You can get a closed form in terms of special functions if you splitinto two cases depending on whether n is even or odd.In[39]:= n = 2*m;In[40]:= InputForm[Sum[Binomial[n,k]*p^(n-k)*(1-p)^k, {k,0,m-1}]]Out[40]//InputForm= p^(2*m)*((p^(-1))^(2*m) - ((-4 + 4/p)^m*Gamma[1/2 + m]* Hypergeometric2F1[1, -m, 1 + m, (-1 + p)/p])/(Sqrt[Pi]*Gamma[1 +m]))In[41]:= n = 2*m+1;In[42]:= InputForm[Sum[Binomial[n,k]*p^(n-k)*(1-p)^k, {k,0,m-1}]]Out[42]//InputForm= p^(1 + 2*m)*((p^(-1))^(1 + 2*m) - (2^(1 + 2*m)*(-1 + p^(-1))^m*Gamma[3/2+ m]* Hypergeometric2F1[1, -1 - m, 1 + m, (-1 + p)/p])/(Sqrt[Pi]*Gamma[2 +m]))Daniel LichtblauWolfram Research ==== I want to get some F and R such that:F[n,p] + R[n,p] = Sum[Binomial[n,k] p^(n-k) (1-p)^k, {k, 0, Floor[n/2] - 1}],when F[n,p] is an approximation to the sum and the R is the remaining error.Constantine.>>I'm looking for a way of finding the approximation for partitial binomial >>sum.>>I'll be pleasant for any hint..>Use the standard add-on package Statistics`NonlinearFit` to do a >NonlinearFit to whatever model you want to use for the approximation.>Bob Hanlon>Chantilly, VA USAConstantine ElsterComputer Science Dept.Technion I.I.T.Office: Taub 411 ==== Constantine,If you break your problem up into two cases, even n and odd n, thenMathematica can sum up your problem and get results, albeit withhypergeometric functions. Consider the following (make sure you look at thiswith a fixed font):In[21]:=evenans = Sum[Binomial[2*n, k]*p^(2*n - k)*(1 - p)^k, {k, 0, n - 1}];In[22]:=PowerExpand[FunctionExpand[FullSimplify[evenans , n [Element] Integers]]]Out[22]= 2 n n n 1 p - 1 2 (1 - p) p Gamma[n + -] Hypergeometric2F1[1, -n, n + 1, -----] 2 p1 - ------------------------------------------------------------- ------- Sqrt[Pi] Gamma[n + 1]In[23]:=oddans = Sum[Binomial[2*n + 1, k]*p^(2*n + 1 - k)*(1 - p)^k, {k, 0, n - 1}];In[24]:=PowerExpand[FunctionExpand[FullSimplify[oddans , n [Element] Integers]]]Out[24]= 2 n + 1 n n + 1 3p - 1 2 (1 - p) p Gamma[n + -] Hypergeometric2F1[1, -n - 1, n +1, -----] 2p1 - ------------------------------------------------------------ -------------------- Sqrt[Pi] Gamma[n + 2]Is this what you were looking for?Carl WollPhysics DeptU of Washington> I want to get some F and R such that: F[n,p] + R[n,p] = Sum[Binomial[n,k] p^(n-k) (1-p)^k, {k, 0, Floor[n/2] -1}],> when F[n,p] is an approximation to the sum and the R is the remainingerror. Constantine.>>I'm looking for a way of finding the approximation for partitialbinomial>>sum.>>I'll be pleasant for any hint..>Use the standard add-on package Statistics`NonlinearFit` to do a>NonlinearFit to whatever model you want to use for the approximation.>>Bob Hanlon> >Chantilly, VA USA Constantine Elster> Computer Science Dept.> Technion I.I.T.> Office: Taub 411 ==== >I want to do these operations:>>eq = m x^2 + 2(m + 1)x + m + 2>l = {}>For[i = -30, i < 30, i++, AppendTo[l, eq /. m -> i]]>Plot[l, {x, -50, 50}]>>but i get this error:>l is not a machine-size real number at x = -49.99999583333334>>Can someone exaplain to me what i am doing wrong here??>eq = m x^2 + 2(m + 1)x + m + 2;l = {};For[i = -30, i < 30, i++, AppendTo[l, eq /. m -> i]];It would be more straightforward and efficient to define the list using Tablel == Table[eq, {m, -30, 29}]TrueTo Plot, use EvaluatePlot[Evaluate[l], {x, -50, 50}];Bob Hanlon ==== >What is a efficient way to use output expression to define new function.Use the menu command Input/Copy Output from Aboveor use % (Out)Bob Hanlon ==== >-----Original Message----->Sent: Wednesday, September 04, 2002 8:57 AM>The problem is that I do not think there is any bounded 3d body >described by your conditions. Anyway, this is how you can use >Mathematica (in general) to solve this sort of problem.>You need two packages:>>In[1]:=><>In[2]:=><< Calculus`Integration`>To see your object use:>InequalityPlot3D[y > 3x && y < 4 - x^2 && z < x^2 + 4, {x}, {y}, {z}]>>However, you will just get some error messages and a picture >that looks >two dimensional. If you wanted the volume, evaluate:>>In[3]:=>Integrate[Boole[y>3x&&y<4-x^2&&z>Out[3]=>-Infinity>>If you impose some limits on z you will get a finite positive answer, >but one that clearly is unbounded:>>In[20]:=>NIntegrate[Boole[y>3x&&y<4-x^2&&z>Out[20]=>2239.58>>In[21]:=>NIntegrate[ Boole[y>3x&&y<4-x^2&&z>Out[21 ]=>20989.6>>Andrzej Sorry, I must have something missing in my previous description. I need to find out the volumn of a 3D object which form by the > equation :> z=x^2 +4 (as bottom surface)> and on the xy plane which bounded by a parabola y=4-x^2 and >y=3x line. How would I use Mathematica to plot out this 3D object or >find out its > volumn with only the equation given? Shz Shz> Andrzej Kozlowski 3/September/2002 04:33pm Mathematica could do this sort of thing if there were a three> dimensional object described by your equations (as boundaries) but> there isn't one. More precisely, the pair of equations {z=x^2 +4,> y=4-x^2} describes a parabola in three space which you can >plot with: g1=ParametricPlot3D[{x, 4 - x^2, x^2 + 4}, {x, -5, 5}] The equation y=3x describes the plane: g2 = ParametricPlot3D[{x, 3x, z}, {x, -5, 5}, {z, -25, 25}]> You can see the two together in> < Show[{g1,g2}] There are clearly two points of intersection. They can be found with:> Solve[{z == x^2 + 4, y == 4 - x^2, y == 3*x}, {x, y, z}]> {{z -> 5, y -> 3, x -> 1}, {z -> 20, y -> -12, x -> -4}} So where is the 3D object whose volume you want to find? Andrzej Kozlowski> Toyama International University> JAPAN Can I use Mathematica to find out the volumn of this 3 dimensional> object from> the equations : z=x^2 +4, y=4-x^2, y=3x> Shz Shz>> >Shz Shz,as Andrzej having said everything about the calculation of the volume, andas still your specification is incomplete, I'm going to show how you canplot at least all you have (communicated).As you said something obout the bottom -- which direction is up? --assuming positive z-direction the I come to the inequalities: y >= 3 x, y <= 4 -x^2, z >= 4 + x^2and arbitrarily, to make the example complete I'll show only parts with z < 20but do not complete the shape.To begin building the graphics we draw the boundary surfaces:[1] the bottom{xmin, xmax} = x /. NSolve[4 - x^2 == 3x, {x}]{-4., 1.}Plot[{4 - x^2, 3x}, {x, xmin, xmax}]{ymin, ymax} = {3 xmin, 4}g3 = Graphics3D[Plot3D[x^2 + 4, {x, xmin, xmax}, {y, ymin, ymax}]][2] the wallsg = ParametricPlot3D[{{x, 4 - x^2, z}, {x, 3x, z}}, {x, xmin, xmax}, {z, 0, 20}]One idea now would be to exploit the Mathematica rendering algorithm, to cutoff the undesired parts of the boundary surfaces:gg = Show[g3, g, PolygonIntersections -> False][Epsilon] = 1. 10^-4g4 = gg /. Line[_] -> {} /. p : Polygon[pts_] :> If[Or @@ (#2 > 4 - #1^2 + [Epsilon] || #2 < 3 #1 - [Epsilon] || #3 < #1^2 + 4 - [Epsilon] &) @@@ pts, {}, p]Show[g4 /. Line[_] -> {}, BoxRatios -> {1, 1, .5}, ViewPoint -> {1.3, -2.4, 5}]But alas! The edges are rather gnawed off. Increasing epsilon doesn't help,unwanted parts will show up and still some polygons are nibbled away.So we have to do it ourselves:We definetag[inequality_][p_] := Block[{x, y, z}, {x, y, z} = p; If[inequality, {p, inside}, {p, outside}]]...tags points whether inside ore outside.sol[equality_, {p1_, _}, {p2_, _}] := Block[{p, d, x, y, z}, {x, y, z} = p = p1 (1 - d) + p2 d; First[Cases[Solve[equality, d], s_ /; NonNegative[d /. s] :> (p /. s)]] ]...computes the cut on the line between to points of opposite sides.sep[equality_][pp1_, pp2_] := If[pp1[[2]] === pp2[[2]], {pp1}, ppx = sol[equality, pp1, pp2]; {pp1, {ppx, pp1[[2]]}, {ppx, pp2[[2]]}}]...effectively cuts a segment crossing the border into two pieces.cutGraphics3D[g_, inequality_] := Module[{equality = Equal @@ inequality}, g /. Polygon[pts_] :> Block[{tagpts = tag[inequality] /@ pts}, With[{t = sep[equality] @@@ Transpose[{tagpts, RotateLeft[tagpts]}]}, With[{newpts = Cases[t, {p_, label_} /; label == inside :> p,2]}, If[Length[newpts] > 2, Polygon[newpts], {}]] ]]]...Polygons crossing a border are cut, only those inside are kept.We first cut the bottom...g3new = Fold[cutGraphics3D, g3, {y >= 3x, y <= 4 - x^2}]...next the sides...gnew = cutGraphics3D[g, z >= 4 + x^2]...and display:Show[gnew, g3new, BoxRatios -> {1, 1, .5}, ViewPoint -> {1.3, -2.4, 5}]Show[gnew, g3new, BoxRatios -> {1, 1, .5}, ViewPoint -> {1.3, -2.4, 15}]Show[gnew, g3new, BoxRatios -> {1, 1, .5}, ViewPoint -> {1.3, -2.4, -2}]Show[gnew, g3new, BoxRatios -> {1, 1, .5}, ViewPoint -> {1.3, 2.5, 2}]Here I have brought a method, I had already posted twice, to a more handyform.To make that reminiscence complete, we could also define a coloring functioncolorGraphics3D[g_, inequality_, colorinside_, coloroutside_] := Module[{equality = Equal @@ inequality}, g /. Polygon[pts_] :> Block[{tagpts = tag[inequality] /@ pts}, With[{t = sep[equality] @@@ Transpose[{tagpts, RotateLeft[tagpts]}]}, {With[{newpts = Cases[t, {p_, label_} /; label == outside :> p, 2]}, If[Length[newpts] > 2,{FaceForm[SurfaceColor[coloroutside]], Polygon[newpts]}, {}]], With[{newpts = Cases[t, {p_, label_} /; label == inside :> p, 2]}, If[Length[newpts] > 2, {FaceForm[SurfaceColor[colorinside]], Polygon[newpts]}, {}]]} ]]]and with...g3D = Graphics3D[Plot3D[E^(-(x^2/2) - y^2/2)/(2*Pi) , {x, -2, 2}, {y, -2,2}]]Show[colorGraphics3D[g3D, With[{e = Take[{x, y, z}, 2] - {0.8, -0.5}}, e.e] <= 1, Hue[0, 0.3, 1], Hue[0.6, 0.3, 1]]]--Hartmut Wolf ==== This result seems to me to be wrong (it's right?):In[1]:=Integrate[Sqrt[a^2*Cos[t]^2 + b^2*Sin[t]^2], {t, 0, 2*Pi}, Assumptions->{Im[a]==0, Im[b]==0}]Out[1]=0....but written in another way, it's right:In[2]:=Integrate[Sqrt[a^2 + (b^2 - a^2)*Sin[t]^2], {t, 0, 2*Pi}, Assumptions -> {Im[a] == 0, Im[b] == 0}]Out[2]=If[a^2/(-a^2 + b^2) >= 0 || b^2/(-a^2 + b^2) <= 0 || Im[a^2/(-a^2 + b^2)] != 0, (4*a^2*Sqrt[-1 + b^2/a^2]*EllipticE[ 1 - b^2/a^2])/Sqrt[-a^2 + b^2], Integrate[Sqrt[a^2 + (-a^2 + b^2)*Sin[t]^2], {t, 0, 2*Pi}]]Someone can explane me this difference?Raf. ==== involved Bessel equations(electro-magnetic field in coaxial) but when I initialize and calculate that2 eqs for the 2nd or 3rd time, Ôcause I changed some parameters, I see astrange Removed[n] instead of the value of n (=0).I use the line Remove[Global`*];at the beginning to clear any variable and I think it could be the problem.So, have I always to quit and restart the program in order to recalculatethat equations?Filippo Sola ==== Can someone provide me more information on how ListIntegrate worksthan what is contained in the version 4 manual ? Is there a websiteperhaps or tutorial ?I would like to know how the beginning and end of a list of {x,y}pairs that is being integrated is dealt with by ListIntegrate.I know a series of polynomials is somehow used but how ? Do theseoverlap ? Are they piecewise continuous ?Are these polynomials available for inspection ? How do they change asa function of k ?optimum k value for a given list ? For any given list, will accuracymonotonically increase with increasing values of k ? Is thereanything in the Option Inspector that could cause unexpected behaivorwith ListIntegrate ?Are there any good rules of thumb or procedures that will help ensurea reasonable answer is produced ?Is there a way of estimating the error or accuracy of integrationperformed by ListIntegrate ? ==== Be sure to note the following and what comes after it near the end of the Help Browser info on ListIntegrate:``This package has been included for compatibility with previous versions of Mathematica. The functionality of this package has been superseded by improvements made to InterpolatingFunction.In other words, ListIntegrate is a dinosaur that you don't need at all!To integrate a list of data with Mathematica, one can proceed in either of two ways: (1) Construct an interpolating function and use NIntegrate (or, better, NIntegrateInterpolatingFunction) on that (which is what ListIntegrate apparently does); or (2) apply a simple routine that implements the trapezoidal rule, Simpson's rule, or maybe some higher order method.Assuming you've chosen to take path #1, you need to realize the following:(a) NIntegrate[Interpolation[data, InterpolationOrder->1][x], {x,a,b}] is equivalent to the trapezoidal rule;(b) NIntegrate[Interpolation[data, InterpolationOrder->2][x], {x,a,b}] is *not* equivalent to Simpson's rule (because of the peculiar way that Interpolation works);(c) Interpolation[data, InterpolationOrder->k] generally does not return a smooth function unless you set InterpolationOrder->n-1, where n is the number of data points, which is the case where a single polynomial of degree n-1 fits the data points.(d) If you want to integrate a smooth interpolant, you can do this: <k][x], {x, 0, 5}, PlotRange -> All], {k, 1, 6}](The last of those plots will give a warning message.)Having said all that, you really should consider path #2 instead. Here are a couple of links to a MathGroup discussion of last July (somehow the thread got split up):http://library.wolfram.com/mathgroup/archive/2002/Jul/ msg00490.htmlhttp://library.wolfram.com/mathgroup/archive/ 2002/Jul/msg00519.htmlI hope all this helps some.----Selwyn Hollis > Can someone provide me more information on how ListIntegrate works > than what is contained in the version 4 manual ? Is there a website > perhaps or tutorial ? > I would like to know how the beginning and end of a list of {x,y} > pairs that is being integrated is dealt with by ListIntegrate. > I know a series of polynomials is somehow used but how ? Do these > overlap ? Are they piecewise continuous ? > Are these polynomials available for inspection ? How do they change as > a function of k ? > optimum k value for a given list ? For any given list, will accuracy > monotonically increase with increasing values of k ? Is there > anything in the Option Inspector that could cause unexpected behaivor > with ListIntegrate ? > Are there any good rules of thumb or procedures that will help ensure > a reasonable answer is produced ? > Is there a way of estimating the error or accuracy of integration > performed by ListIntegrate ? > > ==== I timed Daniel's three solutions and Gary's one (plus a couple of my owna little later):perps1[v_] := If[v[[1]] == v[[2]] == 0, {{1, 0, 0}, {0, 1, 0}}, {{v[[2]], -v[[ 1]], 0}, Cross[v, {v[[2]], -v[[1]], 0}]}]perps2[v_] := With[{vecs = NullSpace[{v}]}, {vecs[[ 1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]}]perps2C = Compile[{{v, _Real, 1}}, Module[{vecs = NullSpace[{v}]}, { vecs[[1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]}]]helzer[v : {a1_, a2_, a3_}] := With[{w = First[Sort[{{a2, -a1, 0}, {a3, 0, -a1}, {0, a3, -a2}}, OrderedQ[{Plus @@ Abs[#2], Plus @@ Abs[#1]}] &]]}, {w, Cross[v, w]}]vecs = Table[Random[], {10000}, {3}]; Timing[perps1 /@ vecs; ]Timing[perps2 /@ vecs; ]Timing[perps2c /@ vecs; ]Timing[helzer /@ vecs; ]{1.7350000000000012*Second, Null}{0.5619999999999994*Second, Null}{0.219 Second, Null}{2.7349999999999994*Second, Null}I made a small change to Daniel's perps1, and got a solution as fast asperps2c, WITHOUT compiling. Compiling tripled the speed again, sotreatC is the fastest solution I've seen so far.treat[{a_, b_, c_}] := If[a == b == 0, {{1, 0, 0}, {0, 1, 0}}, {{b, -a, 0}, {a*c, b*c, -a^2 -b^2}}]treatC = Compile[{{v, _Real, 1}}, If[v[[1]] == v[[2]] == 0, {{1, 0, 0}, {0, 1, 0}}, {{v[[2]], -v[[1]], 0}, {v[[1]]*v[[3]], v[[2]]*v[[3]], -v[[1]]^2 - v[[2]]^2}}]]vecs = Table[Random[], {10000}, {3}]; Timing[perps2c /@ vecs; ]Timing[helzer /@ vecs; ]Timing[treat /@ vecs; ]Timing[treatC /@ vecs;]{0.2190000000000083*Second, Null}{2.7339999999999947*Second, Null}{0.25*Second, Null}{0.07800000000000296*Second, Null}None of these solutions reliably return normalized vectors.Bobby Treat-----Original Message-----> pose the problem to MathGroup. Who has the most elegant Mathematica> routine... OrthogonalUnitVectors::usage = OrthogonalUnitVectors[v:{_,_,_}] willreturn> two unit vectors orthogonal to each other and to v. You can assume that v is nonzero. David Park> djmp@earthlink.net> http://home.earthlink.net/~djmp/Some possibilities:perps1[v_] := If [v[[1]]==v[[2]]==0, {{1,0,0},{0,1,0}}, {{v[[2]],-v[[1]],0}, Cross[v,{v[[2]],-v[[1]],0}]} ]perps2[v_] := With[{vecs=NullSpace[{v}]}, {vecs[[1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]} ]This appears to be 2-3 times faster than perps1 for vectors of machinereals. I get another factor of 2 using Compile, which is appropriate fore.g. graphics use.perps2C = Compile[{{v,_Real,1}}, Module[{vecs=NullSpace[{v}]}, {vecs[[1]], vecs[[2]] - (vecs[[2]].vecs[[1]])*vecs[[1]]} ]]In[61]:= vecs = Table[Random[], {10000}, {3}];In[62]:= Timing[p2 = Map[perps1C,vecs];]Out[62]= {0.49 Second, Null}This is on a 1.5 GHz processor.Daniel LichtblauWolfram Research ==== >I was wondering if there is a method for resizing Raster graphics>(resizing the actual matrix of pixels, not just the display size). I>am processing a large number of JPEG images, and we sometimes need to>reduce the image size to allow data processing algorithms to function>without running out of memory. In the past we simply used a program>such as Photoshop to resize them before importing them into>Mathematica. Due to the number of images we are processing now this>is very inconvenient and it would be very useful if there was a method>for accomplishing it in Mathematica, but I can't find one. I also>thought about doing something simple like sampling every few pixels or>averaging, but I thought there might be a method with more efficacy>than this. I also tried exporting the graphics with the Export command>as new JPEGs and manipulating the ImageResolution and ImageSize>options but this seemed to have no effect. Any help would be much>appreciated.>Aaron UrbasIn 4.2, the following should work (as long as you define newsize).in = Import[file.jpg];Export[newfile.jpg, in, ImageSize->newsize]In 4.1 and earlier, this will not work because an optimization interferes with the ImageSize and rasters are written out with the raster size, not the ImageSize. There is a ConversionOption to force the behavior you want.in = Import[file.jpg];Export[newfile.jpg, Show[in, ImageSize->newsize], ConversionOptions->{RasterExport->Graphics}]-Dale ==== Borrowing liberally from Daniel, I like the following:ClearAll[sumBin, sumBinOdd, sumBinEven, index]sumBinOdd = Sum[Binomial[2index + 1, k]*p^(2index + 1 - k)*(1 - p)^k, { k, 0, index - 1}];sumBinEven = Sum[Binomial[2index, k]* p^(2index - k)*(1 - p)^k, {k, 0, index - 1}];sumBin[n_, Odd] = sumBinOdd /. {index -> (n - 1)/2};sumBin[n_, Even] = sumBinEven /. {index -> n/2};sumBin[n_?EvenQ] = sumBin[n, Even];sumBin[n_?OddQ] = sumBin[n, Odd];It allows you to see the solution symbolically for both odd and even n,and also to calculate it when n is a known integer. We also have theopportunity, for instance, to assume that 3x is even and calculatesumBin[3x, Even]p^(3*x)*((1/p)^(3*x) - ((-4 + 4/p)^((3*x)/2)* Gamma[1/2 + (3*x)/2]*Hypergeometric2F1[1, -((3*x)/2), 1 + (3*x)/2, (-1 + p)/p])/(Sqrt[Pi]*Gamma[1 +(3*x)/2]))or assume 3x is odd and calculatesumBin[3*x, Odd]p^(3*x)*((1/p)^(3*x) - (2^(3*x)*(-1 + 1/p)^((1/2)*(-1 + 3*x))* Gamma[3/2 + (1/2)*(-1 + 3*x)]*Hypergeometric2F1[1, -1 + (1/2)*(1 - 3*x), 1 + (1/2)*(-1 + 3*x), (-1 +p)/p])/(Sqrt[Pi]* Gamma[2 + (1/2)*(-1 + 3*x)]))Bobby Treat-----Original Message-----> when F[n,p] is an approximation to the sum and the R is the remainingerror. Constantine. >In a message dated 8/28/02 4:44:13 AM, celster@cs.technion.ac.il>>I'm looking for a way of finding the approximation for partitialbinomial>>sum.>>I'll be pleasant for any hint..> [...]> Office: Taub 411You can get a closed form in terms of special functions if you splitinto two cases depending on whether n is even or odd.In[39]:= n = 2*m;In[40]:= InputForm[Sum[Binomial[n,k]*p^(n-k)*(1-p)^k, {k,0,m-1}]]Out[40]//InputForm= p^(2*m)*((p^(-1))^(2*m) - ((-4 + 4/p)^m*Gamma[1/2 + m]* Hypergeometric2F1[1, -m, 1 + m, (-1 + p)/p])/(Sqrt[Pi]*Gamma[1 +m]))In[41]:= n = 2*m+1;In[42]:= InputForm[Sum[Binomial[n,k]*p^(n-k)*(1-p)^k, {k,0,m-1}]]Out[42]//InputForm= p^(1 + 2*m)*((p^(-1))^(1 + 2*m) - (2^(1 + 2*m)*(-1 + p^(-1))^m*Gamma[3/2+ m]* Hypergeometric2F1[1, -1 - m, 1 + m, (-1 + p)/p])/(Sqrt[Pi]*Gamma[2 +m]))Daniel LichtblauWolfram Research ==== I've found that Mathematica 4.2 takes too long to finish to generate firsttime help browser. Is this normal? How many time should it take?Calimero ==== Please excuse me if this has been discussed before; I haven't been monitoring the group very much.A long time ago I posted to this group as well as to Wolfram about the errors in many or all functions involving Fourier transforms when values of FourierParameters other than the default values are used. The group responded that indeed there was a problem so I'm wondering if it was ever fixed. My version (then and now) is 4.0.1 for Macintosh.Jerry ==== Could you help me?I'm to solve heat conductivity equation (Laplas equation) - partial differential equation. Are there any ready-to-use packages that will help me do the job?Standard function DSolve cannot! And so do functions from package Calculus`DSolveIntegrals`.More info: the space where the equation is to be solved is cylindre (not infinite). Andrew. ==== When I multiply an expression with 0 it is giving 0.expression whichis creating problem in the further calculation.eg.In[1]:=func1[r_]:=Exp[-r/2]In[2]:coeff[[1,1 ]]=0.0;;;In[34]:=func1[r]coeff[[1,1]]Out[34]:=0.Exp[-r/2]I want anything to be multiply by zero must be zero. How can I do that? Your suggestion will be highly appreciated. Raj ==== > When I multiply an expression with 0 it is giving 0.expression which> is creating problem in the further calculation.No. Your trouble comes when you multiply an expression by theßoating-point number 0.0, rather than by the precise symbolic 0 (having nodecimal point). Merely make the coefficient 0 (rather than 0.0) andeverything should work as you wish.David> eg. In[1]:=func1[r_]:=Exp[-r/2]> In[2]:coeff[[1,1]]=0.0;> ;> ;> In[34]:=func1[r]coeff[[1,1]] Out[34]:=0.Exp[-r/2] I want anything to be multiply by zero must be zero. How can I do that?-- -------------------- http://NewsReader.Com/ -------------------- Usenet Newsgroup ServiceReply-To: kuska@informatik.uni-leipzig.de ==== Unprotect[Times]Times[0 ., __] := 0Protect[Times] Jens > When I multiply an expression with 0 it is giving 0.expression which> is creating problem in the further calculation. eg. In[1]:=func1[r_]:=Exp[-r/2]> In[2]:coeff[[1,1]]=0.0;> ;> ;> In[34]:=func1[r]coeff[[1,1]] Out[34]:=0.Exp[-r/2] I want anything to be multiply by zero must be zero. How can I do that?> Your suggestion will be highly appreciated. > Raj ==== Now I'm trying to calculate this formula:Delta[eq_, x_]:=Coefficient[eq, x]^2 - 4 Coefficient[eq, x^2] Coefficient[eq, x^0]eq has this form a x^2 + b x + cBut there is a problem with the x^0 coefficient!How can I overcome that?CeZaRReply-To: kuska@informatik.uni-leipzig.de ==== andDelta[eq_, x_]:=Coefficient[eq, x]^2 - 4 Coefficient[eq, x,2]Coefficient[eq, x,0]does what you want. Because the Help-Browser say:Coefficient[expr, form, 0] picks out terms that are not proportional toform. Jens > Now I'm trying to calculate this formula: Delta[eq_, x_]:=Coefficient[eq, x]^2 - 4 Coefficient[eq, x^2] Coefficient[eq, x^0] eq has this form a x^2 + b x + c But there is a problem with the x^0 coefficient!> How can I overcome that? > CeZaR ==== I need to fill the space between two contour lines, C1 and C2, withred color, and leave the other place white. What trick I have to use?Any suggestion and advice will be appreciated.Jun LinReply-To: kuska@informatik.uni-leipzig.de ==== gr = ContourPlot[x^2 + y^2, {x, -1, 1}, {y, -1, 1}, Contours -> {0.2, 0.4}, ColorFunction -> (If[# >= 0.2 && # <= 0.4, RGBColor[1, 0, 0], RGBColor[1, 1, 1]] &), ColorFunctionScaling -> False ] Jens I need to fill the space between two contour lines, C1 and C2, with> red color, and leave the other place white. What trick I have to use?> Any suggestion and advice will be appreciated. Jun Lin ==== I am trying to find an example that will demonstrate the difference between$PrePrint and $Post. I found an old thread in this news group where auser wanted to display all matrices using MatrixForm. Some users suggestedthe following: In[1]:= $Post=(#/.mtrx_?MatrixQ:>MatrixForm[mtrx]&);Then Dave Withoff said it's better to assign this to $PrePrint since theobjective here is to adjust the display rather than the result of thecalculation. With the assignment to $Post you could, for example, getunexpected results from calculations using %, since matrices will be wrappedin MatrixForm.--------However, if we use $Post above, the next input will compute the inversethe matrix. I did verify that Inverse can't take a matrix wrapped inMatrixForm. Can somebody give an example where doing this with $PrePrintinstead of $Post gives a different result. In[2]:= m={{2,3},{0,1}}; Inverse[%] Out[3]= (* Inverse of (m) in MatrixForm, not shown. *)------ Ted Ersek Get Mathematica tips, tricks from http://www.verbeia.com/mathematica/tips/Tricks.htmlReply-To : kuska@informatik.uni-leipzig.de ==== yo can just tryIn[]:=$Post = (# /. mtrx_?MatrixQ :> AnyHead[mtrx] &);In[]:=m = {{2, 3}, {0, 1}}In[]:=q=%;In[]:=Head[q]andIn[]:=$PrePost = (# /. mtrx_?MatrixQ :> AnyHead[mtrx] &);In[]:=m = {{2, 3}, {0, 1}}In[]:=q=%;In[]:=Head[q]But you are right -- the behaviour of MatrixForm[] in your exampleis strange. Jens I am trying to find an example that will demonstrate the difference between> $PrePrint and $Post. I found an old thread in this news group where a> user wanted to display all matrices using MatrixForm. Some users suggested> the following: In[1]:= $Post=(#/.mtrx_?MatrixQ:>MatrixForm[mtrx]&); Then Dave Withoff said it's better to assign this to $PrePrint since the> objective here is to adjust the display rather than the result of the> calculation. With the assignment to $Post you could, for example, get> unexpected results from calculations using %, since matrices will be wrapped> in MatrixForm.> -------- However, if we use $Post above, the next input will compute the inverse> the matrix. I did verify that Inverse can't take a matrix wrapped in> MatrixForm. Can somebody give an example where doing this with $PrePrint> instead of $Post gives a different result. In[2]:= m={{2,3},{0,1}};> Inverse[%] Out[3]= (* Inverse of (m) in MatrixForm, not shown. *) ------> Ted Ersek> Get Mathematica tips, tricks from> http://www.verbeia.com/mathematica/tips/Tricks.html ==== myArray is a list of triplets which look like {n,k,p}, where n and kare integers and p is real number. repetitions in n and k are allowedbut no two triplets ar the same.I define a function, myfunc[n_,pthreshold_] as follows.myfunc[n_,pthreshold_] := Module[{t},t=Max[Cases[myArray, {n, k_, p_} /; (p >= pthreshold) -> p, 2]];If[t>=0,t,-1]]It prints the largest k for which p is greater than pthreshold for agiven n and-1 if there is no such k. The function does its job fine. My problem is as follows. I would like to plot a collection of plots,myfunc[n,p], for 1<=n<=21 with respect to p. unfortunately thefollowing does not work as desired.Plot[Table[myfunc[n,p],{n,1,21}],{p,0.0,1.0}] because the Table gets evaluated to numeric value (-1) before Plot isinvoked. ==== I made a mistake in the definition of the question. The function is correctly defined as myfunc[n_,pthreshold_] := Module[{t}, t=Max[Cases[myArray, {n, k_, p_} /; (p >= pthreshold) -> k, 2]]; If[t>=0,t,-1]]Please note that inserting Evaluate before Table does not workcorrectly either, because evaluating myfunc[n,p] when p is not anumber, gives -1.> Plot[Table[myfunc[n,p],{n,1,21}],{p,0.0,1.0}] because the Table gets evaluated to numeric value (-1) before Plot is> invoked.Reply-To: kuska@informatik.uni-leipzig.de ==== could you supply a complete example next time ?Try myfunc[] with a more restrictive pattern and it works:myArray = Flatten[Table[{n, 1, Random[]}, {15}, {n, 1, 21}], 1];myfunc[n_, pthreshold_?NumericQ] := Module[{t}, t = Max[Cases[myArray, {n, k_, p_} /; (p >= pthreshold) -> p, 2]]; If[t >= 0, t, -1]]Plot[Evaluate[Table[myfunc[n, p], {n, 1, 21}]], {p, 0.0, 1.0}] Jens myArray is a list of triplets which look like {n,k,p}, where n and k> are integers and p is real number. repetitions in n and k are allowed> but no two triplets ar the same. I define a function, myfunc[n_,pthreshold_] as follows. myfunc[n_,pthreshold_] := Module[{t},> t=Max[Cases[myArray, {n, k_, p_} /; (p >= pthreshold) -> p, 2]];> If[t>=0,t,-1]] It prints the largest k for which p is greater than pthreshold for a> given n and> -1 if there is no such k. The function does its job fine. My problem is as follows. I would like to plot a collection of plots,> myfunc[n,p], for 1<=n<=21 with respect to p. unfortunately the> following does not work as desired. Plot[Table[myfunc[n,p],{n,1,21}],{p,0.0,1.0}] because the Table gets evaluated to numeric value (-1) before Plot is> invoked. ==== For my PDE class we have been calculating Fourier transforms. The instructorarrived today with a printout of two plots of a certain Fourier transform,done with a different CAS. The first plot was to 30 terms, the second was to120 terms. Curious, I translated the functions into Mathematica (4.0 on Windows2000on a PIII 700) to see how much time this required to process. I wasStaggered at how much time it took. Here's the code:L = 2;f[x_] := UnitStep[x - 1];b[n_] := (2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}];FS[N_, x_] := Sum[b[n]*Sin[n*Pi*x/L], {n, 1, N}];Timing[Plot[FS[30, x], {x, 0, 2}]]Out[23]={419.713 Second, [SkeletonIndicator]Graphics[SkeletonIndicator]}In this case the number of terms is 30.The time required per number of terms seems to fit the following polynomial:y = 0.3926x^2 + 2.2379xThis is a large amount of time. I understand that the code is not optimized,and was more or less copied from the code in the other CAS, but is this areasonable amount of time, or is something going wrong? I don't use Mathematicabecause of the speed, but should it be this slow?Just curious,Steve Story ==== > I don't use Mathematica> because of the speed, but should it be this slow?Mathematica has become more competitive in the speed department in recentyears. See for example the attached comparison (not sent to newsgroup) byStephan Steinhaus (steinhaus-net.de). So when Mathematica takes a very longtime, you should investigate. In this case inserting Evaluate[] in twoplacesIn[91]:=b[n_] := Evaluate[(2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}]];....In[104]:=Timing[Plot[Evaluate[FS[120, x]], {x, 0, 2}]]Out[104]={0.18 Second,[SkeletonIndicator]Graphics[SkeletonIndicator]} speeds the process enormously (18 milliseconds to plot 120 terms on myfeeble old 500MHz PowerBook).Why was it so slow before? When I switch from an ordinary numerical languageto Mathematica, I enter into an implicit bargain withMathematica: the software will go the extra mile to get me a good answer,including (1) using exra precision (sometimes without being asked) and (2)carrying around unevaluated mathematical expressions (usually without beingasked) that could possibly be evaluated more appropriately at a later time.Most tools cannot do either of these things, so I don't have to worry aboutit, except for the bad answers that result now and then. But I need to takecare that Mathematica does not burden itself unnecessarily. That's my sideof the bargain.Number (2) is the issue here. Your definition of b[n] is written so thatMathematica analytically evaluates b separately for each n. But you know inthis case that the integration can be done safely once for all n. So do it!The huge difference, though, comes from pre-evaluating the argument to Plot.Read the on-line help! You should pre-evaluate where possible. In somecases, the most common of which involve branching within the definition offunction to plot, you cannot pre-evaluate so, in keeping with the bargain,Mathematica goes the extra mile and holds back just in case. You need tosteer it into the shortcut when it's OK.Hope this helps,Tom Burton-- Reply-To: kuska@informatik.uni-leipzig.de ==== with your code you compute the expansion coefficents everytime when FS[] is evaluated. Store the values for b[n] withL = 2;f[x_] := UnitStep[x - 1];b[n_] := b[n] = (2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}];FS[N_, x_] := Sum[b[n]*Sin[n*Pi*x/L], {n, 1, N}];Timing[Plot[FS[30, x], {x, 0, 2}]]and you need only a 1-3 seconds (depending on your machine) Jens For my PDE class we have been calculating Fourier transforms. The instructor> arrived today with a printout of two plots of a certain Fourier transform,> done with a different CAS. The first plot was to 30 terms, the second was to> 120 terms.> Curious, I translated the functions into Mathematica (4.0 on Windows2000> on a PIII 700) to see how much time this required to process. I was> Staggered at how much time it took. Here's the code: L = 2;> f[x_] := UnitStep[x - 1];> b[n_] := (2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*x/L], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] Out[23]=> {419.713 Second, [SkeletonIndicator]Graphics[SkeletonIndicator]} In this case the number of terms is 30. The time required per number of terms seems to fit the following polynomial: y = 0.3926x^2 + 2.2379x This is a large amount of time. I understand that the code is not optimized,> and was more or less copied from the code in the other CAS, but is this a> reasonable amount of time, or is something going wrong? I don't use Mathematica> because of the speed, but should it be this slow? Just curious, Steve Story ==== I would like to use mathematica type papers for my math courses, butI'm having trouble formatting documents. Despite searching, I've beenunable to find a complete guide to word processing with mathematica.Does anyone know where such a document could be found? ==== Look at in the Help browser or in the Mathematica Book Style Sheet Also, you have the documentation included for package Author Tools(only in Mathematica 4.2).Guillermo Sanchez > I would like to use mathematica type papers for my math courses, but> I'm having trouble formatting documents. Despite searching, I've been> unable to find a complete guide to word processing with mathematica.> Does anyone know where such a document could be found? ==== Kenny,Sympathy but no solution.I too have been trying to use Mathematica (v4.2 most recently) to typemaths papers and the like but I'm not ready to ditch LaTeX yet. Thereare just too many cases where I cannot figure out how to achieve what Iwant in Mathematica, things like:- left brackets spanning multiple lines for defining hybrid functions;- vertical alignment of equals signs in multi-line equations orderivations;- setting typefaces in tables of material.I figure most of this is do-able, but I don't have the time, orpatience, to spend too much time on it. So, I'll be the first customerwhen you write the guide to math publishing in Mathematica - I just hopeyou won't have to use LaTeX to write it.Mark Westwood I would like to use mathematica type papers for my math courses, but> I'm having trouble formatting documents. Despite searching, I've been> unable