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 411 You can get a closed form in terms of special functions if you split into 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])) ==== 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. > > Constantine Elster Computer Science Dept. Technion I.I.T. Office: Taub 411 ==== Constantine, If you break your problem up into two cases, even n and odd n, then Mathematica can sum up your problem and get results, albeit with hypergeometric functions. Consider the following (make sure you look at this with 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 p 1 - -------------------------------------------------------------------- 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 3 p - 1 2 (1 - p) p Gamma[n + -] Hypergeometric2F1[1, -n - 1, n + 1, -----] 2 p 1 - ------------------------------------------------------------------------ -------- Sqrt[Pi] Gamma[n + 2] Is this what you were looking for? Carl Woll Physics Dept U 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 remaining error. Constantine. >I'm looking for a way of finding the approximation for partitial binomial >sum. Use the standard add-on package Statistics`NonlinearFit` to do a NonlinearFit to whatever model you want to use for the approximation. ==== >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 Table l == Table[eq, {m, -30, 29}] True To Plot, use Evaluate Plot[Evaluate[l], {x, -50, 50}]; ==== >What is a efficient way to use output expression to define new function. Use the menu command Input/Copy Output from Above or use % (Out) ==== >-----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]:= ><<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-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&&z2239.58 In[21]:= >NIntegrate[Boole[y>3x&&y<4-x^2&&z20989.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 bot 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? >> >>> 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 >>> > >> = 3 x, y <= 4 -x^2, z >= 4 + x^2 and arbitrarily, to make the example complete I'll show only parts with z < 20 but do not complete the shape. To begin building the graphics we draw the boundary surfaces: [1] the bot {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 walls g = 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 cut off the undesired parts of the boundary surfaces: gg = Show[g3, g, PolygonIntersections -> False] [Epsilon] = 1. 10^-4 g4 = 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 define tag[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 bot... 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 handy form. To make that reminiscence complete, we could also define a coloring function colorGraphics3D[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 that 2 eqs for the 2nd or 3rd time, 'cause I changed some parameters, I see a strange 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 recalculate that equations? Filippo Sola ==== 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 ? ==== 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.html http://library.wolfram.com/mathgroup/archive/2002/Jul/msg00519.html I hope all this helps some. > 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 own a 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 as perps2c, WITHOUT compiling. Compiling tripled the speed again, so treatC 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. -----Original Message----- > pose the problem to MathGroup. Who has the most elegant Mathematica > routine... OrthogonalUnitVectors::usage = OrthogonalUnitVectors[v:{_,_,_}] will return > two unit vectors orthogonal to each other and to v. You can assume that v is nonzero. 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 machine reals. I get another factor of 2 using Compile, which is appropriate for e.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 Lichtblau Wolfram 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 Urbas In 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 the opportunity, for instance, to assume that 3x is even and calculate sumBin[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 calculate sumBin[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)])) -----Original Message----- > when F[n,p] is an approximation to the sum and the R is the remaining error. 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 partitial binomial >sum. >I'll be pleasant for any hint.. [...] > Office: Taub 411 You can get a closed form in terms of special functions if you split into 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 Lichtblau Wolfram Research ==== I've found that Mathematica 4.2 takes too long to finish to generate first time 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 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 ==== > 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 floating-point number 0.0, rather than by the precise symbolic 0 (having no decimal point). Merely make the coefficient 0 (rather than 0.0) and everything should work as you wish. > 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 Service Reply-To: kuska@informatik.uni-leipzig.de ==== Unprotect[Times] Times[0., __] := 0 Protect[Times] > 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 + c But there is a problem with the x^0 coefficient! How can I overcome that? CeZaR Reply-To: kuska@informatik.uni-leipzig.de ==== and Delta[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 to form. > 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, with red color, and leave the other place white. What trick I have to use? Any suggestion and advice will be appreciated. Jun Lin Reply-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 ] 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 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 Reply-To: kuska@informatik.uni-leipzig.de ==== yo can just try In[]:=$Post = (# /. mtrx_?MatrixQ :> AnyHead[mtrx] &); In[]:=m = {{2, 3}, {0, 1}} In[]:=q=%; In[]:=Head[q] and In[]:=$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 example is strange. 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 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. ==== 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 work correctly either, because evaluating myfunc[n,p] when p is not a number, 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}] 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 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 don't use Mathematica > because of the speed, but should it be this slow? Mathematica has become more competitive in the speed department in recent years. See for example the attached comparison (not sent to newsgroup) by Stephan Steinhaus (steinhaus-net.de). So when Mathematica takes a very long time, you should investigate. In this case inserting Evaluate[] in two places In[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 my feeble old 500MHz PowerBook). Why was it so slow before? When I switch from an ordinary numerical language to Mathematica, I enter into an implicit bargain with Mathematica: 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 being asked) 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 about it, except for the bad answers that result now and then. But I need to take care that Mathematica does not burden itself unnecessarily. That's my side of the bargain. Number (2) is the issue here. Your definition of b[n] is written so that Mathematica analytically evaluates b separately for each n. But you know in this 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 some cases, the most common of which involve branching within the definition of function 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 to steer it into the shortcut when it's OK. Hope this helps, Burton -- Reply-To: kuska@informatik.uni-leipzig.de ==== with your code you compute the expansion coefficents every time when FS[] is evaluated. Store the values for b[n] with L = 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) 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, 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? ==== 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 type maths papers and the like but I'm not ready to ditch LaTeX yet. There are just too many cases where I cannot figure out how to achieve what I want in Mathematica, things like: - left brackets spanning multiple lines for defining hybrid functions; - vertical alignment of equals signs in multi-line equations or derivations; - setting typefaces in tables of material. I figure most of this is do-able, but I don't have the time, or patience, to spend too much time on it. So, I'll be the first cuser when you write the guide to math publishing in Mathematica - I just hope you 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 to find a complete guide to word processing with mathematica. > Does anyone know where such a document could be found? ==== OrthogonalUnitVectors that I sent a few minutes ago doesn't do what I thought it would. After making either definition below I get the following: In[2]:= s1=OrthogonalUnitVectors[{1,0,1/2,1,0},{0,1,-1,1/2,2}] Out[2]= {{0, -2/Sqrt[5], 0, 0, 1/Sqrt[5]}, {-2/3, -1/3, 0, 2/3, 0}, {-1/3, 2/3, 2/3, 0, 0}} The dot products below aren't zero, so the vectors aren't orthogonal. What went wrong? In[3]:= Part[s1,1].Part[s1,2] Out[3]= 2/(3*Sqrt[5]) In[4]:= Part[s1,1].Part[s1,3] Out[3]= -4/(3*Sqrt[5]) -------------- > Hugh Goyder and Park gave a most elegant function to find two > vectors that are orthogonal to one vector in 3D. The key to coming up > with the elegant solution is an understanding of Mathematica's NullSpace > function. We can easily make the version from Hugh and much more > general with the version below. > ------------- > OrthogonalUnitVectors[vect__?VectorQ]:= > #/Sqrt[#.#]&/@NullSpace[{vect}] ------------- > The version above will give a set of unit orthogonal vectors if given any > number of vectors in any dimension. > So besides giving it a 3D vector we can give it the following: > OrthogonalUnitVectors[{2,1,0,-1,1}] > or > OrthogonalUnitVectors[{0,1,0,1/2,1},{1,0,-1,1/2}] ------------ > But the short version above isn't very robust. > (1) Clear[x,y,z];NullSpace[{{x,y,z}}] > returns two vectors orthogonal to {x,y,z}, but the two vectors > NullSpace returns aren't orthogonal to each other. > So (OrthogonalUnitVectors) should only work with numeric vectors. (2) We should ensure all the vectors have the same dimension and length 1. I give a less concise version below that corrects these problems. > ------------ OrthogonalUnitVectors[vect__?(VectorQ[#,NumericQ]&)]/; > (SameQ@@Length/@{vect})&&(Length[First[{vect}]]>1):= > #/Sqrt[#.#]&/@NullSpace[{vect}] -------------- > Ted Ersek > Get Mathematica tips, tricks from > http://www.verbeia.com/mathematica/tips/Tricks.html > ==== Daniel Lichtblau has pointed out that NullSpace does not generally give orthogonal vectors. Therefore the routines that depended upon that were in error. He says that it does give orthogonal vectors when the input vector contains approximate numbers. For graphical purposes this will be good enough for me. Therefore I modify Ted's routine to OrthogonalUnitVectors[vect__?(VectorQ[#, NumericQ] &)] /; (SameQ @@ Length /@ {vect}) && (Length[First[{vect}]] > 1) := #/Sqrt[#.#] & /@ NullSpace[{vect}// N] and the short version for 3D vectors OrthogonalUnitVectors[v : {_, _, _}] := #/Sqrt[#.#] & /@ NullSpace[{v//N}] For exact vectors I might use for 3D OrthogonalUnitVectors[v : {_, _, _}] := #/Sqrt[#.#] & /@ {temp = First[NullSpace[{v}]], v[Cross]temp} I'm still looking for something that is easy to remember. Park djmp@earthlink.net http://home.earthlink.net/~djmp/ The version above will give a set of unit orthogonal vectors if given any number of vectors in any dimension. So besides giving it a 3D vector we can give it the following: OrthogonalUnitVectors[{2,1,0,-1,1}] or OrthogonalUnitVectors[{0,1,0,1/2,1},{1,0,-1,1/2}] ------------ But the short version above isn't very robust. (1) Clear[x,y,z];NullSpace[{{x,y,z}}] returns two vectors orthogonal to {x,y,z}, but the two vectors NullSpace returns aren't orthogonal to each other. So (OrthogonalUnitVectors) should only work with numeric vectors. (2) We should ensure all the vectors have the same dimension and length >1. I give a less concise version below that corrects these problems. ------------ OrthogonalUnitVectors[vect__?(VectorQ[#,NumericQ]&)]/; (SameQ@@Length/@{vect})&&(Length[First[{vect}]]>1):= #/Sqrt[#.#]&/@NullSpace[{vect}] -------------- Ted Ersek Get Mathematica tips, tricks from http://www.verbeia.com/mathematica/tips/Tricks.html ==== I am unsure if my messages are making it through to the usegroup so I will try again. I have several differential equations I would like to plot in Mathematica. I would like to plot the Slope Fields of them though. Can anyone lead me in the right direction? I can solve the equations trivially but I want to display the slope fields. An example follows : y' + 2y = 3 -John ==== John, You can do it from scratch with PlotVectorField from the Graphics`PlotField` package, but you make make your life easier (and get prettier plots) by using my DEGraphics package, which can be found at MathSource (it's part of the DiffEqs suite of packages) or here: http://www.math.armstrong.edu/faculty/hollis/mmade/DiffEqs --- > I am unsure if my messages are making it through to the usegroup so I > will try again. I have several differential equations I would like to > plot in Mathematica. I would like to plot the Slope Fields of them > though. Can anyone lead me in the right direction? I can solve the > equations trivially but I want to display the slope fields. An example > follows : y' + 2y = 3 -John > Reply-To: kuska@informatik.uni-leipzig.de ==== Needs[Graphics`PlotField`] PlotVectorField[{x, 3 - 2 y}, {x, 0, 4}, {y, -1, 4}, Axes -> True] ??? I am unsure if my messages are making it through to the usegroup so I > will try again. I have several differential equations I would like to > plot in Mathematica. I would like to plot the Slope Fields of them > though. Can anyone lead me in the right direction? I can solve the > equations trivially but I want to display the slope fields. An example > follows : y' + 2y = 3 > -John ==== You can proceed like this : delta[poly2_, var_] := Coefficient[poly2, var, 1]^2 - 4*Coefficient[poly2, var, 0]*Coefficient[poly2, var, 2] Meilleures salutations Florian Jaccard -----Message d'origine----- Envoy.8e : ven., 6. septembre 2002 09:17 Ë : mathgroup@smc.vnet.net Objet : Coefficient problem 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 ==== This works: eq = a x^2 + b x + c discriminant[eq_, x_] := Coefficient[eq, x]^2 - 4 Coefficient[eq, x, 2] Coefficient[eq, x, 0] discriminant[eq, x] x^0 is reduced to 1 and the Coefficient of 1 doesn't make sense to Mathematica, because it depends on what the variable is (a, b, c, or x?). So, the other form of the Coefficient call is needed. I used it for the second power too, but that wasn't necessary. I think that form is best, though, since it allows no ambiguity. I renamed the function because that's what the quantity is often called, for a quadratic. -----Original Message----- CeZaR ==== >-----Original Message----- >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. > The problem with your calculation indented is (1) to evaluate Table (within Plot), which gives {myfunc[1, p], ..., myfunc[21, p]} but (2) *then* to prevent further evaluation of myfunc[.., p] to -1. There are several tricks to do so, e.g. ... Hold[Plot[toPlot, {p, 0, 1}]] /. toPlot -> Table[headPlaceholder[n, p], {n, 3}] /. headPlaceholder -> myfunc // ReleaseHold ...but in your case just simply prevent evaluation of myfunc[.., p] with p being a Symbol by defining: myfunc2[n_, pthreshold_?NumericQ] := ... now... Plot[Evaluate[Table[myfunc2[n, p], {n, 21}]], {p, 0, 1}] ...works as expected. BTW, what do you want to read off the plot, that you didn't know from Table[myfunc[n, 0], {n, 21}] ? -- Hartmut Wolf ==== >-----Original Message----- >Sent: Friday, September 06, 2002 9:17 AM 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 > Raj, this simply is, because 0. Exp[-r/2] is not always zero! func1[-Infinity]coeff[[1, 1]] Infinity::indet: Indeterminate expression 0. Infinity encountered. Out[90]= Indeterminate -- Hartmut Wolf ==== I have a mathematica notebook showing using the alternating direction implicit method for solving heat conduction in 2D via finite difference methods. Moving from 2D to 3D is pretty simple although I can't find my notebook on this anymore :( . The 3D method is sometimes called Brian's method. http://mid-ohio.mse.berkeley.edu/scott/projects/index.html Scott > 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. > ==== I am working with Mathematica 4.0, and I would like pass my data to do one graphic in Origin 5.0. How can I do this in more simple manner? ==== > 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. Over the years I have alighted upon the following scheme. Suppose I want to color only between contour levels 1 and 2: In[84]:=Needs[Graphics`Colors`] In[88]:=ContourPlot[x, {x, -1, 3}, {y, 0, 1}, Contours -> {0, 1, 2, 3, 4}, ColorFunction -> (If[1 < # < 2, Red, White] & ), ColorFunctionScaling -> False]; 1. Specify the specific contour levels instead of specifying only the count. 2. Disable ColorFunctionScaling so the argument to ColorFunction corresponds to the contour levels. Hope this helps, Burton ==== Park replied with ---------------- Daniel Lichtblau has pointed out that NullSpace does not generally give orthogonal vectors. Therefore the routines that depended upon that were in error. He says that it does give orthogonal vectors when the input vector contains approximate numbers. For graphical purposes this will be good enough for me. Therefore I modify Ted's routine to OrthogonalUnitVectors[vect__?(VectorQ[#, NumericQ] &)] /; (SameQ @@ Length /@ {vect}) && (Length[First[{vect}]] > 1) := #/Sqrt[#.#] & /@ NullSpace[{vect}// N] ---------------- Lets see what NullSpace does with approximate complex vectors. In[1]:= v1 = {1.0 I, 0.0, 0.5 I, 0.0, 1.0}; v2 = {0.0, 2.0, 1.0 I, 2.0, 0.5}; {v3,v4,v5} = NullSpace[{v1,v2}] Out[3]= {{-0.730153 + 0.*I, 0. - 0.138254*I, 0.250585 + 0.*I, 0. - 0.138254*I, 0. + 0.60486*I}, {0. + 0.*I, -0.515861 + 0.*I, 0. + 0.457321*I, 0.687357 + 0.*I, 0.22866 + 0.*I}, {0. + 0.*I, 0.510406 + 0.*I, 0. + 0.740442*I, -0.23274 + 0.*I, 0.370221 + 0.*I}} -------- In the next line we see NullSpace returned vectors that are orthogonal to the vectors we gave NullSpace. In[4]:= {v1.v3, v1.v4, v1.v5, v2.v3, v2.v4, v2.v5}//Chop Out[4]= {0, 0, 0, 0, 0, 0} ---------- However, the vectors returned aren't orthogonal to each other. In[5]:= {v3.v4, v3.v5, v4.v5}//Chop Out[5]= {0.229195*I, 0.371087*I, -0.677239} --------- I suppose an OrthogonalUnitVectors function that uses NullSpace should (1) Only accept real valued vectors. (2) Ensure NullSpace is given approximate vectors. ------ Ted Ersek ==== >OrthogonalUnitVectors[vect__?(VectorQ[#, NumericQ] &)] /; > (SameQ @@ Length /@ {vect}) && (Length[First[{vect}]] > 1) := > #/Sqrt[#.#] & /@ NullSpace[{vect}// N] ---------------- Lets see what NullSpace does with approximate complex vectors. In[1]:= > v1 = {1.0 I, 0.0, 0.5 I, 0.0, 1.0}; > v2 = {0.0, 2.0, 1.0 I, 2.0, 0.5}; > {v3,v4,v5} = NullSpace[{v1,v2}] Out[3]= > {{-0.730153 + 0.*I, 0. - 0.138254*I, 0.250585 + 0.*I, 0. - 0.138254*I, >0. >+ 0.60486*I}, > {0. + 0.*I, -0.515861 + 0.*I, 0. + 0.457321*I, 0.687357 + 0.*I, 0.22866 >+ 0.*I}, > {0. + 0.*I, 0.510406 + 0.*I, 0. + 0.740442*I, -0.23274 + 0.*I, 0.370221 >+ 0.*I}} -------- >In the next line we see NullSpace returned vectors that are orthogonal to >the vectors we gave NullSpace. In[4]:= > {v1.v3, v1.v4, v1.v5, v2.v3, v2.v4, v2.v5}//Chop Out[4]= > {0, 0, 0, 0, 0, 0} ---------- >However, the vectors returned aren't orthogonal to each other. In[5]:= > {v3.v4, v3.v5, v4.v5}//Chop Out[5]= > {0.229195*I, 0.371087*I, -0.677239} --------- >I suppose an OrthogonalUnitVectors function that uses NullSpace should > (1) Only accept real valued vectors. > (2) Ensure NullSpace is given approximate vectors. ------ > Ted Ersek I think you will find that the output vectors are orthogonal if you use the complex conjugate. for example v4.Conjugate[v5] is zero. Dennis Wangsness ==== >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. Needs[Graphics`FilledPlot`]; Needs[Graphics`Colors`]; FilledPlot[{4 - x^2, 3x}, {x, -4, 1}, Fills -> Red]; ==== >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, Your definition of b recalculates the integral for every call. To evaluate the integral once include Evaluate. Clear[f, b, FS]; L = 2; f[x_] := UnitStep[x - 1]; b[n_] := Evaluate[ (2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}]]; In the case of FS it is best to wait for an integer value of N prior to performing the Sum. FS[N_Integer, x_] := Sum[b[n]*Sin[n*Pi*x/L], {n, 1, N}]; Now Evaluate the argument of the Plot to cause the Sum to be done once. Timing[Plot[Evaluate[FS[120, x]], {x, 0, 2}]][[1]] 0.366667 Second While my computer may be faster than yours, this result for N=120 is 1000 times faster than your result for N=30. ==== > 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 Coefficient cannot figure out who is and is not a variable when a variable of 1 is specified. To work around this you might instead do delta1[poly_, x_] := Coefficient[poly,x]^2 - 4*Coefficient[poly,x,2]*Coefficient[poly,x,0] I prefer instead to use CoefficientList: delta2[poly_, x_] := (#[[2]]^2 - 4*#[[1]]*#[[3]])& [CoefficientList[poly,x]] In[25]:= poly = a*x^2+b*x+c; In[26]:= delta1[poly,x] === delta2[poly,x] Out[26]= True Daniel Lichtblau Wolfram Research ==== 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? see on-line help for Coefficient delta[eq_, x_] := Coefficient[eq, x, 1]^2 - 4 *Coefficient[eq, x, 2] *Coefficient[eq, x, 0]; eq = a*x^2 + b*x + c; delta[eq, x] b^2 - 4*a*c ==== Would someone with a very fast machine and lots of memory be willing to try this Solve for me? It is the inverse of the 20 node quadratic hexahedral mapping used in finite element analysis. None of my computers can handle this - they run out of memory (using (*Hex20 Node definition in global coordinates *) Clear[ x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, x5, y5, z5, x6, y6, z6, x7, y7, z7, x8, y8, z8, x9, y9, z9, x10, y10, z10, x11, y11, z11, x12, y12, z12, x13, y13, z13, x14, y14, z14, x15, y15, z15, x16, y16, z16, x17, y17, z17, x18, y18, z18, x19, y19, z19, x20, y20, z20]; (* local coordinates *) Clear[u, v, w]; (* Global co-ordinates *) Clear[x, y, z]; (* corner nodes *) N1= (1-u)*(1-v)*(1-w)*(-2-u-v-w)/8; N3= (1+u)*(1-v)*(1-w)*(-2+u-v-w)/8; N5= (1+u)*(1+v)*(1-w)*(-2+u+v-w)/8; N7= (1-u)*(1+v)*(1-w)*(-2-u+v-w)/8; N13=(1-u)*(1-v)*(1+w)*(-2-u-v+w)/8; N15=(1+u)*(1-v)*(1+w)*(-2+u-v+w)/8; N17=(1+u)*(1+v)*(1+w)*(-2+u+v+w)/8; N19=(1-u)*(1+v)*(1+w)*(-2-u+v+w)/8; (* to u nodes *) N2= (1-u^2)*(1-v)*(1-w)/4; N6= (1-u^2)*(1+v)*(1-w)/4; N14=(1-u^2)*(1-v)*(1+w)/4; N18=(1-u^2)*(1+v)*(1+w)/4; (* to v nodes *) N4= (1+u)*(1-v^2)*(1-w)/4; N8= (1-u)*(1-v^2)*(1-w)/4; N16=(1+u)*(1-v^2)*(1+w)/4; N20=(1-u)*(1-v^2)*(1+w)/4; (* to w nodes *) N9= (1-u)*(1-v)*(1-w^2)/4; N10=(1+u)*(1-v)*(1-w^2)/4; N11=(1+u)*(1+v)*(1-w^2)/4; N12=(1-u)*(1-v)*(1-w^2)/4; (* solve the inverse transform *) Solve[{ x1*N1+x2*N2+x3*N3+x4*N4+x5*N5+x6*N6+x7*N7+x8*N8+x9*N9+x10*N10+ x11*N11+x12*N12+x13*N13+x14*N14+x15*N15+x16*N16+x17*N17+x18*N18+x19*N19+x20* N20-x==0, y1*N1+y2*N2+y3*N3+y4*N4+y5*N5+y6*N6+y7*N7+y8*N8+y9*N9+y10*N10+ y11*N11+y12*N12+y13*N13+y14*N14+y15*N15+y16*N16+y17*N17+y18*N18+y19*N19+y20* N20-y==0, z1*N1+z2*N2+z3*N3+z4*N4+z5*N5+z6*N6+z7*N7+z8*N8+z9*N9+z10*N10+ z11*N11+z12*N12+z13*N13+z14*N14+z15*N15+z16*N16+z17*N17+z18*N18+z19*N19+z20* N20-z==0}, {u,v,w}] Christopher J. Purcell Defence R&D Canada ö Atlantic 9 Grove St., PO Box 1012 Dartmouth NS Canada B2Y 3Z7 ==== Needs[Graphics`Colors`]; ContourPlot[Sin[x y], {x, -5, 5}, {y, -5, 5}, ColorFunction -> (If[0.5 < # < 0.7, Red, White] &)]; >ContourPlot? >Jun Lin > In a message dated 9/6/02 3:53:58 AM, >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. >> Needs[Graphics`FilledPlot`]; >> Needs[Graphics`Colors`]; >> FilledPlot[{4 - x^2, 3x}, {x, -4, 1}, Fills -> Red]; ==== Jun Lin, Here is an example. Needs[Graphics`Colors`] Let's make a contour plot of this function. f[x_, y_] := Sin[x]Sin[2y] Let's specify the exact contours to use. I got rid of the 0. contour because it is difficult to obtain in this plot. contourvalues = Complement[Range[-1, 1, 0.2], {0.}] {-1, -0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8, 1.} Now we define a ColorFunction for the plot. I actually colored two different bands to show how you can make a general color function to give each band a desired color. cfun[z_] := Which[ -0.6 < z < -0.42, RoyalBlue, 0.4 < z < 0.6, Red, True, White] ContourPlot[f[x, y], {x, 0, Pi}, {y, 0, Pi}, PlotPoints -> 30, ColorFunctionScaling -> False, ColorFunction -> cfun, Contours -> contourvalues]; Using the option ColorFunctionScaling -> False says that the z value will be the actual value of f[x,y]. Otherwise, in general, it will be scaled between 0 and 1. It is easier to write a color function when z is the actual value of the function. Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.net Approved: Steven M. Christensen , Moderator Reply-To: ==== I checked again the six solutions I had previously timed, and they DO give orthogonal results. (None of them depend on NullSpace for that.) By the way, I reused my combinations function (from a recent problem on adding fractions to get 1) to check for orthogonality: ClearAll[orthogonalQ] orthogonalQ[v : {__?VectorQ}] := And @@ (Chop@( Dot @@ #) == 0 & /@ combinations[v, {2}]) << DiscreteMath`Combinatorica`; ClearAll[combinations]; r = Range[1, 9]; combinations::usage = combinations[list,n:{__Integer}] lists the combinations of list taken n at a time; combinations[r_List, n_Integer, {}] := If[n > Length@r, {}, DiscreteMath`Combinatorica`KSubsets[r, n]]; combinations[r_List, n_Integer, e_?VectorQ] := Join[e, #] & /@ DiscreteMath`Combinatorica`KSubsets[Complement[r, e], n]; combinations[r_List, n_Integer, e : {__?VectorQ}] := Flatten[ combinations[r, n, #] & /@ e, 1]; combinations[r_List, n : {__Integer}] := Which[Plus @@ n == Length@r, Join[#, Complement[r, #]] & /@ combinations[r, Drop[n, -1]], Plus @@ n > Length@r, {}, True, Fold[ combinations[r, #2, #1] &, {}, n]] Bobby -----Original Message----- (SameQ @@ Length /@ {vect}) && (Length[First[{vect}]] > 1) := #/Sqrt[#.#] & /@ NullSpace[{vect}// N] and the short version for 3D vectors OrthogonalUnitVectors[v : {_, _, _}] := #/Sqrt[#.#] & /@ NullSpace[{v//N}] For exact vectors I might use for 3D OrthogonalUnitVectors[v : {_, _, _}] := #/Sqrt[#.#] & /@ {temp = First[NullSpace[{v}]], v[Cross]temp} I'm still looking for something that is easy to remember. Park djmp@earthlink.net http://home.earthlink.net/~djmp/ OrthogonalUnitVectors[vect__?VectorQ]:= #/Sqrt[#.#]&/@NullSpace[{vect}] ------------- The version above will give a set of unit orthogonal vectors if given any number of vectors in any dimension. So besides giving it a 3D vector we can give it the following: OrthogonalUnitVectors[{2,1,0,-1,1}] or OrthogonalUnitVectors[{0,1,0,1/2,1},{1,0,-1,1/2}] ------------ But the short version above isn't very robust. (1) Clear[x,y,z];NullSpace[{{x,y,z}}] returns two vectors orthogonal to {x,y,z}, but the two vectors NullSpace returns aren't orthogonal to each other. So (OrthogonalUnitVectors) should only work with numeric vectors. (2) We should ensure all the vectors have the same dimension and length >1. I give a less concise version below that corrects these problems. ------------ OrthogonalUnitVectors[vect__?(VectorQ[#,NumericQ]&)]/; (SameQ@@Length/@{vect})&&(Length[First[{vect}]]>1):= #/Sqrt[#.#]&/@NullSpace[{vect}] -------------- Ted Ersek Get Mathematica tips, tricks from http://www.verbeia.com/mathematica/tips/Tricks.html ==== I have substituted the letters SS, CC, X, Y and Z1 and Z2 for some complicated expressions just to illustrate the form of the function. This function works and all the conditions are necessary but I am sure a more elegant programming solution perhaps using While could be found. Any f[{Sa_, Ca_, Aa_, Sb_, Cb_, Ab_, a_, b_}] := {If[Aa == Ab, a Sa + b Sb, If[Ca == 0 || Cb == 0, a Sa + b Sb, SS]], If[Aa == Ab, a Ca + b Cb, If[Ca == 0 || Cb == 0, a Ca + b Cb, CC]], If[Ca == 0 && Cb == 0, Nil, If[Z1 == Z2, If[Ca < 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} ==== Check out the Help for Coefficient. Coefficient[expr, form, 0] picks out terms that are not proportional to form. Delta[eq_, x_] := Coefficient[eq, x]^2 - 4 Coefficient[eq, x^2] Coefficient[eq, x, 0] Delta[a x^2 + b x + c, x] b^2 - 4*a*c Park djmp@earthlink.net http://home.earthlink.net/~djmp/ ==== Dear all, 1.) How to find the General Solution for below's partial differential equation? (y + u) du/dx + y (du/dy) = x - y ** I use d to represent the partial differential symbol. Can it be solved by function NDSolve in mathematica 4.1? How? 2.) I manage to get the roots of complex equation z^5 = i from Solve[z^5 == i, z]. It gave me straight the 5 roots in the output. Is there any way to view the steps in mathematica? ============================================================================ ====== of the individual or entity to which they are addressed. Any disclosure, copying, distribution and diversion contrary to the applicable export control laws and regulations including US Export Administration Regulations is strictly prohibited. and do not disclose it to others. Please notify the postmaster@hitachi.com.my of the delivery error by replying to this message and then delete it from your ============================================================================ ====== ==== I have been trying to code Sethian's Fast Marching Method in 2D but Mathematica has been very slow (taking something like 1-2 hours for something that should take much less than a second in C++). I am sure part of the problem my time. I looked at the list archives and there was mention of an profiling package for Mathematica but a)I can't find it & b)It may not work with Mathematica 4.*. My questions: 1. Any general suggestions on how to figure out which functions are taking most of the time? I guess I could manually have each function I am interested in monitoring keep a variable that counts the amount of CPU time that has been spent on it by doing something like: function[args_]:=Module[{},functionTimer+=Timing[ .... My actual function ..... ][[1]]] but it would very cumbersome to do this to all of the functions in my program and I am not sure I will get accurate results anyway. 2. Compiling functions is not always that easy. I did read the on-line docs and the archives and it does take some work to make a function compile usefully. Is there an FAQ or a tutorial somewhere? 3. Am I the only one who finds the lack of a profiler really really annoying ? Mathematica is powerful and it is usually easy to ask it to do what u want. The challenge a lot of times is doing so without taking too long. Husain PS: One more quick one: Why does the front end act funny when I have the Reply-To: jmt@dxdydz.net ==== mathematica -primaryModifierMask etc see man mathematica for other (very useful) options PS: One more quick one: Why does the front end act funny when I have the Your X server is set up so that the NumLock key is mapped to Mod2. I > learned by accident that Mod2 is actually quite useful. Mod2-click on a > cell selects all cells of that type in the current notebook. This is an > easy way to delete all the graphics cells and output cells in a notebook to > reduce file size. An annoying aspect is that Mod2-click means press > NumLock, click, then press NumLock again to turn it off. You should be able > to map a different key to Mod2 using xmodmap. Remapping modifier keys in X > is awfully annoying, though. I recommend xkeycaps > (http://www.jwz.org/xkeycaps/). > Reply-To: kuska@informatik.uni-leipzig.de ==== a Mathematica profiler is described in The Mathematica Journal Volume 5, Issue 3, Summer 1995 The Mathematica Toolbox: A Mathematica Profiler by Todd Gayley The electronic material for this issue isn not on MathSource but it may be that Todd has the code some where and can make is acessible. > I have been trying to code Sethian's Fast Marching Method in 2D but Mathematica > has been very slow (taking something like 1-2 hours for something that > should take much less than a second in C++). I am sure part of the problem > my time. I looked at the list archives and there was mention of an profiling > package for Mathematica but a)I can't find it & b)It may not work with Mathematica 4.*. > My questions: > 1. Any general suggestions on how to figure out which functions are taking > most of the time? I guess I could manually have each function I am > interested in monitoring keep a variable that counts the amount of CPU > time that has been spent on it by doing something like: > function[args_]:=Module[{},functionTimer+=Timing[ .... My actual function > ..... ][[1]]] > but it would very cumbersome to do this to all of the functions in my > program and I am not sure I will get accurate results anyway. > 2. Compiling functions is not always that easy. I did read the on-line > docs and the archives and it does take some work to make a function > compile usefully. Is there an FAQ or a tutorial somewhere? > 3. Am I the only one who finds the lack of a profiler really really > annoying ? Mathematica is powerful and it is usually easy to ask it to do what u > want. The challenge a lot of times is doing so without taking too long. Husain > PS: One more quick one: Why does the front end act funny when I have the ==== > PS: One more quick one: Why does the front end act funny when I have the Your X server is set up so that the NumLock key is mapped to Mod2. I learned by accident that Mod2 is actually quite useful. Mod2-click on a cell selects all cells of that type in the current notebook. This is an easy way to delete all the graphics cells and output cells in a notebook to reduce file size. An annoying aspect is that Mod2-click means press NumLock, click, then press NumLock again to turn it off. You should be able to map a different key to Mod2 using xmodmap. Remapping modifier keys in X is awfully annoying, though. I recommend xkeycaps (http://www.jwz.org/xkeycaps/). ==== Caution: the definition b[n_] := Evaluate[(2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}]]; doesn't immediately compute the integral unless f is already defined at this point, and in that case you may as well write b[n_] = (2/L)*Integrate[f[x]*Sin[n*Pi*x/L], {x, 0, L}]; instead. If b will be computed more than once for the same n, its even better to do it THIS way (if f and L do not change): f[x_] = Cos[x] (* for example *) Simplify[Integrate[f[x]*Sin[n*Pi*(x/L)], {x, 0, L}]] (L*((-n)*Pi + n*Pi*Cos[L]*Cos[n*Pi] + L*Sin[L]*Sin[n*Pi]))/ ((L - n*Pi)*(L + n*Pi)) b[n_] := b[n] = (L*((-n)*Pi + n*Pi*Cos[L]*Cos[n*Pi] + L*Sin[L]*Sin[n*Pi]))/ ((L - n*Pi)*(L + n*Pi)) Bobby -----Original Message----- Stephan Steinhaus (steinhaus-net.de). So when Mathematica takes a very long time, you should investigate. In this case inserting Evaluate[] in two places In[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 my feeble old 500MHz PowerBook). Why was it so slow before? When I switch from an ordinary numerical language to Mathematica, I enter into an implicit bargain with Mathematica: 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 being asked) 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 about it, except for the bad answers that result now and then. But I need to take care that Mathematica does not burden itself unnecessarily. That's my side of the bargain. Number (2) is the issue here. Your definition of b[n] is written so that Mathematica analytically evaluates b separately for each n. But you know in this 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 some cases, the most common of which involve branching within the definition of function 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 to steer it into the shortcut when it's OK. Hope this helps, Burton -- ==== With the second method below (mine), the times add up to 40% less than with the first method (-Peer's), for what SEEMS to be exactly the same work. Go figure! Can anybody explain that? (* -Peer *) ClearAll[f, b, FS] L = 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}]] {0.547 Second, .89ª°Graphics.89ª°} (* Treat #1 *) Timing[(2/L)*Integrate[f[x]*Sin[n*Pi*(x/L)], {x, 0, L}]] {0.016000000000000014*Second, (2*(Cos[(n*Pi)/2] - Cos[n*Pi]))/(n*Pi)} ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = (2*(Cos[(n*Pi)/2] - Cos[n*Pi]))/(n*Pi); FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] {0.297 Second, .89ª°Graphics.89ª°} Even stranger, it SLOWS the Plot if we precompute b before Timing starts: (* Treat #2 *) ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = (2*(Cos[(n*Pi)/2] - (-1)^n))/(n*Pi); b /@ Range[30]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] {0.328 Second, .89ª°Graphics.89ª°} But here's a winner: (* Treat #3 *) ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = N[(2*(Cos[(n*Pi)/2] - (-1)^n))/ (n*Pi)]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] {0.204 Second, .89ª°Graphics.89ª°} Apparently, computing b within the Plot causes machine-precision arithmetic to be used, and that saves time. Precomputing b and then converting exact expressions to approximate ones within the Plot seems to take longer. For that to make sense, I think it must be that n is approximate (not Integer) when it is passed to b within Plot. -----Original Message----- and you need only a 1-3 seconds (depending on your machine) 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 ==== Even better -- MUCH better -- add 's Evaluate to the other tricks: ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = N[(2*(Cos[(n*Pi)/2] - (-1)^n))/(n*Pi)]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[Evaluate[FS[30, x]], {x, 0, 2}]] {0.015 Second, .89ª°Graphics.89ª°} ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = N[(2*(Cos[(n*Pi)/2] - (-1)^n))/(n*Pi)]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[Evaluate[FS[120, x]], {x, 0, 2}]] {0.063 Second, .89ª°Graphics.89ª°} -----Original Message----- Timing[Plot[FS[30, x], {x, 0, 2}]] {0.547 Second, .89ª°Graphics.89ª°} (* Treat #1 *) Timing[(2/L)*Integrate[f[x]*Sin[n*Pi*(x/L)], {x, 0, L}]] {0.016000000000000014*Second, (2*(Cos[(n*Pi)/2] - Cos[n*Pi]))/(n*Pi)} ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = (2*(Cos[(n*Pi)/2] - Cos[n*Pi]))/(n*Pi); FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] {0.297 Second, .89ª°Graphics.89ª°} Even stranger, it SLOWS the Plot if we precompute b before Timing starts: (* Treat #2 *) ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = (2*(Cos[(n*Pi)/2] - (-1)^n))/(n*Pi); b /@ Range[30]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] {0.328 Second, .89ª°Graphics.89ª°} But here's a winner: (* Treat #3 *) ClearAll[f, b, FS] L = 2; f[x_] := UnitStep[x - 1]; b[n_] := b[n] = N[(2*(Cos[(n*Pi)/2] - (-1)^n))/ (n*Pi)]; FS[N_, x_] := Sum[b[n]*Sin[n*Pi*(x/L)], {n, 1, N}]; Timing[Plot[FS[30, x], {x, 0, 2}]] {0.204 Second, .89ª°Graphics.89ª°} Apparently, computing b within the Plot causes machine-precision arithmetic to be used, and that saves time. Precomputing b and then converting exact expressions to approximate ones within the Plot seems to take longer. For that to make sense, I think it must be that n is approximate (not Integer) when it is passed to b within Plot. -----Original Message----- and you need only a 1-3 seconds (depending on your machine) 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 ==== How can I plot with Mathematica two function in the same graphic? I explain better: compare a and b). Reply-To: kuska@informatik.uni-leipzig.de ==== Plot[{x,x^2},{x,-1,1}] may do it. How can I plot with Mathematica two function in the same graphic? > I explain better: > compare a and b). > ==== Other than Plot[{x,x^2},...] you can use Show[] as well. plt1=Plot[x,{x,0,1}]; plt2=Plot[x^2,{x,0,1}]; Show[plt1,plt2] Lawrence > How can I plot with Mathematica two function in the same graphic? > I explain better: > compare a and b). > Mario, Here is a fancy version of your plot. I used Text statements within an > Epilog option to label the two curves. The regular plot statement allows you > to plot a series of functions inclosed in a list. Needs[Graphics`Colors`] > Plot[{x, x^2}, {x, 0, 1}, > PlotStyle -> {Black, Blue}, > Frame -> True, > FrameLabel -> {x, y}, > PlotLabel -> Comparison of Two Functions, > Epilog -> {Text[x, {0.5, 0.55}], Blue, Text[x^2, {0.7, 0.4}]}, > Background -> Linen, > ImageSize -> 500]; Park > djmp@earthlink.net > http://home.earthlink.net/~djmp/ How can I plot with Mathematica two function in the same graphic? > I explain better: > compare a and b). > ==== (1) Is there a way in Mathematica 4.2 to put two separate gifs into a single, cell side by side (with some intervening space), without having to combine them in some graphics program first? In particular, I'd like to do that within a text cell. Even in a new Input cell, if I first create a GridBox (via Inut>Create Table/Matrix/Palette) and then try to insert the first gif (via Edit>Insert Object>Create from File ....), Mathematica promptly crashes. (Mathematica 4.2 under Windows 2000.) I just don't see how to get anything other than a single gif into a cell. (2) Is there a way to cause a gif imported into a Mathematica 4.2 notebook to become a hyperlink -- so that when the user clicks on the gif the hyperlink's target is summoned? (My aim in all this is to use Mathematica to create web pages with high mathematical content -- saving the notebook as HTML+MathML -- without having to do any extensive editing of the resulting .xml and related files. That way as the source Mathematica notebook changes, I would need only to re-export without further tinkering with the .xml file, etc.) -- Murray Eisenberg Internet: murray@math.umass.edu Mathematics & Statistics Dept. Voice: 413-545-2859 (W) University of Massachusetts 413-549-1020 (H) Reply-To: kuska@informatik.uni-leipzig.de ==== it works with any graphics NotebookWrite[SelectedNotebook[], GridBox[{Cell[GraphicsData[PostScript, DisplayString[#, MPS]], Graphics] & /@ {leftGraphics, rightGraphics}}]] (1) Is there a way in Mathematica 4.2 to put two separate gifs into a > single, cell side by side (with some intervening space), without > having to combine them in some graphics program first? In particular, I'd like to do that within a text cell. Even in a new Input cell, if I first create a GridBox (via Inut>Create > Table/Matrix/Palette) and then try to insert the first gif (via > Edit>Insert Object>Create from File ....), Mathematica promptly > crashes. (Mathematica 4.2 under Windows 2000.) I just don't see how to get anything other than a single gif into a > cell. (2) Is there a way to cause a gif imported into a Mathematica 4.2 > notebook to become a hyperlink -- so that when the user clicks on the > gif the hyperlink's target is summoned? (My aim in all this is to use Mathematica to create web pages with > high mathematical content -- saving the notebook as HTML+MathML -- > without having to do any extensive editing of the resulting .xml and > related files. That way as the source Mathematica notebook changes, I > would need only to re-export without further tinkering with the .xml > file, etc.) -- > Murray Eisenberg Internet: murray@math.umass.edu > Mathematics & Statistics Dept. Voice: 413-545-2859 (W) > University of Massachusetts 413-549-1020 (H) ==== Your method 1), below, does not produce correct results if one then uses menu>Save As Special>HTML+MathML. I view the resulting .xml file in a MathML-enabled browser (e.g., Mozilla with appropriate TrueType fonts installed -- the four BaKoMa cm fonts, the old Mathematica Math1, etc., fonts, and the two MT fonts. Then the everything in the originally multi-line inline cell appears on one line separate by a ?, as do the alignment markers. In the page source, these unrendered symbols have codes #8289 and #63328, respectively. Is there an issue of the encoding used in the browser here? I did try Unicode-7, Unicode-8, and several of the Western encodings available under Mozilla's Default Character Encoding. Or is it something else? Note that I do also have the new Mathematica TrueType fonts (Mathematica1, etc.). P.S. I don't think Mathematica documents this, but if one wants a file exported from Mathematica as HTML+MathML to be rendered correctly by a browser, it seems to require an extension of .xml rather than .html. > >>When writing a series of equations in mathematica (in a text cell) >>is there any way to align the equations at the = similar to >>what can be done in the equation editor made by math type (used >>in word etc.)? > > > Mike, > Two ways > 1) Start an inline cell in your text cell (menu>Edit>Expression Input>Start > Inline Cell) > Type in your equations with the alignment marker (Esc am Esc) in after each > equal sign. > Close the inline cell (menu>Edit>Expression Input>End Inline Cell). > Select the text cell (or the inline cell, though this is a delicate > operation, do it by repeatedly double clicking in it}. > Use menu> Format>Text Alignment>On Alignment Marker. > -- you will find that all single letters are now italic - you can select > individual letters and change this - alternatively you can select the inline > cell and use the option inspector (menu>Format>Option Inspector) to set > SingleLetterItalics->False) > > -- you can put existing text in an inline cell by selcting it and using > menu>Edit>Expression Input>Start Inline Cell, but you may have to adjust > the line breaks after this.... -- Murray Eisenberg Internet: murray@math.umass.edu Mathematics & Statistics Dept. Voice: 413-545-2859 (W) University of Massachusetts 413-549-1020 (H) ==== Group, Can anyone help me with the following? I guess this is something trivial for most of you, but I have being struggling with it for weeks. The problem is to construct a taxonomic hierarchy whose structure is hereby described by example: << DiscreteMath`Combinatorica` n[1]:= m = 10; (*with the following I get the 1th level (1--element) partition*) in[2]:= ks1 = Partition[Range[m], 1] (*Now the following gives all the 2th-level (?) different partitions that can be set up by joining two elements in the 1-element partition and leaving the rest as they are?*) (*First, these are all possible 'joinings' to generate the two elements 'pieces' from 1th level (1--element) partition*) in[3]:= ks2 = KSubsets[Range[m], 2] (*And this generates all the 2th-level (?) different partitions that can be set up by joining two elements in the 1-element partition and leaving the rest as they are.*) in[4]:= p2 = MapThread[Complement[Append[#1, #2], (Sequence @@ {#} & /@ Partition[#2, 1])] &, {Array[ks1&, Length[ks2]], ks2}]; Length[p2] in[5]:= p2[[1]] out[5]:= {{3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2}} in[6]:= p2[[Length[p2]]] out[6]:= {{1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9, 10}} (*compute in[4] t and see the full output*) (*Now I need to solve the problem of given the ith-level (?) different partitions to set up all the i+1th different partitions by joining two elements in the ith level partition and leaving the rest as they are? *) (*for example given the following 2th-level partition*) {3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2} (*we should get*) {{{4}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2, 3}}, {{3}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2, 4}}, {{3}, {4}, {6}, {7}, {8}, {9}, {10}, {1, 2, 5}}, {{3}, {4}, {5}, {7}, {8}, {9}, {10}, {1, 2, 6}}, {{3}, {4}, {5}, {6}, {8}, {9}, {10}, {1, 2, 7}}, {{3}, {4}, {5}, {6}, {7}, {9}, {10}, {1, 2, 8}}, {{3}, {4}, {5}, {6}, {7}, {8}, {10}, {1, 2, 9}}, {{3}, {4}, {5}, {6}, {7}, {8}, {9}, {1, 2, 10}}} (*The process would end up to *) {1,{2,3,4,5,6,7,8,9,10}} {2,{1,3,4,5,6,7,8,9,10}} {3,{1,2,4,5,6,7,8,9,10}} {4,{1,2,3,5,6,7,8,9,10}} {5,{1,2,3,4,6,7,8,9,10}} {6,{1,2,3,4,5,7,8,9,10}} {7,{1,2,3,4,5,6,8,9,10}} {8,{1,2,3,4,5,6,7,9,10}} {9,{1,2,3,4,5,6,7,8,10}} {10,{1,2,3,4,5,6,7,8,9}} (* and*) {1,2,3,4,5,6,7,8,9,10} Emilio Martin-Serrano ==== > I too have been trying to use Mathematica (v4.2 most recently) to type > maths papers and the like but I'm not ready to ditch LaTeX yet. There > are just too many cases where I cannot figure out how to achieve what I > want in Mathematica, things like: > - left brackets spanning multiple lines for defining hybrid functions; You can accomplish this by doing the following: 1) Put your function braches in the rows of a grid box structure. 2) Add the following options to your cell: ShowAutoStyles -> False SpanMaxSize -> Infinity The following cell snippet demonstrates how this influences the result. To view it, paste the Cell[] expression into a notebook and then click on Yes when you are prompted on whether the front end should interpret the result. Cell[BoxData[ FormBox[ RowBox[{ RowBox[{f, (, x, )}], =, RowBox[{{, GridBox[{ {x, RowBox[{x, , <, 0}]}, { SuperscriptBox[x, 2], RowBox[{0, [LessEqual], x, <, 1}]}, { RowBox[{sin, (, x, )}], RowBox[{1, [LessEqual], x, <, 2}]}, { RowBox[{[CapitalGamma], (, x, )}], RowBox[{x, [GreaterEqual], 2}]} }]}]}], TraditionalForm]], DisplayFormula, ShowAutoStyles->False, SpanMaxSize->Infinity] > - vertical alignment of equals signs in multi-line equations or > derivations; Put your equations in a GridBox and set the ColumnAlignments option to a string containing the equal sign. Cell[BoxData[ FormBox[GridBox[{ { RowBox[{ RowBox[{ RowBox[{3, x}], , +, , RowBox[{4, , y}]}], , =, , 9}]}, { RowBox[{ RowBox[{ RowBox[{2, x}], , -, , RowBox[{7, , y}]}], =, RowBox[{32, , -, , RowBox[{sin, (, x, )}]}]}]} }], TraditionalForm]], DisplayFormula, GridBoxOptions->{ColumnAlignments->{=}}] > - setting typefaces in tables of material. I think the Author Tools material that comes with Mathematica 4.2 might be able to help you do this. -- User Interface Programmer paulh@wolfram.com Wolfram Research, Inc. ==== You have three nonlinear (fourth-order) equations and 23 unknowns. A faster computer won't help any. Bobby -----Original Message----- (*Hex20 Node definition in global coordinates *) Clear[ x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, x5, y5, z5, x6, y6, z6, x7, y7, z7, x8, y8, z8, x9, y9, z9, x10, y10, z10, x11, y11, z11, x12, y12, z12, x13, y13, z13, x14, y14, z14, x15, y15, z15, x16, y16, z16, x17, y17, z17, x18, y18, z18, x19, y19, z19, x20, y20, z20]; (* local coordinates *) Clear[u, v, w]; (* Global co-ordinates *) Clear[x, y, z]; (* corner nodes *) N1= (1-u)*(1-v)*(1-w)*(-2-u-v-w)/8; N3= (1+u)*(1-v)*(1-w)*(-2+u-v-w)/8; N5= (1+u)*(1+v)*(1-w)*(-2+u+v-w)/8; N7= (1-u)*(1+v)*(1-w)*(-2-u+v-w)/8; N13=(1-u)*(1-v)*(1+w)*(-2-u-v+w)/8; N15=(1+u)*(1-v)*(1+w)*(-2+u-v+w)/8; N17=(1+u)*(1+v)*(1+w)*(-2+u+v+w)/8; N19=(1-u)*(1+v)*(1+w)*(-2-u+v+w)/8; (* to u nodes *) N2= (1-u^2)*(1-v)*(1-w)/4; N6= (1-u^2)*(1+v)*(1-w)/4; N14=(1-u^2)*(1-v)*(1+w)/4; N18=(1-u^2)*(1+v)*(1+w)/4; (* to v nodes *) N4= (1+u)*(1-v^2)*(1-w)/4; N8= (1-u)*(1-v^2)*(1-w)/4; N16=(1+u)*(1-v^2)*(1+w)/4; N20=(1-u)*(1-v^2)*(1+w)/4; (* to w nodes *) N9= (1-u)*(1-v)*(1-w^2)/4; N10=(1+u)*(1-v)*(1-w^2)/4; N11=(1+u)*(1+v)*(1-w^2)/4; N12=(1-u)*(1-v)*(1-w^2)/4; (* solve the inverse transform *) Solve[{ x1*N1+x2*N2+x3*N3+x4*N4+x5*N5+x6*N6+x7*N7+x8*N8+x9*N9+x10*N10+ x11*N11+x12*N12+x13*N13+x14*N14+x15*N15+x16*N16+x17*N17+x18*N18+x19*N19+ x20*N20-x==0, y1*N1+y2*N2+y3*N3+y4*N4+y5*N5+y6*N6+y7*N7+y8*N8+y9*N9+y10*N10+ y11*N11+y12*N12+y13*N13+y14*N14+y15*N15+y16*N16+y17*N17+y18*N18+y19*N19+ y20*N20-y==0, z1*N1+z2*N2+z3*N3+z4*N4+z5*N5+z6*N6+z7*N7+z8*N8+z9*N9+z10*N10+ z11*N11+z12*N12+z13*N13+z14*N14+z15*N15+z16*N16+z17*N17+z18*N18+z19*N19+ z20*N20-z==0}, {u,v,w}] Christopher J. Purcell Defence R&D Canada - Atlantic 9 Grove St., PO Box 1012 Dartmouth NS Canada B2Y 3Z7 ==== If I haven't missed a step, the following should work nicely -- best done before defining any of the symbols that appear: Attributes[dummyIf] = HoldAll; Attributes[dummyWhich] = HoldAll; expr = {If[Aa == Ab, a Sa + b Sb, If[Ca == 0 || Cb == 0, a Sa + b Sb, SS]], If[Aa == Ab, a Ca + b Cb, If[Ca == 0 || Cb == 0, a Ca + b Cb, CC]], If[Ca == 0 && Cb == 0, Nil, If[Z1 == Z2, If[Ca < 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[ Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} /. {If -> dummyIf}; rule1 = dummyIf[a_, b_, dummyIf[c_, d_, e_]] -> dummyWhich[a, b, c, d, True, e]; rule2 = dummyIf[a_, dummyIf[ b_, c_, d_], e_] -> dummyWhich[a && b, c, a, d, True, e]; rule3 = dummyIf[a_, dummyIf[b_, c_, d_]] -> dummyWhich[a && b, c, a, d]; rule4 = dummyWhich[a__, b_, dummyIf[c_, d_, e_], f___] -> dummyWhich[a, b && c, d, b, e, f]; rule5 = dummyWhich[a__, b_, dummyWhich[c_, d_, e___], f___] -> dummyWhich[a, b && c, d, b, dummyWhich[e], f]; rule6 = dummyWhich[] :> Null; rule7 = HoldPattern[True && a_] :> a; expr //. {rule1, rule2, rule3, rule4, rule5, rule6, rule7} /. {dummyWhich -> Which The result of that last line is: {Which[Aa == Ab, a Sa + b Sb, Ca == 0 || Cb == 0, a Sa + b Sb, True, SS], Which[Aa == Ab, a Ca + b Cb, Ca == 0 || Cb == 0, a Ca + b Cb, True, CC], Which[Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, X > 180, X - 180, X < 0, X + 180, True, X]} I've used Null where that would be the result of your original logic -- maybe you want it to be Nil. Or maybe you want to use Null instead of Nil. Your choice. Wherever you see Null, there's possibly a case you haven't covered. -----Original Message----- If[Aa == Ab, a Ca + b Cb, If[Ca == 0 || Cb == 0, a Ca + b Cb, CC]], If[Ca == 0 && Cb == 0, Nil, If[Z1 == Z2, If[Ca < 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} ==== Dear Group, I have a program to take the local polynomial non parametric regression of two variables. It uses Compile, unfortunately, and regularly, but inconsistently, causes Mathematica to crash (in Win2K, with I forget what error, and in Win98/Mathematica4.0 with an invalid memory access from MathDLL.dll). I am running Mathematica 4.1, and have this problem consistently on 4 different machines. Here is the code, if it doesn't crash the first time, it will the second or third: (* w = Compile[{{xj, _Real, 0}, {XX, _Real, 1}, {YY, _Real, 1}, {h, _Real, 0}, {nn, _Integer, 0}, {ord, _Integer, 0}}, First[Inverse[Sum[Outer[Times, Table[If[Positive[q], (XX[[i]] - xj)^q, 1], {q, 0, ord}], Table[If[Positive[q], (XX[[i]] - xj)^q, 1.], {q, 0, ord}]*E^(-0.5*((XX[[i]] - xj)/h)^2)], {i, nn}]] . Sum[(Table[If[Positive[q], (XX[[i]] - xj)^q, 1.], {q, 0, ord}]*YY[[i]])* E^(-0.5*((XX[[i]] - xj)/h)^2), {i, nn}]]]; *) (* !(tt = MemoryInUse[]; ListPlot[Table[{((i + 1))^2, (Timing[ nn = ((i + 1))^2; [Epsilon] = Table[Random[]*3, {i, nn}]; testX = Table[i* .3 - Random[]*2, {i, nn}]; testY = Table[Sin[i/3] - 2, {i, nn}] + [Epsilon]; Map[w[#, testX, testY, .1 + Random[], nn, 0] &, testX];])[([1])]/Second}, {i, 15}], PlotJoined -> True, PlotLabel -> ]; MemoryInUse[] - tt) *) I have followed Ted Ersek's tips and tricks (http://www.verbeia.com/mathematica/tips/tip_index.html) about Compile, and have tried changing all variable names (works sometimes), removing any hidden spaces, restructuring the formulas, changing all 0's to 0. 's and all 1's to 1. 's, etcetera etcetera, but I still can't comprehend what the problem might be. Doing w[[-2]] shows a list of op-code numbers, and one function name, Inverse[#1]&, so it seems to me that there are no problems with the use of Compile here. Would anyone have any thoughts??? Bernard Gress burnthebiscuit@netscape.net ==== -Peer Kuska tells me that he gets no crash with Mathematica 4.2, however after working with this program for more than a month, trying every possible permutation and manipulation, and having it still crash, I am still worried that upgrading to 4.2 won't solve my problem. Can anyone else with 4.2 try running this program for me a good number of times, say 3 or 4, and see if they get similar good results? Bernard Gress Dear Group, I have a program to take the local polynomial non parametric regression of two variables. It uses Compile, unfortunately, and regularly, but inconsistently, causes Mathematica to crash (in Win2K, with I forget what error, and in Win98/Mathematica4.0 with an invalid memory access from MathDLL.dll). I am running Mathematica 4.1, and have this problem consistently on 4 different machines. Here is the code, if it doesn't crash the first time, it will the second or third: (* w = Compile[{{xj, _Real, 0}, {XX, _Real, 1}, {YY, _Real, 1}, {h, _Real, 0}, {nn, _Integer, 0}, {ord, _Integer, 0}}, First[Inverse[Sum[Outer[Times, Table[If[Positive[q], (XX[[i]] - xj)^q, 1], {q, 0, ord}], Table[If[Positive[q], (XX[[i]] - xj)^q, 1.], {q, 0, ord}]*E^(-0.5*((XX[[i]] - xj)/h)^2)], {i, nn}]] . Sum[(Table[If[Positive[q], (XX[[i]] - xj)^q, 1.], {q, 0, ord}]*YY[[i]])* E^(-0.5*((XX[[i]] - xj)/h)^2), {i, nn}]]]; *) (* !(tt = MemoryInUse[]; ListPlot[Table[{((i + 1))^2, (Timing[ nn = ((i + 1))^2; [Epsilon] = Table[Random[]*3, {i, nn}]; testX = Table[i* .3 - Random[]*2, {i, nn}]; testY = Table[Sin[i/3] - 2, {i, nn}] + [Epsilon]; Map[w[#, testX, testY, .1 + Random[], nn, 0] &, testX];])[([1])]/Second}, {i, 15}], PlotJoined -> True, PlotLabel -> ]; MemoryInUse[] - tt) *) Reply-To: kuska@informatik.uni-leipzig.de ==== my Mathematica 4.2 does not crash and you may upgrade your version. Dear Group, I have a program to take the local polynomial non parametric regression > of two variables. It uses Compile, unfortunately, and regularly, but > inconsistently, causes Mathematica to crash (in Win2K, with I forget > what error, and in Win98/Mathematica4.0 with an invalid memory access from > MathDLL.dll). I am running Mathematica 4.1, and have this problem consistently > on 4 different machines. Here is the code, if it doesn't crash the first time, it will the second > or third: (* > w = Compile[{{xj, _Real, 0}, {XX, _Real, 1}, {YY, _Real, 1}, {h, _Real, > 0}, > {nn, _Integer, 0}, {ord, _Integer, 0}}, > First[Inverse[Sum[Outer[Times, Table[If[Positive[q], (XX[[i]] - xj)^q, > 1], {q, 0, ord}], > Table[If[Positive[q], (XX[[i]] - xj)^q, 1.], {q, 0, > ord}]*E^(-0.5*((XX[[i]] - xj)/h)^2)], {i, nn}]] . > Sum[(Table[If[Positive[q], (XX[[i]] - xj)^q, 1.], {q, 0, ord}]*YY[[i]])* > E^(-0.5*((XX[[i]] - xj)/h)^2), {i, nn}]]]; > *) (* > !(tt = MemoryInUse[]; > ListPlot[Table[{((i + 1))^2, (Timing[ > nn = ((i + 1))^2; [Epsilon] = Table[Random[]*3, {i, > nn}]; > testX = Table[i* .3 - Random[]*2, {i, nn}]; > testY = Table[Sin[i/3] - 2, {i, nn}] + [Epsilon]; > Map[w[#, testX, testY, .1 + Random[], nn, 0] &, > testX];])[([1])]/Second}, {i, 15}], PlotJoined - True, > PlotLabel -> ]; > MemoryInUse[] - tt) > *) I have followed Ted Ersek's tips and tricks > (http://www.verbeia.com/mathematica/tips/tip_index.html) about Compile, > and have tried changing all variable names (works sometimes), removing > any hidden spaces, restructuring the formulas, changing all 0's to 0. 's > and all 1's to 1. 's, etcetera etcetera, but I still can't comprehend > what the problem might be. Doing w[[-2]] shows a list of op-code > numbers, and one function name, Inverse[#1]&, so it seems to me that > there are no problems with the use of Compile here. Would anyone have > any thoughts??? > Bernard Gress > burnthebiscuit@netscape.net ==== Mr Kuska, 4.2, I thought I had solved this problem so many times and still it would suddenly start crashing again our of the blue, after 20 or 30 executions. Bernard my Mathematica 4.2 does not crash and you may upgrade > your version. > ==== Expanding on my earlier simplification of the Moran Research code, I studied the problem a little more and found another useful Rule and a Rule that should help but doesn't. Maybe somebody can explain what I'm doing wrong on that one. My earlier result (after applying rule1, ... rule7) had the following as third element of a List in the definition of f: Which[ Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, X > 180, X - 180, X < 0, X + 180, True, X] Several things occurred to me in terms of simplifying this. The first was that the last six arguments of Which could be replaced with two, as follows: Which[ Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True,Mod[X,180]] This probably has no advantage other than clarity, but that's worth the others: rule8 = dummyIf[a_ > b_, a_ - b_, dummyIf[a_ < 0, a_ + b_, a_]] :> Mod[a, b]; but it had no effect on the expression. Apparently there was no match, and I can't see why. The second thing I noticed was that the condition Ca == 0 && Cb == 0 occurs twice in the Which statement. The following rule fixes that kind of situation: rule9 = dummyWhich[a___, b_, c_, d__, b_, e_, f__] /; EvenQ[Length[List[a]]] && EvenQ[Length[List[d]]] && EvenQ[Length[List[f]]] :> dummyWhich[a, b, c, d, f]; expr //. {rule1, rule2, rule3, rule4, rule5, rule6, rule7, rule9} The third thing I noticed was the treatment of Aa and Ab. The expressions If[Aa > 90, Aa - 90, Aa + 90] If[Ab > 90, Ab - 90, Ab + 90] in the original expression may be equivalent (depending on what's known about Aa and Ab a priori) to: Mod[Aa+90,180] Mod[Ab+90,180] If that is a valid assumption in the situation (no way for me to know), the expression is now Which[ Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, Z1 == Z2 && Ca > 0, Mod[Aa + 90, 180], Z1 == Z2 && Cb > 0, Mod[Ab + 90, 180], Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True, Mod[X, 180]] The next thing I notice is that the sixth condition, Z1==Z2, which results in Null if found True, isn't needed. If Z1==Z2 then one of the first five conditions is also true, so execution wouldn't get that far. Hence we can delete that condition-response pair. Finally, I'm stopping with this: Which[ Ca == 0 == Cb, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, Z1 == Z2 && Ca > 0, Mod[Aa + 90, 180], Z1 == Z2 && Cb > 0, Mod[Ab + 90, 180], Aa == Ab, Aa, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True, Mod[X, 180]] I have no idea what you're actually doing with this code, but it looks weird even AFTER I've simplified it. -----Original Message----- 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[ Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} /. {If -> dummyIf}; rule1 = dummyIf[a_, b_, dummyIf[c_, d_, e_]] -> dummyWhich[a, b, c, d, True, e]; rule2 = dummyIf[a_, dummyIf[ b_, c_, d_], e_] -> dummyWhich[a && b, c, a, d, True, e]; rule3 = dummyIf[a_, dummyIf[b_, c_, d_]] -> dummyWhich[a && b, c, a, d]; rule4 = dummyWhich[a__, b_, dummyIf[c_, d_, e_], f___] -> dummyWhich[a, b && c, d, b, e, f]; rule5 = dummyWhich[a__, b_, dummyWhich[c_, d_, e___], f___] -> dummyWhich[a, b && c, d, b, dummyWhich[e], f]; rule6 = dummyWhich[] :> Null; rule7 = HoldPattern[True && a_] :> a; expr //. {rule1, rule2, rule3, rule4, rule5, rule6, rule7} /. {dummyWhich -> Which The result of that last line is: {Which[Aa == Ab, a Sa + b Sb, Ca == 0 || Cb == 0, a Sa + b Sb, True, SS], Which[Aa == Ab, a Ca + b Cb, Ca == 0 || Cb == 0, a Ca + b Cb, True, CC], Which[Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, X > 180, X - 180, X < 0, X + 180, True, X]} I've used Null where that would be the result of your original logic -- maybe you want it to be Nil. Or maybe you want to use Null instead of Nil. Your choice. Wherever you see Null, there's possibly a case you haven't covered. -----Original Message----- If[Aa == Ab, a Ca + b Cb, If[Ca == 0 || Cb == 0, a Ca + b Cb, CC]], If[Ca == 0 && Cb == 0, Nil, If[Z1 == Z2, If[Ca < 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} ==== Sorry... ONE MORE simplification: Which[ Ca == 0 == Cb, Nil, Z1 == Z2, Which[Ca < 0, Aa, Cb < 0, Ab, Ca > 0, Mod[Aa + 90, 180], Cb > 0, Mod[Ab + 90, 180]], Aa == Ab, Aa, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True, Mod[X, 180]] I didn't like repeating the Z1==Z2 test. Upon noticing this, I thought at first that my rule5 caused this situation to occur, but actually I have to delete both rule 4 AND rule5 to eliminate it, and the resulting expression wouldn't be as easy to understand. -----Original Message----- Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, X > 180, X - 180, X < 0, X + 180, True, X] Several things occurred to me in terms of simplifying this. The first was that the last six arguments of Which could be replaced with two, as follows: Which[ Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True,Mod[X,180]] This probably has no advantage other than clarity, but that's worth the others: rule8 = dummyIf[a_ > b_, a_ - b_, dummyIf[a_ < 0, a_ + b_, a_]] :> Mod[a, b]; but it had no effect on the expression. Apparently there was no match, and I can't see why. The second thing I noticed was that the condition Ca == 0 && Cb == 0 occurs twice in the Which statement. The following rule fixes that kind of situation: rule9 = dummyWhich[a___, b_, c_, d__, b_, e_, f__] /; EvenQ[Length[List[a]]] && EvenQ[Length[List[d]]] && EvenQ[Length[List[f]]] :> dummyWhich[a, b, c, d, f]; expr //. {rule1, rule2, rule3, rule4, rule5, rule6, rule7, rule9} The third thing I noticed was the treatment of Aa and Ab. The expressions If[Aa > 90, Aa - 90, Aa + 90] If[Ab > 90, Ab - 90, Ab + 90] in the original expression may be equivalent (depending on what's known about Aa and Ab a priori) to: Mod[Aa+90,180] Mod[Ab+90,180] If that is a valid assumption in the situation (no way for me to know), the expression is now Which[ Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, Z1 == Z2 && Ca > 0, Mod[Aa + 90, 180], Z1 == Z2 && Cb > 0, Mod[Ab + 90, 180], Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True, Mod[X, 180]] The next thing I notice is that the sixth condition, Z1==Z2, which results in Null if found True, isn't needed. If Z1==Z2 then one of the first five conditions is also true, so execution wouldn't get that far. Hence we can delete that condition-response pair. Finally, I'm stopping with this: Which[ Ca == 0 == Cb, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, Z1 == Z2 && Ca > 0, Mod[Aa + 90, 180], Z1 == Z2 && Cb > 0, Mod[Ab + 90, 180], Aa == Ab, Aa, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, True, Mod[X, 180]] I have no idea what you're actually doing with this code, but it looks weird even AFTER I've simplified it. -----Original Message----- 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[ Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} /. {If -> dummyIf}; rule1 = dummyIf[a_, b_, dummyIf[c_, d_, e_]] -> dummyWhich[a, b, c, d, True, e]; rule2 = dummyIf[a_, dummyIf[ b_, c_, d_], e_] -> dummyWhich[a && b, c, a, d, True, e]; rule3 = dummyIf[a_, dummyIf[b_, c_, d_]] -> dummyWhich[a && b, c, a, d]; rule4 = dummyWhich[a__, b_, dummyIf[c_, d_, e_], f___] -> dummyWhich[a, b && c, d, b, e, f]; rule5 = dummyWhich[a__, b_, dummyWhich[c_, d_, e___], f___] -> dummyWhich[a, b && c, d, b, dummyWhich[e], f]; rule6 = dummyWhich[] :> Null; rule7 = HoldPattern[True && a_] :> a; expr //. {rule1, rule2, rule3, rule4, rule5, rule6, rule7} /. {dummyWhich -> Which The result of that last line is: {Which[Aa == Ab, a Sa + b Sb, Ca == 0 || Cb == 0, a Sa + b Sb, True, SS], Which[Aa == Ab, a Ca + b Cb, Ca == 0 || Cb == 0, a Ca + b Cb, True, CC], Which[Ca == 0 && Cb == 0, Nil, Z1 == Z2 && Ca < 0, Aa, Z1 == Z2 && Cb < 0, Ab, (Z1 == Z2 && Ca > 0) && Aa > 90, Aa - 90, Z1 == Z2 && Ca > 0, Aa + 90, Z1 == Z2 && (Cb > 0 && Ab > 90), Ab - 90, Z1 == Z2 && Cb > 0, Ab + 90, Z1 == Z2, Null, Aa == Ab, Aa, Ca == 0 && Cb == 0, Nil, Ca == 0, Ab, Cb == 0, Aa, Y == 0, Nil, X > 180, X - 180, X < 0, X + 180, True, X]} I've used Null where that would be the result of your original logic -- maybe you want it to be Nil. Or maybe you want to use Null instead of Nil. Your choice. Wherever you see Null, there's possibly a case you haven't covered. -----Original Message----- If[Aa == Ab, a Ca + b Cb, If[Ca == 0 || Cb == 0, a Ca + b Cb, CC]], If[Ca == 0 && Cb == 0, Nil, If[Z1 == Z2, If[Ca < 0, Aa, If[Cb < 0, Ab, If[Ca > 0, If[Aa > 90, Aa - 90, Aa + 90], If[Cb > 0, If[Ab > 90, Ab - 90, Ab + 90]]]]], If[Aa == Ab, Aa, If[Ca == 0 && Cb == 0, Nil, If[Ca == 0, Ab, If[Cb == 0, Aa, If[Y == 0, Nil, If[X > 180, X - 180, If[X < 0, X + 180, X]]]]]]]]]} ==== If x1 through x20 aren't unknowns, you should have told me their values. If you can't, they are UNKNOWN. True, you're not solving for them... but they affect the dimensionality of the problem just as if you WERE solving for them. You have a 20-dimensional space full of contingencies -- solution forms that depend on the values of x1 through x20. The Solve function won't deal with contingencies even if it could, and if it tried, there would be too many. Bobby -----Original Message----- >Bobby -----Original Message----- >Sent: Saturday, September 07, 2002 1:54 AM >resend >Would someone with a very fast machine and lots of memory be willing to >try >this Solve for me? >It is the inverse of the 20 node quadratic hexahedral mapping used in >finite element analysis. >None of my computers can handle this - they run out of memory (using (*Hex20 Node definition in global coordinates *) >Clear[ >x1, y1, z1, >x2, y2, z2, >x3, y3, z3, >x4, y4, z4, >x5, y5, z5, >x6, y6, z6, >x7, y7, z7, >x8, y8, z8, >x9, y9, z9, >x10, y10, z10, >x11, y11, z11, >x12, y12, z12, >x13, y13, z13, >x14, y14, z14, >x15, y15, z15, >x16, y16, z16, >x17, y17, z17, >x18, y18, z18, >x19, y19, z19, >x20, y20, z20]; (* local coordinates *) >Clear[u, v, w]; (* Global co-ordinates *) >Clear[x, y, z]; (* corner nodes *) >N1= (1-u)*(1-v)*(1-w)*(-2-u-v-w)/8; >N3= (1+u)*(1-v)*(1-w)*(-2+u-v-w)/8; >N5= (1+u)*(1+v)*(1-w)*(-2+u+v-w)/8; >N7= (1-u)*(1+v)*(1-w)*(-2-u+v-w)/8; >N13=(1-u)*(1-v)*(1+w)*(-2-u-v+w)/8; >N15=(1+u)*(1-v)*(1+w)*(-2+u-v+w)/8; >N17=(1+u)*(1+v)*(1+w)*(-2+u+v+w)/8; >N19=(1-u)*(1+v)*(1+w)*(-2-u+v+w)/8; >(* to u nodes *) >N2= (1-u^2)*(1-v)*(1-w)/4; >N6= (1-u^2)*(1+v)*(1-w)/4; >N14=(1-u^2)*(1-v)*(1+w)/4; >N18=(1-u^2)*(1+v)*(1+w)/4; >(* to v nodes *) >N4= (1+u)*(1-v^2)*(1-w)/4; >N8= (1-u)*(1-v^2)*(1-w)/4; >N16=(1+u)*(1-v^2)*(1+w)/4; >N20=(1-u)*(1-v^2)*(1+w)/4; >(* to w nodes *) >N9= (1-u)*(1-v)*(1-w^2)/4; >N10=(1+u)*(1-v)*(1-w^2)/4; >N11=(1+u)*(1+v)*(1-w^2)/4; >N12=(1-u)*(1-v)*(1-w^2)/4; (* solve the inverse transform *) >Solve[{ >x1*N1+x2*N2+x3*N3+x4*N4+x5*N5+x6*N6+x7*N7+x8*N8+x9*N9+x10*N10+ >x11*N11+x12*N12+x13*N13+x14*N14+x15*N15+x16*N16+x17*N17+x18*N18+x19*N19 + >x20*N20-x==0, >y1*N1+y2*N2+y3*N3+y4*N4+y5*N5+y6*N6+y7*N7+y8*N8+y9*N9+y10*N10+ >y11*N11+y12*N12+y13*N13+y14*N14+y15*N15+y16*N16+y17*N17+y18*N18+y19*N19 + >y20*N20-y==0, >z1*N1+z2*N2+z3*N3+z4*N4+z5*N5+z6*N6+z7*N7+z8*N8+z9*N9+z10*N10+ >z11*N11+z12*N12+z13*N13+z14*N14+z15*N15+z16*N16+z17*N17+z18*N18+z19*N19 + >z20*N20-z==0}, >{u,v,w}] >Christopher J. Purcell >Defence R&D Canada - Atlantic >9 Grove St., PO Box 1012 >Dartmouth NS Canada B2Y 3Z7 Christopher J. Purcell Defence R&D Canada - Atlantic 9 Grove St., PO Box 1012 Dartmouth NS Canada B2Y 3Z7 ==== If we do think of x1...x20 as constants, it's still fourth-order, with three simultaneous equations. Evaluate this: Solve[a x^4 + b x^3 + c x^2 + d x + e == 0, x] and look at the four solutions found -- for a single variable and one equation! For your problem, I'm thinking there'd be at least 4^3 solutions, each of them even more complicated than these, and an unknown number of contingencies to deal with, depending on values of the xi. Check each radical for a negative argument, and the special cases multiply quickly. For some values of the xi, there'd be no solutions; for other values, many solutions. If you do manage to solve this, I'd love to see the method! Bobby -----Original Message----- many. Bobby -----Original Message----- >Bobby -----Original Message----- >Sent: Saturday, September 07, 2002 1:54 AM >resend >Would someone with a very fast machine and lots of memory be willing to >try >this Solve for me? >It is the inverse of the 20 node quadratic hexahedral mapping used in >finite element analysis. >None of my computers can handle this - they run out of memory (using (*Hex20 Node definition in global coordinates *) >Clear[ >x1, y1, z1, >x2, y2, z2, >x3, y3, z3, >x4, y4, z4, >x5, y5, z5, >x6, y6, z6, >x7, y7, z7, >x8, y8, z8, >x9, y9, z9, >x10, y10, z10, >x11, y11, z11, >x12, y12, z12, >x13, y13, z13, >x14, y14, z14, >x15, y15, z15, >x16, y16, z16, >x17, y17, z17, >x18, y18, z18, >x19, y19, z19, >x20, y20, z20]; (* local coordinates *) >Clear[u, v, w]; (* Global co-ordinates *) >Clear[x, y, z]; (* corner nodes *) >N1= (1-u)*(1-v)*(1-w)*(-2-u-v-w)/8; >N3= (1+u)*(1-v)*(1-w)*(-2+u-v-w)/8; >N5= (1+u)*(1+v)*(1-w)*(-2+u+v-w)/8; >N7= (1-u)*(1+v)*(1-w)*(-2-u+v-w)/8; >N13=(1-u)*(1-v)*(1+w)*(-2-u-v+w)/8; >N15=(1+u)*(1-v)*(1+w)*(-2+u-v+w)/8; >N17=(1+u)*(1+v)*(1+w)*(-2+u+v+w)/8; >N19=(1-u)*(1+v)*(1+w)*(-2-u+v+w)/8; >(* to u nodes *) >N2= (1-u^2)*(1-v)*(1-w)/4; >N6= (1-u^2)*(1+v)*(1-w)/4; >N14=(1-u^2)*(1-v)*(1+w)/4; >N18=(1-u^2)*(1+v)*(1+w)/4; >(* to v nodes *) >N4= (1+u)*(1-v^2)*(1-w)/4; >N8= (1-u)*(1-v^2)*(1-w)/4; >N16=(1+u)*(1-v^2)*(1+w)/4; >N20=(1-u)*(1-v^2)*(1+w)/4; >(* to w nodes *) >N9= (1-u)*(1-v)*(1-w^2)/4; >N10=(1+u)*(1-v)*(1-w^2)/4; >N11=(1+u)*(1+v)*(1-w^2)/4; >N12=(1-u)*(1-v)*(1-w^2)/4; (* solve the inverse transform *) >Solve[{ >x1*N1+x2*N2+x3*N3+x4*N4+x5*N5+x6*N6+x7*N7+x8*N8+x9*N9+x10*N10+ >x11*N11+x12*N12+x13*N13+x14*N14+x15*N15+x16*N16+x17*N17+x18*N18+x19*N19 + >x20*N20-x==0, >y1*N1+y2*N2+y3*N3+y4*N4+y5*N5+y6*N6+y7*N7+y8*N8+y9*N9+y10*N10+ >y11*N11+y12*N12+y13*N13+y14*N14+y15*N15+y16*N16+y17*N17+y18*N18+y19*N19 + >y20*N20-y==0, >z1*N1+z2*N2+z3*N3+z4*N4+z5*N5+z6*N6+z7*N7+z8*N8+z9*N9+z10*N10+ >z11*N11+z12*N12+z13*N13+z14*N14+z15*N15+z16*N16+z17*N17+z18*N18+z19*N19 + >z20*N20-z==0}, >{u,v,w}] >Christopher J. Purcell >Defence R&D Canada - Atlantic >9 Grove St., PO Box 1012 >Dartmouth NS Canada B2Y 3Z7 Christopher J. Purcell Defence R&D Canada - Atlantic 9 Grove St., PO Box 1012 Dartmouth NS Canada B2Y 3Z7 ==== The documentation you're talking about doesn't even deserve to be called documentation. A few hints; that's it. -----Original Message----- > unable to find a complete guide to word processing with mathematica. > Does anyone know where such a document could be found? ==== This is a question that's been asked and answered a number of times. One answer, by me, can be found by going to www.wolfram.com and selecting first Resource Library and then MathGroup on the Resource Library page; then search for direction field. You'll find, among many others, the URL: http://library.wolfram.com/mathgroup/archive/2002/Jul/msg00163.html (This is not to chastise you for not first looking there, but to suggest a way for you to find answers more efficiently than waiting for the ou need to load the package, then manufacture a vector field whose plot is the slope field (what you call vector flow diagram). To do the latter, let me separately define the function giving the right-hand side of the ODE; you don't really have to do that, but could instead directly use as the first argument to PlotVectorField the right-hand side of my definition of f. << Graphics`PlotField` f[t_, y_] := {1, 1/y} field = PlotVectorField[f[t, y], {t, 0, 3}, {y, 0.05, 3}]; I gave a name to the graphics result just in order to re-display the field below along with the graph of one solution. First, find a general solution: soln = First @ DSolve[{y'[t] == 1/y[t], y[0] == 1}, y[t], t] {y[t] -> Sqrt[2]*Sqrt[1/2 + t]} grf = Plot[y[t] /. soln, {t, 0, 3}]; Now combine the plot with the field: Show[grf, field]; Note that PlotVectorField does not always give quite the result you might want, since the vector lengths are scaled with respect to the magnitude of the actual vectors. You might want just a direction field in which all arrows have the same length. You can produce such a graphic by using the ScaleFunction and ScaleFactor options to PlotVectorField. For example: field = PlotVectorField[f[t, y], {t, 0, 3}, {y, 0.05, 3}, ScaleFunction -> (1 &), ScaleFactor -> 0.175]; ... I have several differential equations I would like to plot in Mathematica. I would like to plot the Slope Fields of them though. Can anyone lead me in the right direction? I can solve the equations trivially but I want to display the slope fields. An example follows : y' + 2y = 3 ... -- Murray Eisenberg murray@math.umass.edu Mathematics & Statistics Dept. Lederle Graduate Research Tower phone 413 549-1020 (H) University of Massachusetts 413 545-2859 (W) 710 North Pleasant Street Amherst, MA 01375 ==== I am trying to implement a very simple sorted tree to quickly store some real numbers I need. I have written an add, delete, minimum, and pop (delete the lowest value) function and they seem to work ok but are very slow. Let's just look @ my implementation of the add part: nums=Null;(*my initial blank Tree) In[326]:= Clear[add] In[327]:= add[Null,x_Real]:=node[x,Null,Null] add[Null,node[x_Real,lower_,higher_]]=node[x,lower,higher] In[328]:= add[node[x_Real,lower_,higher_],y_Real]:= If[x>y,node[x,add[lower,y],higher],node[x,lower,add[higher,y]]] In[288]:= add[node[x_Real,lowerx_,higherx_],node[y_Real,lowery_,highery_]]:=If[x>y, node[x,add[lowerx,node[y,lowery,highery]],higherx], node[x,lowerx,add[higherx,node[y,lowery,highery]]] ] Now this is my attempt to test how fast my add works: SeedRandom[5]; Do[nums=add[nums,Random[]],{5000}];//Timing Out[333]= {13.279 Second,Null} running on an 1.4GHz Athlon with 1GB of ram). Questions: 1. Is this as fast as I can get my code to run? 2. Am I doing something obviously stupid? 3. would Compiling things help? Husain ==== The vectors that NullSpace gives do not need to be orthogonal to each other. What about the following variant, that adds one vector at a time? It may not be fast enough for serious number crunching. OrthogonalComplement[v__ /; MatrixQ[{v}]] := With[{n = Length[NullSpace[{v}]]}, Take[Nest[Join[{First[NullSpace[#]]}, #] &, {v}, n], n]]; OrthonormalComplement[v__ /; MatrixQ[{v}]] := Map[#/Sqrt[#.#] &, OrthogonalComplement[v]]; When the input of OrthogonalComplement is integer, I expect the output to be integer too: OrthogonalComplement[{1, 2, 3}] {{1, -5, 3}, {-3, 0, 1}} Gianluca Gorni OrthogonalUnitVectors that I sent a few minutes ago doesn't do what I > thought it would. > After making either definition below I get the following: In[2]:= > s1=OrthogonalUnitVectors[{1,0,1/2,1,0},{0,1,-1,1/2,2}] Out[2]= > {{0, -2/Sqrt[5], 0, 0, 1/Sqrt[5]}, {-2/3, -1/3, 0, 2/3, 0}, {-1/3, 2/3, > 2/3, 0, 0}} > The dot products below aren't zero, so the vectors aren't orthogonal. What > went wrong? > In[3]:= > Part[s1,1].Part[s1,2] Out[3]= > 2/(3*Sqrt[5]) > In[4]:= > Part[s1,1].Part[s1,3] Out[3]= > -4/(3*Sqrt[5]) > -------------- >> Hugh Goyder and Park gave a most elegant function to find two >> vectors that are orthogonal to one vector in 3D. The key to coming up >> with the elegant solution is an understanding of Mathematica's NullSpace >> function. We can easily make the version from Hugh and much more >> general with the version below. >> ------------- >> OrthogonalUnitVectors[vect__?VectorQ]:= >> #/Sqrt[#.#]&/@NullSpace[{vect}] >> ------------- >> The version above will give a set of unit orthogonal vectors if given any >> number of vectors in any dimension. >> So besides giving it a 3D vector we can give it the following: >> OrthogonalUnitVectors[{2,1,0,-1,1}] >> or >> OrthogonalUnitVectors[{0,1,0,1/2,1},{1,0,-1,1/2}] >> ------------ >> But the short version above isn't very robust. >> (1) Clear[x,y,z];NullSpace[{{x,y,z}}] >> returns two vectors orthogonal to {x,y,z}, but the two vectors >> NullSpace returns aren't orthogonal to each other. >> So (OrthogonalUnitVectors) should only work with numeric vectors. >> (2) We should ensure all the vectors have the same dimension and length > 1. >> I give a less concise version below that corrects these problems. >> ------------ >> OrthogonalUnitVectors[vect__?(VectorQ[#,NumericQ]&)]/; >> (SameQ@@Length/@{vect})&&(Length[First[{vect}]]>1):= >> #/Sqrt[#.#]&/@NullSpace[{vect}] >> -------------- >> Ted Ersek >> Get Mathematica tips, tricks from >> http://www.verbeia.com/mathematica/tips/Tricks.html > ==== I am using MultipleListPlot to plot a range of 2D graphs and have found that the last graph, which has a much greater (by a factor of 4) x range than the others does not plot completely but is truncated in the x direction. By removing the Frame, I can see the points where these lie outside the Frame but these are not joined etc. How do I get over this? - Malcolm Woodruff Reply-To: kuska@informatik.uni-leipzig.de ==== set the PlotRange explicit PlotRange->{{0,largeX},Auatic} I am using MultipleListPlot to plot a range of 2D graphs and have found that > the last graph, which has a much greater (by a factor of 4) x range than the > others does not plot completely but is truncated in the x direction. By > removing the Frame, I can see the points where these lie outside the Frame > but these are not joined etc. How do I get over this? > - > Malcolm Woodruff ==== I noticed the same thing when I attempted to plot data with disjoint x-ranges. One way is not to use the MultipleListPlot[] function. You can plot both the lines and the markers by applying Line and Point graphics to your datasets. And use Show[] to display everything on one plot. You can use the PlotSymbol[] function instead of Point. The PlotSymbol function is available in the MultipleListPlot package to get the diamond, star, square, etc. Another way is to modify the MultipleListPlot[] function. I modified the function to do multiple plots on a polar grid. WARNING - I advise you not to change any of the standpackages. Instead, first make a copy then modify it. One of the modifications I made was to change 'pts = First[Transpose[data]]' in the handleset function to 'pts=data'. This seems to work for my application. Hope this helps. Lawrence > I am using MultipleListPlot to plot a range of 2D graphs and have found that > the last graph, which has a much greater (by a factor of 4) x range than the > others does not plot completely but is truncated in the x direction. By > removing the Frame, I can see the points where these lie outside the Frame > but these are not joined etc. How do I get over this? > - > Malcolm Woodruff ==== Group, Sorry for the previous posting. The following solves (to some extent my problem) < (1) Is there a way in Mathematica 4.2 to put two separate gifs into a > single, cell side by side (with some intervening space), without > having to combine them in some graphics program first? In particular, I'd like to do that within a text cell. Even in a new Input cell, if I first create a GridBox (via Inut>Create > Table/Matrix/Palette) and then try to insert the first gif (via > Edit>Insert Object>Create from File ....), Mathematica promptly > crashes. (Mathematica 4.2 under Windows 2000.) You could do something like this: grlist = Map[Import[#, GIF]&, {, }] celist = Map[Cell[GraphicsData[PostScript, DisplayString[#]]]&, grlist] CellPrint[Cell[BoxData[GridBox[{celist}]], Text]] This embeds the PostScript versions of these graphics each as inline cells in a GridBox[] structure. > (2) Is there a way to cause a gif imported into a Mathematica 4.2 > notebook to become a hyperlink -- so that when the user clicks on the > gif the hyperlink's target is summoned? Again, we embed the PostScript graphic within an inline cell. Here is a function that does what you describe: createGIFHyperlink[gifImage_String, url_String] := NotebookWrite[EvaluationNotebook[], Cell[BoxData[ ButtonBox[ Cell[GraphicsData[PostScript, DisplayString[Import[gifImage, GIF]]]], ButtonStyle -> Hyperlink, ButtonData :> {URL[url], None}]], Text]] -- User Interface Programmer paulh@wolfram.com Wolfram Research, Inc. ==== > 1.) How to find the General Solution for below's partial differential > equation? > (y + u) du/dx + y (du/dy) = x - y > ** I use d to represent the partial differential symbol. > Can it be solved by function NDSolve in mathematica 4.1? How? It's been a while, but I suspect that the presence of the nonlinear convective term u du/dx alone makes a general solution unlikely. It's easy to get a numerical solution for a particular set of boundary conditions, but you are not guaranteed a solution. For instance, the following choice of a and b are skirting failure, as the diagnostic contour plot shows. Burton Needs[Graphics`Colors`] !((soln = With[{a = (-1), b = 1.27, e = 0.001}, solutions = NDSolve[{((y + u[x, y])) [PartialD]_x u[x, y] + y [PartialD]_y u[x, y] == x - y, u[x, b] == 0, u[a, y] == 0}, u, {x, a, 10}, {y, b, 10}]; [IndentingNewLine]Plot3D[ Evaluate[u[x, y] /. [InvisibleSpace]First[solutions]], {x, a, 10}, {y, b, 10}, PlotPoints [Rule] 50]; [IndentingNewLine]ContourPlot[ Evaluate[y + u[x, y] /. [InvisibleSpace]First[solutions]], {x, a, 10}, {y, b, 10}, PlotPoints [Rule] 50, ColorFunction [Rule] ((If[#1 > e, Green, Red] &))][IndentingNewLine]];)) ==== I'd be happier if I could just control pagination, or if auatic pagination were done sensibly. Yes, I can make a page break where I want one, but I can't eliminate page-breaks that leave tons of white-space at the end of some pages. In some of my notebooks, Mathematica puts the first cell on a page by itself, using less than an inch --- and NOT because the next cell takes a whole page. If I force the first two cells to stay together, the problem just moves to the next page. It's enough to make Mathematica useless as a word processor all by itself. -----Original Message----- derivations; - setting typefaces in tables of material. I figure most of this is do-able, but I don't have the time, or patience, to spend too much time on it. So, I'll be the first cuser when you write the guide to math publishing in Mathematica - I just hope you 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 to find a complete guide to word processing with mathematica. > Does anyone know where such a document could be found? ==== Bobby, Have a look at the Glue palette found here: --- > I'd be happier if I could just control pagination, or if auatic > pagination were done sensibly. Yes, I can make a page break where I want one, but I can't eliminate > page-breaks that leave tons of white-space at the end of some pages. In > some of my notebooks, Mathematica puts the first cell on a page by > itself, using less than an inch --- and NOT because the next cell takes > a whole page. If I force the first two cells to stay together, the > problem just moves to the next page. It's enough to make Mathematica useless as a word processor all by > itself. -----Original Message----- > Kenny, Sympathy but no solution. I too have been trying to use Mathematica (v4.2 most recently) to type > maths papers and the like but I'm not ready to ditch LaTeX yet. There > are just too many cases where I cannot figure out how to achieve what I > want in Mathematica, things like: > - left brackets spanning multiple lines for defining hybrid functions; > - vertical alignment of equals signs in multi-line equations or > derivations; > - setting typefaces in tables of material. I figure most of this is do-able, but I don't have the time, or > patience, to spend too much time on it. So, I'll be the first cuser > when you write the guide to math publishing in Mathematica - I just hope > you 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 to find a complete guide to word processing with mathematica. >>Does anyone know where such a document could be found? > ==== Wouldn't it be nice if NullSpace's behavior were DOCUMENTED? Otherwise, > it's futile to give it approximate numbers expecting any particular > behavior. Even if it always works, it may not work in the next version > of Mathematica. Bobby The expected, and documented, behavior is that the output should be a basis for the null space, that is, solutions of the homogeneous matrix equation A.x==0. If this were to stop working then that would be a serious bug. Is this the behavior you mean? The implementation notes of the manual mention that approximate NullSpace is based on a singular values decomposition. This in fact gives resulting vectors that are orthonormal by the usual conjugate-symmetric inner product on C (though these are now not normal to the original vector in this same inner product, unless they are real-valued). But this basis-orthogonality is not part of the mission of NullSpace and moreover should not become part of it. Hence that particular (and implementation dependent) aspect of NullSpace should not become documented. Daniel Lichtblau Wolfram Research ==== Wouldn't it be nice if NullSpace's behavior were DOCUMENTED? Otherwise, it's futile to give it approximate numbers expecting any particular behavior. Even if it always works, it may not work in the next version of Mathematica. Bobby -----Original Message----- enough for me. Therefore I modify Ted's routine to OrthogonalUnitVectors[vect__?(VectorQ[#, NumericQ] &)] /; (SameQ @@ Length /@ {vect}) && (Length[First[{vect}]] > 1) := #/Sqrt[#.#] & /@ NullSpace[{vect}// N] ---------------- Lets see what NullSpace does with approximate complex vectors. In[1]:= v1 = {1.0 I, 0.0, 0.5 I, 0.0, 1.0}; v2 = {0.0, 2.0, 1.0 I, 2.0, 0.5}; {v3,v4,v5} = NullSpace[{v1,v2}] Out[3]= {{-0.730153 + 0.*I, 0. - 0.138254*I, 0.250585 + 0.*I, 0. - 0.138254*I, 0. + 0.60486*I}, {0. + 0.*I, -0.515861 + 0.*I, 0. + 0.457321*I, 0.687357 + 0.*I, 0.22866 + 0.*I}, {0. + 0.*I, 0.510406 + 0.*I, 0. + 0.740442*I, -0.23274 + 0.*I, 0.370221 + 0.*I}} -------- In the next line we see NullSpace returned vectors that are orthogonal to the vectors we gave NullSpace. In[4]:= {v1.v3, v1.v4, v1.v5, v2.v3, v2.v4, v2.v5}//Chop Out[4]= {0, 0, 0, 0, 0, 0} ---------- However, the vectors returned aren't orthogonal to each other. In[5]:= {v3.v4, v3.v5, v4.v5}//Chop Out[5]= {0.229195*I, 0.371087*I, -0.677239} --------- I suppose an OrthogonalUnitVectors function that uses NullSpace should (1) Only accept real valued vectors. (2) Ensure NullSpace is given approximate vectors. ------ Ted Ersek ==== Plot[{x,x^2},{x,0,10}] as Garza Mexico City ----- Original Message ----- ==== Mario, Here is a fancy version of your plot. I used Text statements within an Epilog option to label the two curves. The regular plot statement allows you to plot a series of functions inclosed in a list. Needs[Graphics`Colors`] Plot[{x, x^2}, {x, 0, 1}, PlotStyle -> {Black, Blue}, Frame -> True, FrameLabel -> {x, y}, PlotLabel -> Comparison of Two Functions, Epilog -> {Text[x, {0.5, 0.55}], Blue, Text[x^2, {0.7, 0.4}]}, Background -> Linen, ImageSize -> 500]; Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.net Approved: Steven M. Christensen , Moderator ==== >How can I plot with Mathematica two function in the same graphic? I >explain better: I'd like to draw a) y=x and b) y=x^2 in the same axes Plot[{x, x^2},{x,-2,2}] ==== This might be a convenient way of defining color functions: ClearAll[cfun] cfun[colors_List, brkPts_List] /; Length@colors == Length@brkPts := Function[z, Which @@ Sequence@Flatten@Transpose[{Less[z, #] & /@ brkPts, colors}]] colors = {White, RoyalBlue, White, Red, White}; brkPts = {-0.6, -0.42, 0.4, 0.6, Infinity}; ContourPlot[f[x, y], { x, 0, Pi}, {y, 0, Pi}, PlotPoints -> 30, ColorFunctionScaling -> False, ColorFunction -> cfun[colors, brkPts], Contours -> contourvalues]; colors = {Yellow, Peru, Salmon, Apricot, HotPink, Linen}; brkPts = {-0.6, -0.42, 0.2, 0.4, 0.6, Infinity}; Timing[ContourPlot[f[x, y], {x, 0, Pi}, {y, 0, Pi}, PlotPoints -> 30, ColorFunctionScaling -> False, ColorFunction -> cfun[colors, brkPts], Contours -> contourvalues];] -----Original Message----- it is difficult to obtain in this plot. contourvalues = Complement[Range[-1, 1, 0.2], {0.}] {-1, -0.8, -0.6, -0.4, -0.2, 0.2, 0.4, 0.6, 0.8, 1.} Now we define a ColorFunction for the plot. I actually colored two different bands to show how you can make a general color function to give each band a desired color. cfun[z_] := Which[ -0.6 < z < -0.42, RoyalBlue, 0.4 < z < 0.6, Red, True, White] ContourPlot[f[x, y], {x, 0, Pi}, {y, 0, Pi}, PlotPoints -> 30, ColorFunctionScaling -> False, ColorFunction -> cfun, Contours -> contourvalues]; Using the option ColorFunctionScaling -> False says that the z value will be the actual value of f[x,y]. Otherwise, in general, it will be scaled between 0 and 1. It is easier to write a color function when z is the actual value of the function. Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.net Approved: Steven M. Christensen , Moderator ==== For what it's worth, here's a version of Simpson's rule that works with both equally spaced and unequally spaced points. simp = Compile[{x1,y1,x2,y2,x3,y3}, With[{h13=x1-x3, h12=x1-x2, h23=x2-x3}, (-y1*h23*(2x1-3x2+x3) - y2*h13^2 + y3*h12*(x1-3x2+2x3))*h13/(6h12*h23)]]; ListSimpson[vals_?VectorQ, dx_Real]:= With[{n=Length[vals]}, If[OddQ[n], (dx/3)*(First[vals] + Last[vals] + 4*Plus@@vals[[Range[2,n-1,2]]] + 2*Plus@@vals[[Range[3,n-2,2]]]), $Failed]]; ListSimpson[datapts_?MatrixQ] := If[OddQ[Length[datapts]], Plus@@(simp@@Flatten[#]&/@ Partition[datapts,3,2]), $Failed] (* equally spaced points*) data = Table[100*Sin[x], {x, 0, 100, 0.001}]; ListSimpson[data, .001] // Timing {0.21 Second, 13.7681} (* randomly spaced points *) Table[ {data = Sort[Table[With[{x = Random[Real,{0,100}]},{x,100*Sin[x]}],{100001}]]; ListSimpson[data] // Timing, With[{xvals = Transpose[data][[1]]}, gaps = Drop[RotateLeft[xvals] - xvals, -1]; {Min[gaps], Max[gaps]}]}, {4}] {{{1.62 Second, 13.8348}, {1.32223*10^-8, 0.0114454}}, {{1.59 Second, 13.8261}, {2.42163*10^-8, 0.0127803}}, {{1.62 Second, 13.8455}, {1.08923*10^-9, 0.0119059}}, {{1.54 Second, 13.7987}, {9.66814*10^-9, 0.0106988}}} --- > Mathew, Some possibilities < Integrate[ > Interpolation[data, InterpolationOrder[Rule]1][x], > {x,0,100}]//Timing {4.56 Second,13.768} Trapezium rule with equal steps: #[[1]]+#[[-1]]+ 2 Tr[Take[#,{2,-2}]]&[data[[All,2]]] 0.01/2//Timing {0.22 Second,13.768} Trapezium rule with possibly unequal steps (Drop[#1,1] - Drop[#1,-1]).(Drop[#2,-1] + Drop[#2,1])&[ > data[[All,1]], data[[All,2]]]/2//Timing {0.83 Second,13.768} -- > Allan --------------------- > Allan Hayes > Mathematica Training and Consulting > Leicester UK > www.haystack.demon.co.uk > hay@haystack.demon.co.uk > Voice: +44 (0)116 271 4198 I've tracked down the slow operation of my Mathematica simulation > code to lie in the ListIntegrate command: G[n_] := ListIntegrate[xsec Phi[n], 1] where both xsec and Phi[n] are 400 values long. Is there a way to speed up ListIntegrate via Compile or a similar > technique? Matt --- Matthew Rosen Harvard-Smithsonian Center for Astrophysics Mail Stop 59 60 Garden Street Cambridge, MA 02138 e: mrosen@cfa.harvard.edu o: (617) 496-7614 ==== P.J., I'd be very interested in how you created these fancy cells. Do you open a cell, push Ctrl-Shift-e, type in all that text, then push Ctrl-Shift-e again? -----Original Message----- 1) Put your function braches in the rows of a grid box structure. 2) Add the following options to your cell: ShowAutoStyles -> False SpanMaxSize -> Infinity The following cell snippet demonstrates how this influences the result. To view it, paste the Cell[] expression into a notebook and then click on Yes when you are prompted on whether the front end should interpret the result. Cell[BoxData[ FormBox[ RowBox[{ RowBox[{f, (, x, )}], =, RowBox[{{, GridBox[{ {x, RowBox[{x, , <, 0}]}, { SuperscriptBox[x, 2], RowBox[{0, [LessEqual], x, <, 1}]}, { RowBox[{sin, (, x, )}], RowBox[{1, [LessEqual], x, <, 2}]}, { RowBox[{[CapitalGamma], (, x, )}], RowBox[{x, [GreaterEqual], 2}]} }]}]}], TraditionalForm]], DisplayFormula, ShowAutoStyles->False, SpanMaxSize->Infinity] > - vertical alignment of equals signs in multi-line equations or > derivations; Put your equations in a GridBox and set the ColumnAlignments option to a string containing the equal sign. Cell[BoxData[ FormBox[GridBox[{ { RowBox[{ RowBox[{ RowBox[{3, x}], , +, , RowBox[{4, , y}]}], , =, , 9}]}, { RowBox[{ RowBox[{ RowBox[{2, x}], , -, , RowBox[{7, , y}]}], =, RowBox[{32, , -, , RowBox[{sin, (, x, )}]}]}]} }], TraditionalForm]], DisplayFormula, GridBoxOptions->{ColumnAlignments->{=}}] > - setting typefaces in tables of material. I think the Author Tools material that comes with Mathematica 4.2 might be able to help you do this. -- User Interface Programmer paulh@wolfram.com Wolfram Research, Inc. ==== > I'd be very interested in how you created these fancy cells. Do you > open a cell, push Ctrl-Shift-e, type in all that text, then push > Ctrl-Shift-e again? Surely you have grossly overestimated my expertise with Mathematica typesetting language syntax. :-) Here is the walkthrough for creating each sample cell. Example 1: 1) Set down a cell insertion point in your notebook. 3) Click on the cell's bracket. 4) Click on the front end menu command sequence: Format -> Cell Style -> DisplayFormula 5) Click on the front end menu command sequence: Cell -> Display As -> TraditionalForm 6) Click on the front end menu command: Format -> Option Inspector... 7) Make sure that the Option Inspector scope indicator is set to selection. 8) Enter the name of the option ShowAutoStyles into the search field and press the Lookup button. 9) Change the value of the option ShowAutoStyles from True to False. This prevents the front end from showing the left brace with the umatched syntax coloring. 10) Look up the option SpanMaxSize as was done in step (8). 11) Change the value of this option to Infinity. 12) Close the Option Inspector dialog. 13) Click within the cell you just created to create an editor caret. 14) Type in the text f(x) = { 4 x 2. 17) Enter in the function branch definitions in the left column, and the domains of definition in the right column. You can use the Tab key to navigate between grid elements. Example 2: 1) Repeat Steps (1) - (11) in Example 1. 2) With the cell bracket still selected, look up the option named ColumnAlignments. Change it from {Center} to {=}. 3) Close the Option Inspector dialog. 4) Click within the cell you just created to create an editor caret. 6) Enter in the equations, with one equation per grid element. I used the front end menu command Edit -> Copy As -> Cell Expression to get the underlying expressions. Hope that demystifies the procedure for you. -- User Interface Programmer paulh@wolfram.com Wolfram Research, Inc. ==== You're not doing anything dumb as far as I can see, but you're using far more memory than necessary. Here are statistics for your algorithm on my machine: SeedRandom[5]; nums = Null Do[nums = add[nums, Random[]], {5000}]; // Timing nums // LeafCount nums // ByteCount Count[nums, Null, Infinity] Count[nums, node[_, Null, _], Infinity] Count[nums, node[_, _, Null], Infinity] Count[nums, node[_, Null, Null], Infinity] Count[nums, node[___], Infinity] Count[nums, node[_, _node | _?NumericQ, _node | _?NumericQ], Infinity] Depth[nums] {4.406000000000001*Second, Null} 15001 240000 5001 (* Null values in the tree *) 2458 (* nodes without left offspring *) 2543 (* nodes without right offspring *) 1669 (* leaf nodes *) 4999 (* total nodes *) 1667 (* nodes with both left and right offspring *) 28 (* minus one, makes 27 levels in the tree *) It's clearly not ideal to store 5001 Nulls for 5000 actual values! Here's a method that doesn't change the algorithm much, but changes the storage method a great deal: First of all, I'm lazy, so I'll redefine <,>,<=,>=, etc. to make the code simpler: Unprotect[Less, Greater, LessEqual, GreaterEqual]; Less[a : _node ..] := Less @@ (First /@ {a}) LessEqual[a : _node ..] := LessEqual @@ (First /@ {a}) Greater[a : _node ..] := Greater @@ (First /@ {a}) GreaterEqual[a : _node ..] := GreaterEqual @@ (First /@ {a}) Protect[Less, Greater, LessEqual, GreaterEqual]; Again because I'm lazy, I'll define node so that I don't have to consciously avoid leaf nodes and unnecessary nesting: ClearAll[node] node[a___, node[b_], c___] = node[a, b, c]; node[node[a__]] = node[a]; Here's my add function: Clear[add] add[Null, x_] = node[x]; add[x_, (y_)?NumericQ] := add[node[x], node[y]] add[(x_)?NumericQ, y_] := add[node[x], node[y]] add[node[x_], node[y_]] := If[x > y, node[x, y], node[y, x]] (* <-- THERE! *) add[node[x_, lower_], y_node] := Which[y >= node[x], node[x, lower, y], y <= node[lower], node[lower, y, x], True, node[y, lower, x]] add[node[x_, lower_, higher_], y_node] /; node[y] <= node[x] := node[x, add[lower, y], higher] add[node[x_, lower_, higher_], y_node] := node[x, lower, add[higher, y]] Where I have a pointer is the code that prevents adding right-offspring to a leaf node. That avoids nodes like your node[a,Null,c]. When you would have node[a,b,Null], I have node[a,b]. When you would have node[a,Null,Null], I use a itself. Timing and space requirements are much better now: SeedRandom[5]; nums = Null Timing[Do[nums = add[nums, Random[]], {5000}]; ] {2.109*Second, Null} LeafCount[nums] ByteCount[nums] Count[nums, Null, Infinity] Count[nums, node[_], Infinity] Count[nums, node[_, _], Infinity] Count[nums, node[_, _, _], Infinity] Count[nums, node[___], Infinity] Depth[nums] 7852 (* 48% fewer leaf expressions, mostly due to Nulls and right-offspring eliminated *) 165624 (* 31% less storage *) 0 (* no Nulls, versus 5001 *) 0 (* no leaf nodes, versus 1669 *) (* no nodes without left-offspring, versus 2458 *) 705 (* nodes without right offspring, versus 2543 *) 2146 (* nodes with both left and right offspring, versus 1667 *) 2851 (* total nodes *) 21 (* 21 tree levels versus 27 *) The logical structure is identical to yours except that there are no nodes with only right-offspring. This incidentally tends to reduce tree depth. The storage format is far different: there are no trivial (childless) nodes, and nodes with left-offspring but no right-offspring are stored without a Null on the right. I could make this quite a bit faster, I think, if I spent time on the code to eliminate changes to Less, Greater, etc. and the frequent nesting followed by unnesting that goes on in the algorithm. It might take twice as much code, however, so I like it as is. You're probably better off storing these things in a heap, of course. Or -- better yet -- use the built-in Sort and Ordering functions. -----Original Message----- In[327]:= add[Null,x_Real]:=node[x,Null,Null] add[Null,node[x_Real,lower_,higher_]]=node[x,lower,higher] In[328]:= add[node[x_Real,lower_,higher_],y_Real]:= If[x>y,node[x,add[lower,y],higher],node[x,lower,add[higher,y]]] In[288]:= add[node[x_Real,lowerx_,higherx_],node[y_Real,lowery_,highery_]]:=If[x>y , node[x,add[lowerx,node[y,lowery,highery]],higherx], node[x,lowerx,add[higherx,node[y,lowery,highery]]] ] Now this is my attempt to test how fast my add works: SeedRandom[5]; Do[nums=add[nums,Random[]],{5000}];//Timing Out[333]= {13.279 Second,Null} RH7.3 running on an 1.4GHz Athlon with 1GB of ram). Questions: 1. Is this as fast as I can get my code to run? 2. Am I doing something obviously stupid? 3. would Compiling things help? Husain ==== > I am using MultipleListPlot to plot a range of 2D graphs and > have found that > the last graph, which has a much greater (by a factor of 4) x > range than the > others does not plot completely but is truncated in the x > direction. By > removing the Frame, I can see the points where these lie > outside the Frame > but these are not joined etc. How do I get over this? The easiest way is to expand your PlotRange to include the errant points. Dave. ==== Probably a rather simple question: what is the easiest way to open a notebook, execute it, then quit, from the command line? I would like to be able to do this from a Makefile. Sidney Reply-To: kuska@informatik.uni-leipzig.de ==== what may be in the *.m files? Mathematica commands ? that you can load into the kernel ? Makr all your input cells in the notebook as initialization, save it as an *.m file, add a Quit[] as the last command to the *.m file and run it with math << yourJustGeneratedMFile Regardds > Probably a rather simple question: what is the easiest way to open a > notebook, execute it, then quit, from the command line? I would like to be > able to do this from a Makefile. > Sidney ==== I made the algorithm more than twice as fast with very little change, while eliminating modifications to <,>,<=, and >=. Most of the improvement, in fact, was due to eliminating those changes AFTER the code was changed so they wouldn't be used anyway. That's because the comparisons appear everywhere, so adding rules to those functions really slows down pattern matching. Clear[add, node] node[a___, node[b_], c___] = node[a, b, c]; node[node[a__]] = node[a]; add[Null, x_] = node[x]; add[x_, (y_)?NumericQ] := add[node[x], node[y]] add[(x_)?NumericQ, y_] := add[node[x], node[y]] add[node[x_], node[y_]] := If[x >= y, node[x, y], node[y, x]] add[node[x_, lower_], y_node] := Which[First[y] >= x, node[x, lower, y], First[y] <= lower, node[lower, y, x], True, node[y, lower, x]] add[node[x_, lower_, higher_], y_node] /; First[y] <= x := node[x, add[lower, y], higher] add[node[x_, lower_, higher_], y_node] := node[x, lower, add[higher, y]] SeedRandom[5]; nums = Null Timing[Do[nums = add[nums, Random[]], {5000}]; ] {0.968*Second, Null} None of the other stats changed. -----Original Message----- Count[nums, Null, Infinity] Count[nums, node[_, Null, _], Infinity] Count[nums, node[_, _, Null], Infinity] Count[nums, node[_, Null, Null], Infinity] Count[nums, node[___], Infinity] Count[nums, node[_, _node | _?NumericQ, _node | _?NumericQ], Infinity] Depth[nums] {4.406000000000001*Second, Null} 15001 240000 5001 (* Null values in the tree *) 2458 (* nodes without left offspring *) 2543 (* nodes without right offspring *) 1669 (* leaf nodes *) 4999 (* total nodes *) 1667 (* nodes with both left and right offspring *) 28 (* minus one, makes 27 levels in the tree *) It's clearly not ideal to store 5001 Nulls for 5000 actual values! Here's a method that doesn't change the algorithm much, but changes the storage method a great deal: First of all, I'm lazy, so I'll redefine <,>,<=,>=, etc. to make the code simpler: Unprotect[Less, Greater, LessEqual, GreaterEqual]; Less[a : _node ..] := Less @@ (First /@ {a}) LessEqual[a : _node ..] := LessEqual @@ (First /@ {a}) Greater[a : _node ..] := Greater @@ (First /@ {a}) GreaterEqual[a : _node ..] := GreaterEqual @@ (First /@ {a}) Protect[Less, Greater, LessEqual, GreaterEqual]; Again because I'm lazy, I'll define node so that I don't have to consciously avoid leaf nodes and unnecessary nesting: ClearAll[node] node[a___, node[b_], c___] = node[a, b, c]; node[node[a__]] = node[a]; Here's my add function: Clear[add] add[Null, x_] = node[x]; add[x_, (y_)?NumericQ] := add[node[x], node[y]] add[(x_)?NumericQ, y_] := add[node[x], node[y]] add[node[x_], node[y_]] := If[x > y, node[x, y], node[y, x]] (* <-- THERE! *) add[node[x_, lower_], y_node] := Which[y >= node[x], node[x, lower, y], y <= node[lower], node[lower, y, x], True, node[y, lower, x]] add[node[x_, lower_, higher_], y_node] /; node[y] <= node[x] := node[x, add[lower, y], higher] add[node[x_, lower_, higher_], y_node] := node[x, lower, add[higher, y]] Where I have a pointer is the code that prevents adding right-offspring to a leaf node. That avoids nodes like your node[a,Null,c]. When you would have node[a,b,Null], I have node[a,b]. When you would have node[a,Null,Null], I use a itself. Timing and space requirements are much better now: SeedRandom[5]; nums = Null Timing[Do[nums = add[nums, Random[]], {5000}]; ] {2.109*Second, Null} LeafCount[nums] ByteCount[nums] Count[nums, Null, Infinity] Count[nums, node[_], Infinity] Count[nums, node[_, _], Infinity] Count[nums, node[_, _, _], Infinity] Count[nums, node[___], Infinity] Depth[nums] 7852 (* 48% fewer leaf expressions, mostly due to Nulls and right-offspring eliminated *) 165624 (* 31% less storage *) 0 (* no Nulls, versus 5001 *) 0 (* no leaf nodes, versus 1669 *) (* no nodes without left-offspring, versus 2458 *) 705 (* nodes without right offspring, versus 2543 *) 2146 (* nodes with both left and right offspring, versus 1667 *) 2851 (* total nodes *) 21 (* 21 tree levels versus 27 *) The logical structure is identical to yours except that there are no nodes with only right-offspring. This incidentally tends to reduce tree depth. The storage format is far different: there are no trivial (childless) nodes, and nodes with left-offspring but no right-offspring are stored without a Null on the right. I could make this quite a bit faster, I think, if I spent time on the code to eliminate changes to Less, Greater, etc. and the frequent nesting followed by unnesting that goes on in the algorithm. It might take twice as much code, however, so I like it as is. You're probably better off storing these things in a heap, of course. Or -- better yet -- use the built-in Sort and Ordering functions. -----Original Message----- In[327]:= add[Null,x_Real]:=node[x,Null,Null] add[Null,node[x_Real,lower_,higher_]]=node[x,lower,higher] In[328]:= add[node[x_Real,lower_,higher_],y_Real]:= If[x>y,node[x,add[lower,y],higher],node[x,lower,add[higher,y]]] In[288]:= add[node[x_Real,lowerx_,higherx_],node[y_Real,lowery_,highery_]]:=If[x>y , node[x,add[lowerx,node[y,lowery,highery]],higherx], node[x,lowerx,add[higherx,node[y,lowery,highery]]] ] Now this is my attempt to test how fast my add works: SeedRandom[5]; Do[nums=add[nums,Random[]],{5000}];//Timing Out[333]= {13.279 Second,Null} RH7.3 running on an 1.4GHz Athlon with 1GB of ram). Questions: 1. Is this as fast as I can get my code to run? 2. Am I doing something obviously stupid? 3. would Compiling things help? Husain ==== Daniel, I don't mean to be overly critical of WRI's documentation -- it's very > good, as such things go. Nor do I like to overlook chances to make it > better. Do you? I'll try to address your specific remarks below and then return to this. > Your own comments below point out that we should expect vectors > resulting from NullSpace to be orthonormal by the usual > conjugate-symmetric inner product on C. But that's not spelled out > (documented) for dummies like myself who don't know that's a natural > result of singular values decomposition. My comments indicate that they will be orthonormal, for approximate numeric matrix input, so long as the implementation is based on SVD. But that's not by any means part of what NullSpace is required to do, and it is quite clearly implementation dependent. > The mission of NullSpace (or any function) is to adhere to > documentation, so reasonable persons may differ on whether orthogonality > is a feature we should depend on. Well... the name NullSpace gives a pretty good idea of what it needs to do. When you stop to consider all the matrix types (approximate numbers, exact numbers, symbolic, numeric but with non-number elements e.g. Pi and Sqrt[2]) you quickly realize that different types will most likely be handled by different methods, and these may have different characteristics of output e.g. normalized vectors, vectors with last nonzero component equal to one, or other quirks. But these are artifacts of implementation, not part of the requirements of a null space. > The mission of documentation is to tell us what to expect. When it > doesn't, the result is that we spend all this time discussing issues > online, trying to figure things out. A simple don't depend on > orthogonal results would be nice, if that's the intent. Perhaps. If there is strong reason to expect that people might depend on such results in the first place. I don't think there is. > In any case, I just spent ten minutes LOOKING for implementation notes > for NullSpace, and have not found any. Searching for implementation > notes doesn't help and there's no link from NullSpace. What use is > documentation I can't find? In general, I don't like Mathematica's quirky Help Browser, in which I > cannot search for anything that's not indexed. Every other help engine > on my computer (and there are hundreds) allows me to search for words, > and that's exactly what I need in order to find all mentions of > NullSpace. I found it in the back of the book, appendix A.9, Some Notes on Internal Implementation. See A.9.4, Approximate numerical linear algebra. In the Help Browser, click on The Mathematica Book and then type implementation in the window. The second URL it gives is A.9. I'll grant you one might reasonably desire a good search engine. I just want to indicate that what is there now is not entirely useless in this instance provided you at least know about the existance of these notes in the book. [...] To respond to your question at the beginning, yes, I also like to see improvements in documentation. I certainly agree that a good search engine would be useful in this case. Though note that I also may not be using the Help Browser to its fullest capabilities (I claim no great familiarity with the Mathematica user interface). As for the specific issue of what to document about NullSpace, suffice it to say that I do not favor making claims beyond what it is required to do. Daniel Lichtblau Wolfram Research ==== Daniel, I don't mean to be overly critical of WRI's documentation -- it's very good, as such things go. Nor do I like to overlook chances to make it better. Do you? Your own comments below point out that we should expect vectors resulting from NullSpace to be orthonormal by the usual conjugate-symmetric inner product on C. But that's not spelled out (documented) for dummies like myself who don't know that's a natural result of singular values decomposition. The mission of NullSpace (or any function) is to adhere to documentation, so reasonable persons may differ on whether orthogonality is a feature we should depend on. The mission of documentation is to tell us what to expect. When it doesn't, the result is that we spend all this time discussing issues online, trying to figure things out. A simple don't depend on orthogonal results would be nice, if that's the intent. In any case, I just spent ten minutes LOOKING for implementation notes for NullSpace, and have not found any. Searching for implementation notes doesn't help and there's no link from NullSpace. What use is documentation I can't find? In general, I don't like Mathematica's quirky Help Browser, in which I cannot search for anything that's not indexed. Every other help engine on my computer (and there are hundreds) allows me to search for words, and that's exactly what I need in order to find all mentions of NullSpace. -----Original Message----- > of Mathematica. Bobby The expected, and documented, behavior is that the output should be a basis for the null space, that is, solutions of the homogeneous matrix equation A.x==0. If this were to stop working then that would be a serious bug. Is this the behavior you mean? The implementation notes of the manual mention that approximate NullSpace is based on a singular values decomposition. This in fact gives resulting vectors that are orthonormal by the usual conjugate-symmetric inner product on C (though these are now not normal to the original vector in this same inner product, unless they are real-valued). But this basis-orthogonality is not part of the mission of NullSpace and moreover should not become part of it. Hence that particular (and implementation dependent) aspect of NullSpace should not become documented. Daniel Lichtblau Wolfram Research ==== I've this problem with NIntegrate: ********************************************** In[1]=f[y_,z_]: = (y^4/( (1+(8.44*10^-4)^2 * (1+z)^2 y^2) (Exp[z]+1) ) In[2]=F[z_]: = NIntegrate[f[y,z], {z, 0, Infinity} ********************************************** Mathematica 4.0 says: **************************************************************************** * NIntegrate: : inum : Integrand 1.07577/( 1+7.12336*10^-7 (1. + z)^2 ) is not numerical at {y}={1.}. **************************************************************************** ** What is it?? -- Rob_ Reply-To: kuska@informatik.uni-leipzig.de ==== do you realy think that 1.07577/( 1+7.12336*10^-7 (1. + z)^2 ) with the z inside is numerical ? or do you mean F[z_?NumericQ]: = NIntegrate[f[y,z], {y, 0, Infinity}] > I've this problem with NIntegrate: > ********************************************** > In[1]=f[y_,z_]: = (y^4/( (1+(8.44*10^-4)^2 * (1+z)^2 y^2) (Exp[z]+1) ) In[2]=F[z_]: = NIntegrate[f[y,z], {z, 0, Infinity} > ********************************************** Mathematica 4.0 says: **************************************************************************** > * > NIntegrate: : inum : Integrand 1.07577/( 1+7.12336*10^-7 (1. + z)^2 ) is not > numerical at {y}={1.}. > **************************************************************************** > ** What is it?? > -- > Rob_ ==== There is an error in the previous message. This is the just msg: ********************************************** In[1]=f[y_,z_]: = (y^4/( (1+(8.44*10^-4)^2 * (1+z)^2 y^2) (Exp[y]+1) ) In[2]=F[z_]: = NIntegrate[f[y,z], {y, 0, Infinity} ********************************************** Mathematica 4.0 says: **************************************************************************** NIntegrate: : inum : Integrand 1.07577/( 1+7.12336*10^-7 (1. + z)^2 ) is not numerical at {y}={1.}. **************************************************************************** What is it?? Moreover, I have used the following procedure: ****************************************************** If z<<1, f[y,z]=Sum[(-1)^n*a^n*(1+z)^2n *(y^(2n+4)/(Exp[y]+1),{n, 0, N} ] here a=(8.44*10^-4)^2 therefore: F[z]:=Sum[(-1)^n*a^n*(1+z)^2n *(y^(2n+4)/(Exp[y]+1) * Gamma[2 n+5]*Zeta[2 n+5],{n, 0, N} ] this series (N<+oo) only approximates the function for small z, and me it interests the behavior of F for z in the range 800-1300. Thx in advance. Rob_ ==== I've this problem with NIntegrate: > ********************************************** > In[1]=f[y_,z_]: = (y^4/( (1+(8.44*10^-4)^2 * (1+z)^2 y^2) (Exp[z]+1) ) In[2]=F[z_]: = NIntegrate[f[y,z], {z, 0, Infinity} > ********************************************** Mathematica 4.0 says: **************************************************************************** > * > NIntegrate: : inum : Integrand 1.07577/( 1+7.12336*10^-7 (1. + z)^2 ) is not > numerical at {y}={1.}. > **************************************************************************** > ** What is it?? > I doubt that you could input your equations as you stated: The first input has a missing ) The second one a missing ] In both equations there should be no space in := Apart from that the second definitions doesn't seem to make much sense, maybe it should be F[y_] := ... Have you assigned y some value before (I get a different error message, when I mix up the y and z as you did.) Good luck Alois -- Vienna University of Technology, A-1040 Wiedner Hauptstr. 8-10 ==== Daniel, We mostly agree. I simply look at things from a slightly different (user, non-WRI) perspective. My original point was -- in line with your own comments -- that people shouldn't depend on NullSpace to return orthogonal vectors even when it seems to do so, because it's not documented behavior. (Although YOU knew about it and could point to documentation of it! Sort of.) Unlike you, perhaps, I consider the lack of search capabilities in Help an EGREGIOUS failing -- a mistake no other Help engine makes (so far as I know). Fix that, and 90% of my complaints about documentation would go away; I could have found the implementation notes on NullSpace, for instance, without already KNOWING where they are. That is precisely what annoys me most of the time; the information is there, but I can't find it. (If searching for implementation works, why doesn't implementation notes yield anything at all?) Now that I've found this page, I'll add it to all the other information locations I'm trying to remember. (There's not much there anyway, though.) Bobby -----Original Message----- I'll try to address your specific remarks below and then return to this. > Your own comments below point out that we should expect vectors > resulting from NullSpace to be orthonormal by the usual > conjugate-symmetric inner product on C. But that's not spelled out > (documented) for dummies like myself who don't know that's a natural > result of singular values decomposition. My comments indicate that they will be orthonormal, for approximate numeric matrix input, so long as the implementation is based on SVD. But that's not by any means part of what NullSpace is required to do, and it is quite clearly implementation dependent. > The mission of NullSpace (or any function) is to adhere to > documentation, so reasonable persons may differ on whether orthogonality > is a feature we should depend on. Well... the name NullSpace gives a pretty good idea of what it needs to do. When you stop to consider all the matrix types (approximate numbers, exact numbers, symbolic, numeric but with non-number elements e.g. Pi and Sqrt[2]) you quickly realize that different types will most likely be handled by different methods, and these may have different characteristics of output e.g. normalized vectors, vectors with last nonzero component equal to one, or other quirks. But these are artifacts of implementation, not part of the requirements of a null space. > The mission of documentation is to tell us what to expect. When it > doesn't, the result is that we spend all this time discussing issues > online, trying to figure things out. A simple don't depend on > orthogonal results would be nice, if that's the intent. Perhaps. If there is strong reason to expect that people might depend on such results in the first place. I don't think there is. > In any case, I just spent ten minutes LOOKING for implementation notes > for NullSpace, and have not found any. Searching for implementation > notes doesn't help and there's no link from NullSpace. What use is > documentation I can't find? In general, I don't like Mathematica's quirky Help Browser, in which I > cannot search for anything that's not indexed. Every other help engine > on my computer (and there are hundreds) allows me to search for words, > and that's exactly what I need in order to find all mentions of > NullSpace. I found it in the back of the book, appendix A.9, Some Notes on Internal Implementation. See A.9.4, Approximate numerical linear algebra. In the Help Browser, click on The Mathematica Book and then type implementation in the window. The second URL it gives is A.9. I'll grant you one might reasonably desire a good search engine. I just want to indicate that what is there now is not entirely useless in this instance provided you at least know about the existance of these notes in the book. [...] To respond to your question at the beginning, yes, I also like to see improvements in documentation. I certainly agree that a good search engine would be useful in this case. Though note that I also may not be using the Help Browser to its fullest capabilities (I claim no great familiarity with the Mathematica user interface). As for the specific issue of what to document about NullSpace, suffice it to say that I do not favor making claims beyond what it is required to do. Daniel Lichtblau Wolfram Research ==== -----Original Message----- Yes, I can make a page break where I want one, but I can't eliminate > page-breaks that leave tons of white-space at the end of some pages. In > some of my notebooks, Mathematica puts the first cell on a page by > itself, using less than an inch --- and NOT because the next cell takes > a whole page. If I force the first two cells to stay together, the > problem just moves to the next page. It's enough to make Mathematica useless as a word processor all by > itself. -----Original Message----- > Kenny, Sympathy but no solution. I too have been trying to use Mathematica (v4.2 most recently) to type > maths papers and the like but I'm not ready to ditch LaTeX yet. There > are just too many cases where I cannot figure out how to achieve what I > want in Mathematica, things like: > - left brackets spanning multiple lines for defining hybrid functions; > - vertical alignment of equals signs in multi-line equations or > derivations; > - setting typefaces in tables of material. I figure most of this is do-able, but I don't have the time, or > patience, to spend too much time on it. So, I'll be the first cuser > when you write the guide to math publishing in Mathematica - I just hope > you 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 to find a complete guide to word processing with mathematica. >>Does anyone know where such a document could be found? > ==== Group, Sorry for the previous posting. The following solves (to some extent my problem) <1 0.02 times2=First[Timing[fLPartition[3,20];]]/. Second->1 0.18 times3=First[Timing[fLPartition[3,30];]]/. Second->1 0.701 times4=First[Timing[fLPartition[3,40];]]/. Second->1 1.913 times5=First[Timing[fLPartition[3,50];]]/. Second->1 4.236 times6=First[Timing[fLPartition[3,60];]]/. Second->1 8.222 times7=First[Timing[fLPartition[3,70];]]/. Second->1 14.39 times8=First[Timing[fLPartition[3,80];]]/. Second->1 23.563 So far so good. But for times9=First[Timing[fLPartition[3,90];]]/. Second->1 The computation does not always end. Some times it ends some times it does not (I aborted the process four times after 1 hour of computing and running out of, 2Gbites of, virtual memory) Emilio Martin-Serrano ___________________________________ Emilio Martin-Serrano Schlumberger Oil & Gas Business Consulting Principal 1325 South Dairy Ashford Road Houston, TX 77077 ==== I'm going to try using Mathematica to typeset some mathematics because I know that my students would benefit from being able to do the same., and they'll need someone who can answer their questions, etc. etc. I'm pretty comfortable with some things (e.g., Mathematica's recognition of latex commands for most characters), and now I'm ready to push the envelope a bit. Q1: How can I typeset a table or array of information (e.g., list of variables with definitions) so that information in one column is centered and information in another column is aligned flush left? I know exactly how I would do this in LaTeX, but I find no information on how to do this in Mathematica. Q2. How can I embed a Quicktime movie into a Mathematica notebook (e.g., so someone with MathReader can play the animation)? That's all for now. -- Jason Miller, Ph.D. Division of Mathematics and Computer Science Truman State University 100 East Normal St. Kirksville, MO 63501 http://vh216801.truman.edu 660.785.7430 ==== Jason, On item #1 I can help you. I have a notebook that explains the basics of creating tables using either TableForm or Gridboxes. Should get you started. It can be downloaded from my Mathematica course web site: http://www.higgins.ucdavis.edu/ECH198.php Use the link for lecture 3. Brian > I'm going to try using Mathematica to typeset some mathematics > because I know that my students would benefit from being able to do > the same., and they'll need someone who can answer their questions, > etc. etc. I'm pretty comfortable with some things (e.g., Mathematica's > recognition of latex commands for most characters), and now I'm ready > to push the envelope a bit. Q1: How can I typeset a table or array of information (e.g., list of > variables with definitions) so that information in one column is > centered and information in another column is aligned flush left? I > know exactly how I would do this in LaTeX, but I find no information > on how to do this in Mathematica. Q2. How can I embed a Quicktime movie into a Mathematica notebook > (e.g., so someone with MathReader can play the animation)? That's all for now. ==== It seems there is an error in the BinCounts function of Mathematica's standard <Left, FontSize->18], Cell[< Sept. 9, 2002 Mark D. Normand UMass Dept. of Food Science mnormand@foodsci.umass.edu >, Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[< Load the Data Manipulation functions from the Statistics Standard Package. >, Text, TextAlignment->Left], Cell[BoxData[ (<< Statistics`DataManipulation`)], Input] }, Open ]], Cell[CellGroupData[{ Cell[Assign data to the vector (1D list) dat., Text, TextAlignment->Left], Cell[BoxData[ ((dat = {63, 184, 23, 14, 17, 32, 26, 666, 27, 11, 28, 53, 37, 29, 4, 60, 7, 23, 94, 18, 43, 15, 74, 42, 81, 8, 7, 19, 0, 27, 87, 35, 3, 24, 94, 42, 4, 7, 18, 7, 38, 1, 0, 61, 3, 40, 22, 8, 3, 2, 0, 2, 5, 5, 7, 0, 9, 2, 11, 13, 105, 51, 36, 149, 147, 12, 1, 1, 7, 8, 10, 1, 7, 25, 38, 142, 15, 8, 9, 16, 0, 1, 4, 3, 31, 30, 10};))], Input] }, Open ]], Cell[CellGroupData[{ Cell[Count the number of dat data rows., Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ (datRows = Length[dat])], Input], Cell[BoxData[ (87)], Output] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[Find the minimum value in the dat data rows., Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ (datMin = Min[dat])], Input], Cell[BoxData[ (0)], Output] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[Find the maximum value in the dat data rows., Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ (datMax = Max[dat])], Input], Cell[BoxData[ (666)], Output] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[Calculate the number of histogram bins to use., Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ (datBins = Ceiling[Log[2, datRows]])], Input], Cell[BoxData[ (7)], Output] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[Set the width of a histogram bin., Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ (datBinWidth = N[((datMax - datMin))/datBins])], Input], Cell[BoxData[ (95.14285714285714`)], Output] }, Open ]] }, Open ]], Cell[CellGroupData[{ TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ Cell[BoxData[ ({76, 5, 0, 0, 0, 0, 0, 1})], Output] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[< For the above dat list, if datBins equals, for example, 5,7,10,14 of histogram bins generated) is one more than the assigned datBins value. Yet 13 or 16. Why!? (Is this a bug?) >, Text, TextAlignment->Left], Cell[CellGroupData[{ Cell[BoxData[ Cell[BoxData[ (8)], Output] }, Open ]] }, Open ]] }, ScreenRectangle->{{0, 1024}, {0, 748}}, CellGrouping->Manual, WindowSize->{520, 626}, WindowMargins->{{-1, Auatic}, {Auatic, -1}}, MacintoshSystemPageSetup->< 0080001804P000000_82@?okonh34`9B;@<5:0?l00;m009H0UP0000068dB`0B` 02d5X5k/02H80@4101P00BL?00400@00000000000000P0010000020D00000000 0000000000000000000004T400002004> ] ==== On my machine, replacing Null with made Husain's code 8 times faster on the same test: 0.592 seconds versus 4.688. I suppose making the same change in my code makes no difference because there are no Nulls in the tree (except one, very briefly), or perhaps because I have fewer patterns that involve Nulls. By rearranging the order of my rules, and deleting a couple that were never used, I made my algorithm about 15% faster, getting the time down to 0.906. Mihajlo's amazing discovery beats that, of course -- it's 50% faster. I suppose what we've found is that (1) my algorith makes the code four to five times faster by reducing expression complexity, tree depth, and memory use, but it's a slower algorithm otherwise (too many rules?); and (2) Mathematica pattern matching and code execution has a real problem with Null used as if it means something. Bobby -----Original Message----- add[node[x_Real, lower_, higher_], y_Real] := node[x, lower, add[higher, y]]; nums = NULL; SeedRandom[5]; Do[nums = add[nums, Random[]], {5000}] // Timing Out {5.11 Second, Null} Clear[add]; add[Null, x_Real] := node[x, Null, Null]; add[node[x_Real, lower_, higher_], y_Real] /; x > y := node[x, add[lower, y], higher]; add[node[x_Real, lower_, higher_], y_Real] := node[x, lower, add[higher, y]]; nums = Null; SeedRandom[5]; Do[nums = add[nums, Random[]], {5000}] // Timing Out {41.03 Second, Null} Also: Simplify[1==Null] 1==Null doesn't evaluate to false!! Btw. changing Null to NULL in your code doesn't infulence the speed... Mihajlo Vanevic mvane@EUnet.yu 2002-09-09 ************************************************************** * At 2002-09-09, 00:29:00 ************************************************************** >You're not doing anything dumb as far as I can see, but you're using far >more memory than necessary. Here are statistics for your algorithm on >my machine: SeedRandom[5]; >nums = Null >Do[nums = add[nums, Random[]], {5000}]; // Timing >nums // LeafCount >nums // ByteCount >Count[nums, Null, Infinity] >Count[nums, node[_, Null, _], Infinity] >Count[nums, node[_, _, Null], Infinity] >Count[nums, node[_, Null, Null], Infinity] >Count[nums, node[___], Infinity] >Count[nums, node[_, _node | _?NumericQ, _node | _?NumericQ], Infinity] >Depth[nums] {4.406000000000001*Second, Null} >15001 >240000 >5001 (* Null values in the tree *) >2458 (* nodes without left offspring *) >2543 (* nodes without right offspring *) >1669 (* leaf nodes *) >4999 (* total nodes *) >1667 (* nodes with both left and right offspring *) >28 (* minus one, makes 27 levels in the tree *) It's clearly not ideal to store 5001 Nulls for 5000 actual values! Here's a method that doesn't change the algorithm much, but changes the >storage method a great deal: First of all, I'm lazy, so I'll redefine <,>,<=,>=, etc. to make the >code simpler: Unprotect[Less, Greater, LessEqual, GreaterEqual]; >Less[a : _node ..] := Less @@ (First /@ {a}) >LessEqual[a : _node ..] := LessEqual @@ (First /@ {a}) >Greater[a : _node ..] := Greater @@ (First /@ {a}) >GreaterEqual[a : _node ..] := GreaterEqual @@ (First /@ {a}) >Protect[Less, Greater, LessEqual, GreaterEqual]; Again because I'm lazy, I'll define node so that I don't have to >consciously avoid leaf nodes and unnecessary nesting: ClearAll[node] >node[a___, node[b_], c___] = node[a, b, c]; >node[node[a__]] = node[a]; Here's my add function: Clear[add] >add[Null, x_] = node[x]; >add[x_, (y_)?NumericQ] := add[node[x], node[y]] >add[(x_)?NumericQ, y_] := add[node[x], node[y]] >add[node[x_], node[y_]] := If[x > y, node[x, y], node[y, x]] (* <-- >THERE! *) >add[node[x_, lower_], y_node] := > Which[y >= node[x], node[x, lower, y], > y <= node[lower], node[lower, y, x], True, node[y, lower, x]] >add[node[x_, lower_, higher_], y_node] /; > node[y] <= node[x] := node[x, add[lower, y], higher] >add[node[x_, lower_, higher_], y_node] := node[x, lower, add[higher, >y]] Where I have a pointer is the code that prevents adding right-offspring >to a leaf node. That avoids nodes like your node[a,Null,c]. When you >would have node[a,b,Null], I have node[a,b]. When you would have >node[a,Null,Null], I use a itself. Timing and space requirements are much better now: SeedRandom[5]; >nums = Null >Timing[Do[nums = add[nums, Random[]], {5000}]; ] {2.109*Second, Null} LeafCount[nums] >ByteCount[nums] >Count[nums, Null, Infinity] >Count[nums, node[_], Infinity] >Count[nums, node[_, _], Infinity] >Count[nums, node[_, _, _], Infinity] >Count[nums, node[___], Infinity] >Depth[nums] 7852 (* 48 fewer leaf expressions, > mostly due to Nulls and right-offspring eliminated *) >165624 (* 31 less storage *) >0 (* no Nulls, versus 5001 *) >0 (* no leaf nodes, versus 1669 *) > (* no nodes without left-offspring, versus 2458 *) >705 (* nodes without right offspring, versus 2543 *) >2146 (* nodes with both left and right offspring, versus 1667 *) >2851 (* total nodes *) >21 (* 21 tree levels versus 27 *) The logical structure is identical to yours except that there are no >nodes with only right-offspring. This incidentally tends to reduce tree >depth. The storage format is far different: there are no trivial (childless) >nodes, and nodes with left-offspring but no right-offspring are stored >without a Null on the right. I could make this quite a bit faster, I think, if I spent time on the >code to eliminate changes to Less, Greater, etc. and the frequent >nesting followed by unnesting that goes on in the algorithm. It might >take twice as much code, however, so I like it as is. You're probably better off storing these things in a heap, of course. Or -- better yet -- use the built-in Sort and Ordering functions. -----Original Message----- So Slow? I am trying to implement a very simple sorted tree to quickly store some >real numbers I need. I have written an add, delete, minimum, and pop >(delete the lowest value) function and they seem to work ok but are very >slow. Let's just look @ my implementation of the add part: >nums=Null;(*my initial blank Tree) >In[326]:= >Clear[add] In[327]:= >add[Null,x_Real]:=node[x,Null,Null] add[Null,node[x_Real,lower_,higher_]]=node[x,lower,higher] In[328]:= >add[node[x_Real,lower_,higher_],y_Real]:= > If[x>y,node[x,add[lower,y],higher],node[x,lower,add[higher,y]]] In[288]:= >add[node[x_Real,lowerx_,higherx_],node[y_Real,lowery_,highery_]]:=If[x> y >, > node[x,add[lowerx,node[y,lowery,highery]],higherx], > node[x,lowerx,add[higherx,node[y,lowery,highery]]] > ] Now this is my attempt to test how fast my add works: SeedRandom[5]; >Do[nums=add[nums,Random[]],{5000}];//Timing Out[333]= >{13.279 Second,Null} RH7.3 >running on an 1.4GHz Athlon with 1GB of ram). Questions: >1. Is this as fast as I can get my code to run? >2. Am I doing something obviously stupid? >3. would Compiling things help? >Husain ************************************************************** ==== The following is a button that will export a selected graphic as a JPEG, ButtonBox[JPEG, ButtonFunction :> Module[{sel}, SelectionMove[InputNotebook[], All, Cell]; sel = NotebookRead[InputNotebook[]]; If[Head[sel] === Cell, If[!ValueQ[dpi], dpi = Auatic]; ImageSize -> Auatic, ImageResolution -> dpi, ConversionOptions -> {Quality -> 75}]]], ButtonEvaluator -> Auatic, Active -> True]//DisplayForm It works fine as long as the graphic is no wider than about 330 pixels. If the graphic is wider than 330 pixels, then the size of exported graphic is correct, but everything on the right beyond the first 330 pixels is blank. There seems to be no problem with the height though. Any ideas on how to fix or get around this problem? --- Reply-To: kuska@informatik.uni-leipzig.de ==== not with Mathematica 4.2. But the problem can be the page width that may be used to convert the Cell[], if you can you should use the Graphics[]/Graphics3D[] .. inside the Cell or the PostScript string in the GraphicsData[] The following is a button that will export a selected graphic as a JPEG, ButtonBox[JPEG, > ButtonFunction : Module[{sel}, > SelectionMove[InputNotebook[], All, Cell]; > sel = NotebookRead[InputNotebook[]]; > If[Head[sel] === Cell, > If[!ValueQ[dpi], dpi = Auatic]; > ImageSize -> Auatic, ImageResolution -> dpi, > ConversionOptions -> {Quality -> 75}]]], > ButtonEvaluator -> Auatic, Active -> True]//DisplayForm It works fine as long as the graphic is no wider than about 330 pixels. > If the graphic is wider than 330 pixels, then the size of exported > graphic is correct, but everything on the right beyond the first 330 > pixels is blank. There seems to be no problem with the height though. Any ideas on how to fix or get around this problem? --- > ==== I am having some trouble with numerical integration. !(((Timing[ ansOUT = {}; [IndentingNewLine]For[EbNo = 0, EbNo <= 2, EbNo += 10/10, [IndentingNewLine]linEbNo = 10^(EbNo/10); [IndentingNewLine]M = 64; [IndentingNewLine]K = Log[2, M]; [IndentingNewLine]p = (M/2)/(M - 1) Sum[ Binomial[M - 1, n] (((-1)))^(n + 1)/(1 + n) Exp[(- linEbNo) K n/(1 + n)], {n, 1, M - 1}]; [IndentingNewLine]Pm = (1/@(2 [Pi])) NIntegrate[ 1 - ((((1/@(2 [Pi])) NIntegrate[ Exp[(-x^2)/2], {x, (-[Infinity]), y}, WorkingPrecision -> 32 ]))^(M - 1)) Exp[((-1)/2) ((y - @(2 K linEbNo)))^2], {y, (- [Infinity]), [Infinity]}, WorkingPrecision -> 32]; [IndentingNewLine]Pb = ((M/2)/(M - 1)) Pm; [IndentingNewLine]AppendTo[ ansOUT, {N[a, 4], N[p, 8], N[Pb, 8]}];[IndentingNewLine]];])([IndentingNewLine]))) TableForm[ansOUT] The code can also be downloaded from http://www.itee.uq.edu.au/~dsalman/_downloads/Proakis.nb Col. 2(symblic formula) and 3 (result of numerical integration) of the variable ansOUT should be identical. However the numerical integration is giving a 0 answer. Can anyone please suggest how a solution to this problem. Salman ==== Reply-To: kuska@informatik.uni-leipzig.de ==== Needs[Graphics`Graphics`] LogPlot[Exp[-x],{x,0,10}] ?? ==== <>ac.dat all done in an instant (relatively speaking) and without the transient burst of kernel memory usage. Ditto Import and Get. Can anyone explain this? Mike ==== Dear DrBob, I fear that your diagnosis of the efficiency problems in Husain's binary tree implementation is not correct. Furthermore your code is inefficient and awkward. Consider the following implementation using list to represent the tree. It is three lines long and is three times faster than yours on my computer though it uses twice as much memory as your representation. myadd[{},x_]:={x,{},{}}; myadd[{v_,a_,b_},x_]:=If[v>x, {v,myadd[a,x],b}, {v,a,myadd[b,x]}] t={}; Timing[Do[t=myadd[t,Random[]],{5000}];] {1.98 Second,Null} The real problem with the implementation is probably that the tree structure has to be 'copied' each time when add is called. You have to write tree = add[tree, x] t2 = add[t1, x] then you have *two* trees t1 and t2 which differ by the one element inserted by add. To the best of my knowledge, the internal representation of t1 and t2 in Mathematica is a tree too. But in principle the user code is duplicated in Mathematica. A further problem of all direct implementations is that they cannot be compiled (since they use pattern matching and the tree is a nested structure). Thus I guess that the efficient way to implement tree data structures is the indirect one using heaps (i.e. a combination of contiguous arrays with length known in advance and hashing). Johannes Ludsteck > You're not doing anything dumb as far as I can see, but you're using far > more memory than necessary. Here are statistics for your algorithm on > my machine: SeedRandom[5]; > nums = Null > Do[nums = add[nums, Random[]], {5000}]; // Timing > nums // LeafCount > nums // ByteCount > Count[nums, Null, Infinity] > Count[nums, node[_, Null, _], Infinity] > Count[nums, node[_, _, Null], Infinity] > Count[nums, node[_, Null, Null], Infinity] > Count[nums, node[___], Infinity] > Count[nums, node[_, _node | _?NumericQ, _node | _?NumericQ], Infinity] > Depth[nums] {4.406000000000001*Second, Null} > 15001 > 240000 > 5001 (* Null values in the tree *) > 2458 (* nodes without left offspring *) > 2543 (* nodes without right offspring *) > 1669 (* leaf nodes *) > 4999 (* total nodes *) > 1667 (* nodes with both left and right offspring *) > 28 (* minus one, makes 27 levels in the tree *) It's clearly not ideal to store 5001 Nulls for 5000 actual values! Here's a method that doesn't change the algorithm much, but changes the > storage method a great deal: First of all, I'm lazy, so I'll redefine <,>,<=,>=, etc. to make the > code simpler: Unprotect[Less, Greater, LessEqual, GreaterEqual]; > Less[a : _node ..] := Less @@ (First /@ {a}) > LessEqual[a : _node ..] := LessEqual @@ (First /@ {a}) > Greater[a : _node ..] := Greater @@ (First /@ {a}) > GreaterEqual[a : _node ..] := GreaterEqual @@ (First /@ {a}) > Protect[Less, Greater, LessEqual, GreaterEqual]; Again because I'm lazy, I'll define node so that I don't have to > consciously avoid leaf nodes and unnecessary nesting: ClearAll[node] > node[a___, node[b_], c___] = node[a, b, c]; > node[node[a__]] = node[a]; Here's my add function: Clear[add] > add[Null, x_] = node[x]; > add[x_, (y_)?NumericQ] := add[node[x], node[y]] > add[(x_)?NumericQ, y_] := add[node[x], node[y]] > add[node[x_], node[y_]] := If[x > y, node[x, y], node[y, x]] (* <-- > THERE! *) > add[node[x_, lower_], y_node] := > Which[y >= node[x], node[x, lower, y], > y <= node[lower], node[lower, y, x], True, node[y, lower, x]] > add[node[x_, lower_, higher_], y_node] /; > node[y] <= node[x] := node[x, add[lower, y], higher] > add[node[x_, lower_, higher_], y_node] := node[x, lower, add[higher, > y]] Where I have a pointer is the code that prevents adding right-offspring > to a leaf node. That avoids nodes like your node[a,Null,c]. When you > would have node[a,b,Null], I have node[a,b]. When you would have > node[a,Null,Null], I use a itself. Timing and space requirements are much better now: SeedRandom[5]; > nums = Null > Timing[Do[nums = add[nums, Random[]], {5000}]; ] {2.109*Second, Null} LeafCount[nums] > ByteCount[nums] > Count[nums, Null, Infinity] > Count[nums, node[_], Infinity] > Count[nums, node[_, _], Infinity] > Count[nums, node[_, _, _], Infinity] > Count[nums, node[___], Infinity] > Depth[nums] 7852 (* 48% fewer leaf expressions, > mostly due to Nulls and right-offspring eliminated *) > 165624 (* 31% less storage *) > 0 (* no Nulls, versus 5001 *) > 0 (* no leaf nodes, versus 1669 *) > (* no nodes without left-offspring, versus 2458 *) > 705 (* nodes without right offspring, versus 2543 *) > 2146 (* nodes with both left and right offspring, versus 1667 *) > 2851 (* total nodes *) > 21 (* 21 tree levels versus 27 *) The logical structure is identical to yours except that there are no > nodes with only right-offspring. This incidentally tends to reduce tree > depth. The storage format is far different: there are no trivial (childless) > nodes, and nodes with left-offspring but no right-offspring are stored > without a Null on the right. I could make this quite a bit faster, I think, if I spent time on the > code to eliminate changes to Less, Greater, etc. and the frequent > nesting followed by unnesting that goes on in the algorithm. It might > take twice as much code, however, so I like it as is. You're probably better off storing these things in a heap, of course. Or -- better yet -- use the built-in Sort and Ordering functions. -----Original Message----- I am trying to implement a very simple sorted tree to quickly store some > real numbers I need. I have written an add, delete, minimum, and pop > (delete the lowest value) function and they seem to work ok but are very > slow. Let's just look @ my implementation of the add part: > nums=Null;(*my initial blank Tree) > In[326]:= > Clear[add] In[327]:= > add[Null,x_Real]:=node[x,Null,Null] add[Null,node[x_Real,lower_,higher_]]=node[x,lower,higher] In[328]:= > add[node[x_Real,lower_,higher_],y_Real]:= > If[x>y,node[x,add[lower,y],higher],node[x,lower,add[higher,y]]] In[288]:= > add[node[x_Real,lowerx_,higherx_],node[y_Real,lowery_,highery_]]:=If[x>y > , > node[x,add[lowerx,node[y,lowery,highery]],higherx], > node[x,lowerx,add[higherx,node[y,lowery,highery]]] > ] Now this is my attempt to test how fast my add works: SeedRandom[5]; > Do[nums=add[nums,Random[]],{5000}];//Timing Out[333]= > {13.279 Second,Null} RH7.3 > running on an 1.4GHz Athlon with 1GB of ram). Questions: > 1. Is this as fast as I can get my code to run? > 2. Am I doing something obviously stupid? > 3. would Compiling things help? > Husain > <><><><><><><><><><><><> Johannes Ludsteck Economics Department University of Regensburg Universitaetsstrasse 31 93053 Regensburg ==== Is it possible to configure a style sheet so that functions or data shown in graphs are coloured (colored :)) in working but black in printing...without actually re-executing the function that created the graphics? thanks Mike Reply-To: kuska@informatik.uni-leipzig.de ==== no because the PostScript creation know nothing about the notebook frontend. Mathematica does its graphics in PostScript since version 1.0 but the notebook frontend was introduced in version 2.x and the actual typsetting/markup format was introduced in version 3.0. Is it possible to configure a style sheet so that functions or data shown in > graphs are coloured (colored :)) in working but black in printing...without > actually re-executing the function that created the graphics? ==== A commonly used symbol for the Floor function is a square bracket with the upper indents removed. Is this symbol part of Mathematica's built-in symbols? I think not, so then, can it be constructed by hand. ==== It turns out that Mathematica does have the floor symbols! Use esc l f esc for the left floor and esc r f esc for the right floor. In fact, if you type in esc l f esc a esc r f esc Mathematica will interpret the input as Floor[a]. Carl Woll Physics Dept U of Washington A commonly used symbol for the Floor function is a square bracket with the > upper indents removed. Is this symbol part of Mathematica's built-in symbols? I > think not, so then, can it be constructed by hand. ==== > A commonly used symbol for the Floor function is a square bracket with > the upper indents removed. Is this symbol part of Mathematica's built-in > symbols? Sure. Just look for Floor in A.10 Listing of Major Built-In Mathematica Objects! For example, you can use [LeftFloor] and [RightFloor] to get the symbols you want. ==== what may [LeftFloor] and [RightFloor] show on screen ? It's easy to find in the complete characters palette | Operators | General > A commonly used symbol for the Floor function is a square bracket with the > upper indents removed. Is this symbol part of Mathematica's built-in symbols? I > think not, so then, can it be constructed by hand. ==== > A commonly used symbol for the Floor function is a square bracket with the > upper indents removed. Is this symbol part of Mathematica's built-in symbols? I > think not, so then, can it be constructed by hand. [LeftFloor] and [RightFloor] can be used for this. Look at Chapter Advanced Mathematics/Mathematical and Other Notations/ Operators. -- / Alexander Dreyer, Dipl.-Math. - Abteilung Adaptive Systeme / Fraunhofer Institut fuer Techno- und Wirtschaftsmathematik (ITWM) Gottlieb-Daimler-Strasse, Geb. 7^2=49/313 D-67663 Kaiserslautern / ==== Needs[Graphics`Graphics`] LogLinearPlot[Sqrt[x], {x, 1., 1000}, ImageSize -> 500]; ==== > I'm pretty comfortable with some things (e.g., Mathematica's > recognition of latex commands for most characters), and now I'm ready > to push the envelope a bit. Q1: How can I typeset a table or array of information (e.g., list of > variables with definitions) so that information in one column is > centered and information in another column is aligned flush left? I > know exactly how I would do this in LaTeX, but I find no information > on how to do this in Mathematica. See Section 2.8.11 of _The Mathematica Book_. http://documents.wolfram.com/v4/MainBook/2.8.11.html At the beginning of that section, you will find a description of the options that one may apply to GridBox[] objects, which are used to construct tables, matrices, and button palettes. More information on the values that these options may take can be found in the front end option documentation. http://documents.wolfram.com/v4/OtherInformation/GridBoxOptions.html With the current selection in your notebook being the table itself, you can use the Option Inspector dialog (with its scope set to selection) to change the ColumnAlignments option to {Center, Left}, which means that first column elements should be centered, and that subsequent columns should be left justified. -- ==== I don't know about elegance or simplicity, but I do have a FASTER solution, incidentally using combinations code that often comes in handy for me. Here's that code: << DiscreteMath`Combinatorica`; ClearAll[combinations]; r = Range[1, 9]; combinations::usage = combinations[list,n:{__Integer}] lists the combinations of list taken n at a time; combinations[r_List, n_Integer, {}] := If[n > Length@r, {}, DiscreteMath`Combinatorica`KSubsets[r, n]]; combinations[r_List, n_Integer, e_?VectorQ] := Join[e, #] & /@ DiscreteMath`Combinatorica`KSubsets[Complement[r, e], n]; combinations[r_List, n_Integer, e : {__?VectorQ}] := Flatten[combinations[r, n, #] & /@ e, 1]; combinations[r_List, n : {__Integer}] := Which[Plus @@ n == Length@r, Join[#, Complement[r, #]] & /@ combinations[r, Drop[n, -1]], Plus @@ n > Length@r, {}, True, Fold[combinations[r, #2, #1] &, {}, n]] The following duplicates your fLPartitions function, but is about 40% faster and uses about 15% less memory (for this example): Quit << DiscreteMath`Combinatorica` fLPartition[ksN_, r_] := With[{ks = KSubsets[ Range[r], ksN], ks1 = Partition[Range[r], 1]}, MapThread[ Complement[Append[#1, #2], (Sequence @@ {#1} & ) /@ Partition[#2, 1]] & , {Array[ks1 & , Length[ ks]], ks}]] Module[{ksN = 2}, While[ksN < r + 1, s = fLPartition[ksN, r]; ksN = ksN + 1; Print[s]]] Timing[fLPartition[3, 80]; ] {11.389999999999999*Second, Null} MaxMemoryUsed[] 233956264 Quit wantCombinations (* this loads my cominations function definition *) ClearAll[mine] mine[k_, n_] := (Join[List /@ Complement[ Range[n], #1], {#1}] & ) /@ combinations[Range[n], {k}] Timing[mine[3, 80]; ] {6.734999999999999*Second, Null} MaxMemoryUsed[] 195828264 However, I would change the storage structure to replace {{1},{2},...,{7}} with just {1,2,...7}. That cuts memory usage by 85% from your method for this example: Quit wantCombinations ClearAll[mySecond] mySecond[k_, n_] := (Join[Complement[Range[n], #1], {#1}] & ) /@ combinations[Range[n], {k}] Timing[mySecond[3, 80]; ] {2.922*Second, Null} (* versus 11.485 with your code, 6.7 with my other code *) MaxMemoryUsed[] 43996408 (* 44MB versus 234MB with your code, 196 MB with my other code *) Those comparisons include memory used during code execution; byte counts for that example are 134,742,424 for my storage method and 287,560,024 for yours. If you sometimes have to have a partition in your format, its easy enough to temporarily convert it from mine with the following function: f = Join[List /@ Drop[#1,-1], {Last[#1]}] & ; First[mySecond[2, 10]] f[%] {3, 4, 5, 6, 7, 8, 9, 10, {1, 2}} {{3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2}} First[fLPartition[2, 10]] {{3}, {4}, {5}, {6}, {7}, {8}, {9}, {10}, {1, 2}} -----Original Message----- [ksN<(r+1),s=fLPartition[ksN,r];ksN=ksN+1;Print[s]]] However, I do not like this solution since to my taste is to complex and rather inelegant. Nest or Fold would do better. Any Idea? Just an additional remark. Besides, I do not need the whole taxonomy, but a subset satisfying a minimizing condition on connectivity implied by some relationships among the elements. This condition is satisfied by one and only one element, or class of equivalence in each level of the hierarchy. So it is that very class which matter and consequently the procedure should give control on each class as the total structure is being generated. Once it is known which class satisfies the condition it is possible to avoid generating the whole structure and the subsequent combinatorial explosion. With just This an addition to the previous posting In fact, the thing is even worst. These are the CPU times for the following computations times1=First[Timing[fLPartition[3,10];]]/. Second->1 0.02 times2=First[Timing[fLPartition[3,20];]]/. Second->1 0.18 times3=First[Timing[fLPartition[3,30];]]/. Second->1 0.701 times4=First[Timing[fLPartition[3,40];]]/. Second->1 1.913 times5=First[Timing[fLPartition[3,50];]]/. Second->1 4.236 times6=First[Timing[fLPartition[3,60];]]/. Second->1 8.222 times7=First[Timing[fLPartition[3,70];]]/. Second->1 14.39 times8=First[Timing[fLPartition[3,80];]]/. Second->1 23.563 So far so good. But for times9=First[Timing[fLPartition[3,90];]]/. Second->1 The computation does not always end. Some times it ends some times it does not (I aborted the process four times after 1 hour of computing and running out of, 2Gbites of, virtual memory) ==== I wasn't trying to write a truly efficient algorithm -- I consciously tried not to change the algorithm (or even the storage method) a lot, so that Husain could move along from where he was. Of course a heap would be more efficient. In fact, given the built-in Sort and Ordering functions, I see no reason to build sorted binary trees in Mathematica at all (as I said in my first post) -- except as a step in the process of learning. Your solution (on my machine) is equivalent in speed to Mihajlo's simple modification of Husain's original -- just by replacing Null with anything else -- and using instead of {} in your code makes no difference. The following code is Husain's, but uses List rather than node and {} rather than Null, and builds precisely the same tree as your code. Yet it's 10% slower than your code. I think that's simply because there are more rules -- the same reason my code is slower, because I have even more rules, so that more time is spent on pattern-matching. Clear[add] add[{}, x_Real] := {x, {}, {}} add[{}, {x_Real, lower_, higher_}] = {x, lower, higher}; add[{x_Real, lower_, higher_}, y_Real] := If[x > y, {x, add[ lower, y], higher}, {x, lower, add[higher, y]}] add[{x_Real, lowerx_, higherx_}, { y_Real, lowery_, highery_}] := If[x > y, {x, add[lowerx, {y, lowery, highery}], higherx}, {x, lowerx, add[higherx, {y, lowery, highery}]}] In fact, using List instead of node slowed the algorithm compared to Mihajlo's solution, by the same 10%. I think that's because it requires more pattern matching, to distinguish {} from other List forms in the tree at each call. As for the inefficiency inherent in copying trees, how would you avoid that, other than using something like a heap? Even with a heap, it has to be a global variable (not a function argument) to avoid copying. Bobby -----Original Message----- than yours on my computer though it uses twice as much memory as your representation. myadd[{},x_]:={x,{},{}}; myadd[{v_,a_,b_},x_]:=If[v>x, {v,myadd[a,x],b}, {v,a,myadd[b,x]}] t={}; Timing[Do[t=myadd[t,Random[]],{5000}];] {1.98 Second,Null} The real problem with the implementation is probably that the tree structure has to be 'copied' each time when add is called. You have to write tree = add[tree, x] t2 = add[t1, x] then you have *two* trees t1 and t2 which differ by the one element inserted by add. To the best of my knowledge, the internal representation of t1 and t2 in Mathematica is a tree too. But in principle the user code is duplicated in Mathematica. A further problem of all direct implementations is that they cannot be compiled (since they use pattern matching and the tree is a nested structure). Thus I guess that the efficient way to implement tree data structures is the indirect one using heaps (i.e. a combination of contiguous arrays with length known in advance and hashing). Johannes Ludsteck > You're not doing anything dumb as far as I can see, but you're using far > more memory than necessary. Here are statistics for your algorithm on > my machine: SeedRandom[5]; > nums = Null > Do[nums = add[nums, Random[]], {5000}]; // Timing > nums // LeafCount > nums // ByteCount > Count[nums, Null, Infinity] > Count[nums, node[_, Null, _], Infinity] > Count[nums, node[_, _, Null], Infinity] > Count[nums, node[_, Null, Null], Infinity] > Count[nums, node[___], Infinity] > Count[nums, node[_, _node | _?NumericQ, _node | _?NumericQ], Infinity] > Depth[nums] {4.406000000000001*Second, Null} > 15001 > 240000 > 5001 (* Null values in the tree *) > 2458 (* nodes without left offspring *) > 2543 (* nodes without right offspring *) > 1669 (* leaf nodes *) > 4999 (* total nodes *) > 1667 (* nodes with both left and right offspring *) > 28 (* minus one, makes 27 levels in the tree *) It's clearly not ideal to store 5001 Nulls for 5000 actual values! Here's a method that doesn't change the algorithm much, but changes the > storage method a great deal: First of all, I'm lazy, so I'll redefine <,>,<=,>=, etc. to make the > code simpler: Unprotect[Less, Greater, LessEqual, GreaterEqual]; > Less[a : _node ..] := Less @@ (First /@ {a}) > LessEqual[a : _node ..] := LessEqual @@ (First /@ {a}) > Greater[a : _node ..] := Greater @@ (First /@ {a}) > GreaterEqual[a : _node ..] := GreaterEqual @@ (First /@ {a}) > Protect[Less, Greater, LessEqual, GreaterEqual]; Again because I'm lazy, I'll define node so that I don't have to > consciously avoid leaf nodes and unnecessary nesting: ClearAll[node] > node[a___, node[b_], c___] = node[a, b, c]; > node[node[a__]] = node[a]; Here's my add function: Clear[add] > add[Null, x_] = node[x]; > add[x_, (y_)?NumericQ] := add[node[x], node[y]] > add[(x_)?NumericQ, y_] := add[node[x], node[y]] > add[node[x_], node[y_]] := If[x > y, node[x, y], node[y, x]] (* <-- > THERE! *) > add[node[x_, lower_], y_node] := > Which[y >= node[x], node[x, lower, y], > y <= node[lower], node[lower, y, x], True, node[y, lower, x]] > add[node[x_, lower_, higher_], y_node] /; > node[y] <= node[x] := node[x, add[lower, y], higher] > add[node[x_, lower_, higher_], y_node] := node[x, lower, add[higher, > y]] Where I have a pointer is the code that prevents adding right-offspring > to a leaf node. That avoids nodes like your node[a,Null,c]. When you > would have node[a,b,Null], I have node[a,b]. When you would have > node[a,Null,Null], I use a itself. Timing and space requirements are much better now: SeedRandom[5]; > nums = Null > Timing[Do[nums = add[nums, Random[]], {5000}]; ] {2.109*Second, Null} LeafCount[nums] > ByteCount[nums] > Count[nums, Null, Infinity] > Count[nums, node[_], Infinity] > Count[nums, node[_, _], Infinity] > Count[nums, node[_, _, _], Infinity] > Count[nums, node[___], Infinity] > Depth[nums] 7852 (* 48% fewer leaf expressions, > mostly due to Nulls and right-offspring eliminated *) > 165624 (* 31% less storage *) > 0 (* no Nulls, versus 5001 *) > 0 (* no leaf nodes, versus 1669 *) > (* no nodes without left-offspring, versus 2458 *) > 705 (* nodes without right offspring, versus 2543 *) > 2146 (* nodes with both left and right offspring, versus 1667 *) > 2851 (* total nodes *) > 21 (* 21 tree levels versus 27 *) The logical structure is identical to yours except that there are no > nodes with only right-offspring. This incidentally tends to reduce tree > depth. The storage format is far different: there are no trivial (childless) > nodes, and nodes with left-offspring but no right-offspring are stored > without a Null on the right. I could make this quite a bit faster, I think, if I spent time on the > code to eliminate changes to Less, Greater, etc. and the frequent > nesting followed by unnesting that goes on in the algorithm. It might > take twice as much code, however, so I like it as is. You're probably better off storing these things in a heap, of course. Or -- better yet -- use the built-in Sort and Ordering functions. -----Original Message----- So Slow? I am trying to implement a very simple sorted tree to quickly store some > real numbers I need. I have written an add, delete, minimum, and pop > (delete the lowest value) function and they seem to work ok but are very > slow. Let's just look @ my implementation of the add part: > nums=Null;(*my initial blank Tree) > In[326]:= > Clear[add] In[327]:= > add[Null,x_Real]:=node[x,Null,Null] add[Null,node[x_Real,lower_,higher_]]=node[x,lower,higher] In[328]:= > add[node[x_Real,lower_,higher_],y_Real]:= > If[x>y,node[x,add[lower,y],higher],node[x,lower,add[higher,y]]] In[288]:= > add[node[x_Real,lowerx_,higherx_],node[y_Real,lowery_,highery_]]:=If[x>y > , > node[x,add[lowerx,node[y,lowery,highery]],higherx], > node[x,lowerx,add[higherx,node[y,lowery,highery]]] > ] Now this is my attempt to test how fast my add works: SeedRandom[5]; > Do[nums=add[nums,Random[]],{5000}];//Timing Out[333]= > {13.279 Second,Null} RH7.3 > running on an 1.4GHz Athlon with 1GB of ram). Questions: > 1. Is this as fast as I can get my code to run? > 2. Am I doing something obviously stupid? > 3. would Compiling things help? ==== I am interested in creating a slide presentation using slide view mode in Mathematica 4.2. I would like to have hyperlinks between different slides ( e.g., a hyperlink in say slide 6 takes me back to slide 3). I have no difficulty in creating such a hyperlink (using tags) when I am in the author mode. But the hyperlink does not work when a revert back to slide view mode. I guess it has somthing to do with the fact that in author view mode your slides are a subset of a single notebook but in slide view mode each slide becomes a distinct notebook. So my question: Is there a way to reference the individual slides using the hyperlink command? Brian ==== This graphs f from 0 to 8 Pi on a logarithmic horizontal scale: f[x_] := 1 + Abs[Sin[x]/x] Plot[f@Exp@x, {x, 0, Log[8Pi]}, PlotRange -> All]; If you want both scales to be logarithmic, use Plot[Log@f@Exp@x, {x, 0, Log[8Pi]}, PlotRange -> All]; The first is always possible if you only have positive x to deal with, and the second is possible if x and y are BOTH positive. ==== > This graphs f from 0 to 8 Pi on a logarithmic horizontal scale: f[x_] := 1 + Abs[Sin[x]/x] > Plot[f@Exp@x, {x, 0, Log[8Pi]}, PlotRange -> All]; The graph itself is indeed what is desired, but the labelling of the x-axis is then incorrect for a logarithmic scale. For example, regardless of the scale used, any proper graph must show that the first minimum occurs at x = Pi. But the method you suggest makes it seem that it occurs at approx. 1.14 (precisely Log[Pi]). > -----Original Message----- > ==== To find n digits before the decimal point of a number, we can proceed in the > following way. We compute the number in sufficiently many digits, then take > the Floor of the result (i.e. we round it downwards to an integer) and > finally take the result modulo 10^n. Mod[Floor[ N[(Sqrt[2] + Sqrt[3])^2002 , 1000]], 10^2] results in 9, so the last two digits before the decimal point are 09. With a slight modification we can find the first n digits after the decimal > point. Simply find the last n digits before the decimal point of 10^n times > the number. Mod[Floor[ N[10^2 (Sqrt[2] + Sqrt[3])^2002 , 1000]], 10^2] results in 99, so these are the digits you are interested in. But there is something curious about this number. Mod[Floor[N[10^1000 (Sqrt[2] + Sqrt[3])^2002 , 2000]], 10^1000] results in 996 digits 9 followed by 7405. You can also play with the following command, resulting in the digits around > the decimal point: Mod[ N[(Sqrt[2] + Sqrt[3])^2002, 2300], 10^6] The decimal expansion of (Sqrt[2]+Sqrt[3])^2002 contains a sequence of > 997 consecutive digits 9. Do you have any idea why? Not unusual. Try other even exponents and there should be large blocks of 9s . Smaller number for small exponents. Some kind of propagation of 9s as the exponent grows. Try other odd values for 3 and the same thing happens for some values. Also for other values for 2. ==== >-----Original Message----- A commonly used symbol for the Floor function is a square >bracket with the >upper indents removed. Is this symbol part of Mathematica's >built-in symbols? I >think not, so then, can it be constructed by hand. yes it is already built-in! Bring up the palette with menu: File > Palettes > CompleteCharachters; there is a section Operators > General, where you'll find what you want. Alternatively type 'esc' l f 'esc' 'esc' r f 'esc', or instead of the esc-sequence use the corresponding mark-ups. For output try TraditionalForm. If you don't like that for all of your output, you may just add a formatting rule for Floor: In[7]:= Unprotect[Floor] In[8]:= Floor /: MakeBoxes[Floor[expr_], StandardForm] := RowBox[{[LeftFloor], MakeBoxes[expr, StandardForm], [RightFloor]}] In[9]:= Protect[Floor] In[11]:= Floor /@ ([Pi] + [Lambda]) Out[11]= 3 + [LeftFloor][Lambda][RightFloor] Same thing with Ceiling, BTW. ==== >-----Original Message----- >Is it possible to configure a style sheet so that functions or >data shown in >graphs are coloured (colored :)) in working but black in >printing...without >actually re-executing the function that created the graphics? thanks Mike Mike, if for some reason you can't configure your color printer when printing, you may rerender your graphics in Mathematica with different options, such avoiding costly recalculation. You may either use something special, as in... In[1]:= ContourPlot[Sin[x]^2 Sin[y]^2, {x, 0, Pi}, {y, 0, Pi}, ColorFunction -> Hue, ContourStyle -> GrayLevel[1], PlotPoints -> 100] Out[1]= [SkeletonIndicator]ContourGraphics[SkeletonIndicator] In[2]:= Show[%, ColorFunction -> (GrayLevel[1 - #] &)] ...using a B/W color function. Or in general you may use the option ColorOutput: In[3]:= Plot3D[Sin[x]^2 Sin[y]^2, {x, 0, Pi}, {y, 0, Pi}] Out[3]= [SkeletonIndicator]SurfaceGraphics[SkeletonIndicator] ... In[5]:= Show[%%, ColorOutput -> GrayLevel] -- ==== I'm trying to write a fast empirical cummulative distribution function (CDF). Empirical CDFs are step functions that can be expressed in terms of a Which statement. For example, given the list of observations {1, 2, 3}, f = Which[# < 1, 0, # < 2, 1/3, # < 3, 2/3, True, 1]& is the empirical CDF. Note that f /@ {1, 2, 3} returns {1/3, 2/3, 1} and f is continuous from the right. When the number of observations is large, the Which statement evaluates fairly slowly (even if it has been Compiled). Since InterpolationFunction evaluates so much faster in general, I've tried to use Interpolation with InterpolationOrder -> 0. The problem is that the resulting InterpolatingFunction doesn't behave the way (I think) it ought to. For example, let g = Interpolation[{{1, 1/3}, {2, 2/3}, {3, 1}}, InterpolationOrder -> 0] Then, g /@ {1, 2, 3} returns {2/3, 2/3, 1} instead of {1/3, 2/3, 1}. In addition, g is continuous from the left rather than from the right. Obviously I am not aware of the considerations that went into determining the behavior of InterpolationFunction when InterpolationOrder -> 0. So I have two questions: (1) Does anyone have any opinions about how InterpolatingFunction ought to behave with InterpolationOrder -> 0? (2) Does anyone have a faster way to evaluate an empirical CDF than a compiled Which function? By the way, here's my current version: CompileEmpiricalCDF[list_?(VectorQ[#, NumericQ] &)] := Block[{x}, Compile[{{x, _Real}}, Evaluate[ Which @@ Flatten[ Append[ Transpose[{ Thread[x < Sort[list]], Range[0, 1 - 1/#, 1/#] & @ Length[list] }], {True, 1}]] ]]] --Mark ==== I am new to numerical computing. I have a equation dQ/dz = exp(-i*4*a^2*z - i*pi/4)*f(z,t)/sqrt(4*pi*z) where z is the position range from 0 to 10, t is the time range from 0 to 10 f(z,t) = integration of exp(i*(t-t1-z)^2/z)*g(t1) from t1=-100 to 100 g(t1) = d^3 sech(t1) / dt^3 I would like to solve Q(z,t), however, I am a new to mathematica. Can anyone help me? ==== I know that there is an input form. I want the output to have the same form rather than look like Floor[x]. Sorry, A commonly used symbol for the Floor function is a square bracket with the > upper indents removed. Is this symbol part of Mathematica's built-in symbols? I > think not, so then, can it be constructed by hand. Reply-To: kuska@informatik.uni-leipzig.de ==== does MakeBoxes[Floor[x_], fmt_:StandardForm] := RowBox[{[LeftFloor], ToBoxes[x, fmt], [RightFloor]}] help you ? I know that there is an input form. I want the output to have the same > form rather than look like Floor[x]. Sorry, A commonly used symbol for the Floor function is a square bracket with the upper indents removed. Is this symbol part of Mathematica's built-in symbols? I think not, so then, can it be constructed by hand. ==== Recent threads about word processing and typesetting have got me dreaming again... Wouldn't it be nice if someone developed a package that provided a function like this: DisplayTeX[ some TeX code ] which would cause the kernel to generate PostScript for the typeset text, to be displayed by the front end? I doubt that WRI will ever consider this, but surely there's someone out there who's both a TeXpert and a Mathematica guru who can do it. No doubt there's money to be made. Probably just a pipe dream... --- Reply-To: kuska@informatik.uni-leipzig.de ==== you mean the inverse of TeXForm[] ? For what ? A TeXpert && Mathematica guru && MathLink expert has a TeX frontend, that translate a TeXNotebook (that is a TeX file with some Mathematica Input/Output environments) and send the input cells to a kernel and insert the output into the final TeX file ... The most of the remaining work like creating hyperlinks to references and figures is done my LaTeX and pdfLaTeX, with help of CTAN and some macros one generate almost every layout *and* I have more than 20 books about TeX/LaTeX (including Don Knuths excelent manuals) but I have not a single book about the Mathematica Frontend. The way is not to teach TeX to the frontend, the way is to teach TeX a bit Mathematica. And a TeXpert will never switch from his beloved TeX to the Frontend and it's typesetting -- that's why he is a TeXpert. Recent threads about word processing and typesetting have got me > dreaming again... Wouldn't it be nice if someone developed a package that provided a > function like this: DisplayTeX[ some TeX code ] which would cause the kernel to generate PostScript for the typeset > text, to be displayed by the front end? I doubt that WRI will ever consider this, but surely there's someone out > there who's both a TeXpert and a Mathematica guru who can do it. No > doubt there's money to be made. Probably just a pipe dream... --- > ==== Mike, I don't think so. Style sheets do not control the styles used within an output graphics cell. I think you will have to set the plot style as a statement in the notebook, change it and re-evaluate for printing. This might not be all that inconvenient because you probably won't be printing that often. Park djmp@earthlink.net http://home.earthlink.net/~djmp/ Sender: steve@smc.vnet.net Approved: Steven M. Christensen , Moderator ==== Dear MathGroup, To plot FrameTicks outward with the same style as the Frame, I use ExtendGraphics in the following way. Needs[ExtendGraphics`Ticks`]; tf[min_,max_] := TickFunction[min, max, MajorStyle -> {Thickness[0.003]}, MajorLength - > {0, 0.01}, MinorStyle -> {Thickness[0.003]}, MinorLength -> {0, 0.005}] Plot [-2.4 (x+6) (x-8), {x, -6.2, 9.2}, FrameTicks -> {tf, tf, None, None}, AspectRatio -> 18/25, Frame -> True, DefaultFont -> {Helvetica-Bold, 12}, Axes -> None, FrameStyle -> Thickness[0.003], FrameLabel -> {x, y, None, None}, PlotStyle -> Thickness[0.006], ImageSize -> 504]; (instead of the same Plot statement with FrameTicks -> {Auatic, Auatic, None, None}) This arrangement of major and minor ticks is clearly unacceptable. The minor ticks do not divide intervals of adjacent labelled major ticks into equal parts. Are there other solutions? I have noticed that some major statistical (plotting) programs don't do minor ticks at all. Aldenberg RIVM Bilthoven Netherlands ==== I need to create an adjacency matrix from my data, which is currently in the form of a .txt file and is basically a two column incidence list. For example: 1 A 1 B 2 B 3 C . . . . . . m n Where 1 to m represent actors and A to n represent events. My goal is to have an (m x m) matrix where cell i,j equals 1 if two actors are incident to the same event (in the sample above, 1 and 2 are both incident to B) and 0 otherwise (w/ zeros on the diagonal). I'm new to Mathmatica, and so I'm on the steep part of the learning curve ... All I've been able to figure out so far is how to get my incidence list into the program using Import[filename.txt]. But then what? How do I convert to the adjacency matrix? I've found the ToAdjacencyMatrix[] command in DiscreteMath`Combinatorica`, but I can't seem to get it to work ... ==== with In[]:=lst = {{1, A}, {1, B}, {2 , B}, {3, C}, {3, D}, {1, D}}; and In[]:= AdjacenceMatrix[lst : {{_, _} ..}] := Module[ {actors,events adj}, {events, actors} = Union /@ Transpose[lst]; adj = Table[0, {Length[events]}, {Length[actors]}]; Scan[(Part[adj, Sequence @@ #] = 1) &, lst /. MapIndexed[Rule[#1, First[#2]] &, actors]]; adj ] you get In[]:=AdjacenceMatrix[lst] Out[]={{1, 1, 0, 1}, {0, 1, 0, 0}, {0, 0, 1, 1}} I need to create an adjacency matrix from my data, which is currently in > the form of a .txt file and is basically a two column incidence list. > For example: 1 A > 1 B > 2 B > 3 C > . . > m n Where 1 to m represent actors and A to n represent events. My goal is to > have an (m x m) matrix where cell i,j equals 1 if two actors are > incident to the same event (in the sample above, 1 and 2 are both > incident to B) and 0 otherwise (w/ zeros on the diagonal). I'm new to Mathmatica, and so I'm on the steep part of the learning > curve ... All I've been able to figure out so far is how to get my > incidence list into the program using Import[filename.txt]. But then > what? How do I convert to the adjacency matrix? I've found the > ToAdjacencyMatrix[] command in DiscreteMath`Combinatorica`, but I can't > seem to get it to work ... ==== >Probably a rather simple question: what is the easiest way to open a >notebook, execute it, then quit, from the command line? I would like to be >able to do this from a Makefile. actually, this is not as simple as it should be and one needs the extremely useful JLink to really manipulate Notebooks from the similar under MacOSX and even Windows. Maybe others can try this. ================================================= Execute this in the FrontEnd to create a test notebook, then quit the FrontEnd completely: NotebookSave[NotebookPut[Notebook[{ Cell[< Plot3D[Sin[x]*Cos[y], {x, 0, 10}, {y, 0, 10}, PlotLabel -> >, Input], Cell[Integrate[Log[1 - x]^2/(1 + x), {x, 0, 1}], Input] } ]],/tmp/test.nb] (* **************************** *) Next, create a file t.m : [rolf@uranus tmp]$ cat t.m nb = /tmp/test.nb; (* ***************************** *) <, Input], Cell[NotebookSave[en, Interactive -> False], Input]}]]; SelectionMove[m, All, Notebook]; SelectionEvaluate[m]; ]; Pause[1]; (* **************************************** *) Now you can run it like this in the background: [rolf@uranus tmp]$ math < t.m > t.out & This will open /tmp/test.nb , evaluate it and save it under the same name. (* ******************************* *) Maybe there is an easier and cleaner way, but the problem here is that the kernel and the FrontEnd are really two different programs, i.e., the kernel is not waiting until the FrontEnd finished certain tasks (like NotebookSave), i.e., think of it as threads. Therefore I used the trick with a second steering notebook. (One usually has similar problems when writing more intricate buttons.) If you are not, e.g. you work remotely at a Unix-box without an virtual frame buffer (Xvfb), like in webMathematica. ==== > I'm curious about the different efficiencies of Export vs. Put; Import vs. > Get. Eg. Saving a list called c with approx. 16,000 elements, each element > being a three element list of reals. Export[ac.dat,c] takes forever and consumes heaps of kernel memory c>>ac.dat all done in an instant (relatively speaking) and without the > transient burst of kernel memory usage. Ditto Import and Get. Can anyone explain this? The comparison is that of apples and oranges. Get[] and Put[] deal with reading and writing Mathematica expressions in InputForm syntax. They are implemented in C, so they are fairly fast, and the target format is human readable. However the output might not be readily parsable by other computer programs. The functions may be used with Mathematica expressions of all kinds. In the event that human readability and platform independence are not important, one can use DumpSave[] in lieu of Put[]. Export[] and Import[] for the cases you describe are using the Table target format. They are implemented in top-level Mathematica code, so they are not as fast as Get[] and Put[]. The target is a regularly formatted array of data, so not all Mathematica expressions are appropriate for this kind of conversion. In the case of Export[], the expression is formatted completely in memory using TableForm[] before being saved out to file. In the case of Import[], there are a number of heuristics that are applied to each field to determine whether the field should be converted to a number or a string. This should explain why they are slower and require more memory. If you are dealing with a large collection of machine precision numbers, then neither of these approaches may be suitable for your purposes. You might be better off using the free add-on FastBinaryFiles, which is available up on MathSource at this URL: readable by C programs. -- User Interface Programmer paulh@wolfram.com Wolfram Research, Inc. ==== > A commonly used symbol for the Floor function is a square bracket with > the upper indents removed. Is this symbol part of Mathematica's > built-in symbols? I think not, so then, can it be constructed by > hand. Is this what you're describing? In[1]:= [LeftFloor]x[RightFloor] Out[1]= Floor[x] This notation is documented in Section A.2.6 of _The Mathematica Book_ (Fourth Edition). -- ==== , Check pages 959 and 960 in Section 3.10.4 of The Mathematica Book. You can use LeftFloor and RightFloor as bracketing operators, which in turn are intrepreted as Floor. [LeftFloor] 3.2 [RightFloor] 3 ==== Adding my two cents to: > There are many cases in graphics, and otherwise, where it is useful to > obtain two orthogonal unit vectors to a given vector. I know a number of > ways to do it, but they all seem to be slightly inelegant. I thought I would > pose the problem to MathGroup. Who has the most elegant Mathematica > routine... To this I would like to add a criterion of smoothness. Armed with a second vector b not parallel to the given vector a, it's a trivial matter to orthogonalize b WRT a by Gram-Schmidt and then form the third vector c = a x b. (Normalize as needed.) I don't need more elegance that this, but I would like a scheme to select the vector b that results in a triad {a,b,c} that various smoothly as a varies over all possible directions. Each of my attempts to date involve a branched algorithm and jumps in the resulting triad for certain small changes in a. To 's call for elegance I add a call for smoothness. Burton ==== Does anybody have Benchmarks for Mathematica 4.2 to share? If possible the same or similar machine and different OSes. to the Windows version. I tried http://www2.staff.fh-vorarlberg.ac.at/~ku/karl/timings40.html value of 6.4. ---more details-- Times = [InvisibleSpace]{1.18, 1.14, 1.08, 0.55, 1.18, 2.48, 0.74, 0.79, 0.41, 0.11, 0.38, 0.96, 0.97, 1.26, 1.11} Time = [InvisibleSpace]17. Benchmark = [InvisibleSpace]6.40784 ------------------------------------------------ I would like to see how Mathematica 4.2 improved in speed compared to Mathematica 4 or 3. ==== Use Text statements. Click of a coordinate point on the plot, copy it and paste it into the Text statement. Here is an example. Show[Graphics[ {Circle[{1, 1}, 2], Text[A, {1.75659, 2.03292}], Circle[{0, 0}, 2], Text[B, {-1.38818, -0.0197372}], Circle[{1, -1}, 2], Text[C, {1.57238, -1.99344}]}], AspectRatio -> Auatic, ImageSize -> 400]; To click off a coordinate: 1) Select the plot by putting the cursor in it and clicking. 2) Press and hold Ctrl. When you move the cursor you will obtain cross hairs. 3) Move the cursor to where you want to put the text and left click. 4) Right click and use Copy. (Or use Ctrl-C or the menu). 5) Paste the copied coordinate as the second argument in the Text statement. Park djmp@earthlink.net http://home.earthlink.net/~djmp/ ==== Mark, Just add a PlotRange option and that each penino will be plotted to the same scale. Show[ GraphicsArray[ Partition[ Table[ Graphics[ Map[Rectangle[#, (# + {1, 1})] &, peninoes[[i]]], AspectRatio -> Auatic, PlotRange -> {{0, 5}, {0, 5}} ], {i, 12}], 6] ] ] Park djmp@earthlink.net http://home.earthlink.net/~djmp/ {{0, 0}, {1, 0}, {2, 0}, {3, 0}, {4, 0}}, (*I*) {{0, 0}, {1, 0}, {2, 0}, {3, 0}, {0, 1}}, (*L*) {{0, 0}, {1, 0}, {2, 0}, {2, 1}, {3, 1}}, (*N*) {{0, 0}, {1, 0}, {2, 0}, {0, 1}, {1, 1}}, (*P*) {{0, 0}, {1, 0}, {2, 0}, {1, 1}, {1, 2}}, (*T*) {{0, 0}, {1, 0}, {2, 0}, {0, 1}, {2, 1}}, (*U*) {{0, 0}, {1, 0}, {2, 0}, {0, 1}, {0, 2}}, (*V*) {{0, 0}, {1, 0}, {1, 1}, {2, 1}, {2, 2}}, (*W*) {{1, 0}, {0, 1}, {1, 1}, {2, 1}, {1, 2}}, (*X*) {{0, 0}, {1, 0}, {2, 0}, {3, 0}, {1, 1}}, (*Y*) {{0, 0}, {1, 0}, {1, 1}, {1, 2}, {2, 2}} (*Z*) } When I display them with the code below, the aspect ratio is correct *within each penino* but they are not scaled equivalently. I would appreciate a pointer to what I have done wrong. The only solution I have been able to develop so far effectively places each of the 12 peninoes at a different location in the plane, and then draws the whole region in one go. Show[ GraphicsArray[ Partition[ Table[ Graphics[ Map[Rectangle[#, (# + {1, 1})] &, peninoes[[i]]], AspectRatio -> Auatic ], {i, 12}], 6] ] ] -- ==== In my previous post, I proposed OrthogonalUnitVectors[v:{_, _, _}] := With[{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 ] }, #/Sqrt[#.#]& /@ {u, Cross[u,v]}] The trouble with this is that w ends up being a global variable. The only way I see around that is to use Module instead of With. (May as well put in a Return[$Failed] too.) 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]} ] ==== I have spent an enormous amount of time (far too much) on this question. Indeed, I have just completed a program that handles all sorts of piecewise defined functions. I am in the process of writing it up to offer it to all of you. In my opinion the best, cleansest method is to invoke the UnitStep function as illustrated in Park's solutions. To make everything cleaner, define a characteristic function Chi[x_,a_,b_] := UnitStep[x-a]-UnitStep[x-b] This function is continuous on the right, vanishes outside of [a,b) and is one in the half-open interval [a,b). Now your answer is x^2 Chi[x,0,6] + (x+1) Chi[x,6,Infinity] The integrator handles UnitSteps easily and smoothly. My program handles these cases easily and many more, including such oddities as Integrate[Abs[x],x] and (for amusement sake only) UnitStep[ Abs[x]-Sign[x]+3] Another virtue of this approach (also handled in my program) is the peculiar error messages and poor answers given by NIntegrate for simple piecewise continuous functions at jump discontinuites. The limit function also fails dismally on some examples where it really shouldn't I am rather naive about how to transmit programs to interested users, so if you would like a copy of my program and about 100 worked examples, I will be glad to send them to you IF YOU TELL ME HOW TO DO IT! Goldberg > In a message dated 8/31/02 1:58:36 AM, berlusconi_pagliusi@fis.unical.it I'd like to use Mathematica 4.0 to write a function having different expressions in different domain's intervals. Let's say: F[x_]= x^2 if 0=6 I know It's a stupid syntax problem, but I really do not know how/where to search the solution on the Mathematica Book Just for grins here are several methods: f1[x_/;x<=0] := 0 ; f1[x_/;0=6] := x+1; > f2[x_] := x^2*UnitStep[x]+(x+1-x^2)*UnitStep[x-6]; > f3[x_] := Which[ x<=0, 0, 0=6, x+1]; > f4[x_] := If[0=6, x+1,0]]; > f5[x_] := Switch[x, _?(#<=0&), 0, _?(0<#<6&), x^2, _?(#>=6&), x+1 ] > f6[x_?NumericQ] := {0,x^2,x+1}[[Position[ {x<=0,0=6}, True][[1,1]]]]; > Needs[Calculus`Integration`]; (* needed for definition of Boole *) > f7a[x_?NumericQ] := Evaluate[ (Boole /@ {0=6}). {x^2,x+1}]; > Off[Part::pspec]; f7b[x_?NumericQ] := Evaluate[ {0,x^2,x+1}[[1+Tr[Boole /@ {x>0, x>=6}]]]]; > f8[x_?NumericQ] := Cases[ {{x<=0, 0}, {0=6, x+1}}, {True, z_} :>z][[1]]; > f9[x_?NumericQ] := DeleteCases[ {{x<=0, 0}, {0=6, x+1}}, {False, z_}][[1,2]]; > f10[x_?NumericQ] := Select[ {{x<=0, 0}, {0=6, x+1}}, First[#]&][[1,2]]; > f11[x_?NumericQ] := Last[Sort[ {{x<=0, 0}, {0=6, x+1}}]][[2]]; > f12[x_?NumericQ] := Module[{n=1}, While[{x<6,x<0, False}[[n]], n++]; {x+1,x^2,0}[[n]]]; > Generating some test points: ts = {Random[Real,{-5,0}],0, Random[Real,{0,6}],6,Random[Real,{6,15}]}; > Checking the different representations Equal[(# /@ pts)& /@ {f1,f2,f3,f4,f5,f6,f7a,f7b,f8,f9,f10,f11,f12}] > True To pick a favorite, look at how the different definitions behave. > Of the definitions that evaluate with symbolic input, only f2 and f4 > simplify with assumptions > . For example, FullSimplify[#[x]& /@ {f2,f4}, 1 f2 through f5 respond immediately to differentiation > #'[x]& /@ {f2,f3,f4,f5} // Simplify //ColumnForm > Only f2 responds immediately to integration > Integrate[f2[x],x]//Simplify Consequently, f2 (UnitStep) appears to be the most versatile. ==== 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 Reply-To: kuska@informatik.uni-leipzig.de ==== the volume is 0 because you get two points Solve[{y == 4 - x^2, y == 3x}, {x, y}] {{y -> -12, x -> -4}, {y -> 3, x -> 1}} with {{y -> -12, x -> -4,z->20}, {y -> 3, x -> 1,z -> 5}} To get a nonzero volume you neeed a singel implicit equation like 0==f[x,y,z] or you need 3 parametric equations { x==f[1][u,v], y==f[2][u,v], z==f[3][u,v]} If you have an implicit equation you can compute the surface with MathGL3d's MVContourPlot3D[] (not with the ContourPlot3D[] from the standard addons, because you need a consistent oriented surface) and than you can triangulate the surface an sum the signed voulmes of the tetrahedrons build from the triangles and a fixed point in space. For parametric surfaces you can do the same with the polygons. 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 ==== Hugh Goyder and Park gave a most elegant function to find two vectors that are orthogonal to one vector in 3D. The key to coming up with the elegant solution is an understanding of Mathematica's NullSpace function. We can easily make the version from Hugh and much more general with the version below. ------------- OrthogonalUnitVectors[vect__?VectorQ]:= #/Sqrt[#.#]&/@NullSpace[{vect}] ------------- The version above will give a set of unit orthogonal vectors if given any number of vectors in any dimension. So besides giving it a 3D vector we can give it the following: OrthogonalUnitVectors[{2,1,0,-1,1}] or OrthogonalUnitVectors[{0,1,0,1/2,1},{1,0,-1,1/2}] ------------ But the short version above isn't very robust. (1) Clear[x,y,z];NullSpace[{{x,y,z}}] returns two vectors orthogonal to {x,y,z}, but the two vectors NullSpace returns aren't orthogonal to each other. So (OrthogonalUnitVectors) should only work with numeric vectors. (2) We should ensure all the vectors have the same dimension and length >1. I give a less concise version below that corrects these problems. ------------ OrthogonalUnitVectors[vect__?(VectorQ[#,NumericQ]&)]/; (SameQ@@Length/@{vect})&&(Length[First[{vect}]]>1):= #/Sqrt[#.#]&/@NullSpace[{vect}] -------------- Get Mathematica tips, tricks from http://www.verbeia.com/mathematica/tips/Tricks.html ==== First of all, I would like to thanks to all for trying to help. 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 bot 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? ============================================================================ ====== of the individual or entity to which they are addressed. Any disclosure, copying, distribution and diversion contrary to the applicable export control laws and regulations including US Export Administration Regulations is strictly prohibited. and do not disclose it to others. Please notify the postmaster@hitachi.com.my of the delivery error by replying to this message and then delete it from your ============================================================================ ==== you mean: In[]:=Needs[Calculus`Integration`] Integrate[Boole[x^2 + 4 - z > 0 && 4 - x^2 - y > 0 && 3 - x - y < 0 && z > 0], {x, -1, 2}, {y, 1, 5}, {z, -1, 10}]] Out[]=(15*Sqrt[5])/4 and to plot the volume: In[]:=Get[MathGL3d`OpenGLViewer`] In[]:=MVContourPlot3D[ If[1 == Boole[x^2 + 4 - z > 0 && 4 - x^2 - y > 0 && 3 - x - y < 0 && z > 0], 1., -1.], {x, -1, 2}, {y, 1, 5}, {z, -0.1, 7}, Contours -> {0.}, PlotPoints -> 64, MVNewScene -> True] First of all, I would like to thanks to all for trying to help. > 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 bot 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? ============================================================================= ===== of the individual or entity to which they are addressed. Any disclosure, copying, > distribution and diversion contrary to the applicable export control laws and > regulations including US Export Administration Regulations is strictly prohibited. and do not disclose it to others. Please notify the postmaster@hitachi.com.my > of the delivery error by replying to this message and then delete it from your > ============================================================================= ===== ==== 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 3 notation is composed of k ones. By doubling those counterweights, you can precisely express all the EVEN integers up to TWICE that limit, or -1 + 3^k. If you know that unknowns are limited to the integers from 1 to 3^k, this allows you to precisely weight the even numbers and bracket the odd numbers. Hence counterweights 2, 6, ... 2*3^(k-1) are sufficient for unknowns up to 3^k. Notice that this is an efficient coding scheme: Each counterweight is multiplied by -1, 0, or 1, so there can be at most 3^k distinct counterweight sums. More than half are negative or zero, however; the positive sums and differences number at most (3^k - 1)/2. Since we need to enumerate (and CAN enumerate with even weights) only the even numbers and because we get the maximum unknown 3^k by elimination, we cover up to twice as many as we can enumerate, plus one. Twice (3^k - 1)/2 plus one is 3^k, so we're getting the maximum coverage possible for a given number of counterweights. -----Original Message----- For these 40 positive integers to match 1 to 40 correspondingly, the biggest should 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 a nature > 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 to solve > 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 this math > 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 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}} > ==== >>But my actual output reproduces the input form, i.e., it is simply f[...listA,newEntry,...]. That means your arguments don't fit a pattern for which the function is defined. 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 of arguments passed to the function. You haven't told us what the pattern is 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 exhausive tests, each test with its own action that modifies the concrete list subsituted for listA_. For instance, if there were only two tests the rhs above would be one that puts the newEntry at 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, ...]. The output should be the appropriately modified list. But my actual output reproduces the input form, i.e., it is simply f[...listA,newEntry,...]. Yet, the program does do something, as the use of evaluate in placeshows: Suppose that test1 is satisfied by the particular values of the arguments. Then when I apply evaluate in place to the lhs of the action that is supposed to occur when test1 is true, in fact the value of Insert[listA,newEntry,1] is correct, but the rhs gives only the initial value of listA and not the modified value that is on the rhs. Concretely, suppose that in the function argument listA is {121} and newEntry is 200. Then after evaluation of the function (input to the kernel), that line reads, according to evaluate in place: true, {121} = {200, 121} So it seems that only part of the correct action was taken: 200 was put at the start of listA. But the assignment of this value -- the extended list -- to be the new value of listA did not take place. I then tried using a new name for the modified list, substituting this program line for test1: test1, newListA = Insert[listA, newEntry, 1] What I get here, applying evaluate in place after evaluating the function with concrete values as above is: true, newListA = {200, 121} Trace doesn't help because all I get back is the function name with its concrete arguments. ==== I have installed 128 MB RAM extra in my computer. I had 64 before, so now I have 3 times as much. I expected some calculations that took very long before, to be executed in less time. To the contrary, my computer still goes to the hard drive fro memory for its calculations, and Mathematica stops with an out of memory message. Windows is recognizing this new memory. I dont know if Mathematica is. How can I find out about this? ==== > It seems there is an error in the BinCounts function of Mathematica's > standard < returns one more than the expected number of histogram bins. An example > of this is shown in the notebook code below where the comment in the > last notebook cell explicitly describes the problem. If anyone can > explain this behavior it would be appreciated, otherwise Wolfram > Research should be made aware of this feature of the function when it > is called as: BinCounts[{x1, x2, ...}, {xmin, xmax, dx}] A developer of Mathematica's Standard Packages has requested that this response be posted to the forum: [begin developer's comments] Note that if you don't apply N when computing datBinWidth, you will always get the expected number of bins; i.e., datBinWidth = (datMax - datMin)/datBins (assuming datMax and datMin are also exact numbers). Given input {xmin, xmax, dx}, BinCounts determines the number of bins to be computed via Ceiling[(xmax - xmin)/dx]. For numericalized values, it's possible for dx to end up such that the Ceiling forces an additional bin to be added; presumably, it's better to have too many rather than too few bins -- for example, if the dx was deliberately given such that the xmax is not at an integer bin boundary... If you use exact values for bin bounds and increment, then there will be no problem. [end developer's comments] ==== I would like to create a 3D chart from an ASCII file using Mathematica 3.0. Can anyone give me some help? ==== You might want to tell us how to find A.10. Anyway, it's easier to simply go to Floor in the help browser (in version 4.2, anyway). -----Original Message----- ==== I doubt there's any money to be made here. Lots of people have created very valuable additions to Mathematica, but they're all being given away, so far as I can tell. It's too bad, too. And here we are, giving away our time. Sigh... Bobby -----Original Message----- there who's both a TeXpert and a Mathematica guru who can do it. No doubt there's money to be made. Probably just a pipe dream... --- ==== OK, I give up; how would I find Chapter Advanced Mathematics/Mathematical and Other Notations/Operators in the Help Browser? -----Original Message----- > think not, so then, can it be constructed by hand. [LeftFloor] and [RightFloor] can be used for this. Look at Chapter Advanced Mathematics/Mathematical and Other Notations/ Operators. ==== > You might want to tell us how to find A.10. Open your fourth edition of _The Mathematica Book_ to page 1065. That's where part 10 of the Appendix starts. > Anyway, it's easier to simply go to Floor in the help browser (in > version 4.2, anyway). I suppose you're right. But I -- call me old-fashioned if you wish -- prefer to use the book. ==== I'll assume you have the information in a matrix like x: TableForm[x = {{1, A }, {1, B }, {2, B}, {3, D}, {4, D}, {5, C}}] First find the highest actor number (or use what you've input elsewhere): m = Last[Union[x[[All,1]]]] 5 Define the incidence function: f[x_, a_, b_] /; a == b := 0 f[x_, a_, b_] := If[{} $B!b(B Intersection[Cases[x, {a, y_} -> y], Cases[x, {b, y_} -> y]], 1, 0] Here's the incidence matrix: Array[f[x, #1, #2] &, {m, m}] {{0, 1, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 1, 0, 0}, { 0, 0, 0, 0, 0}} Once you have the incidence function, the only reason to store the incidence matrix is to avoid computing over again the same answers, and that can be accomplished in other ways. If you won't have other x matrices, it's convenient to define f this way: f[a_, b_] /; a == b := f[a, b] = 0 f[a_, b_] /; a < b := f[a, b] = If[{} $B!b(B Intersection[Cases[x, {a, y_} -> y], Cases[x, {b, y_} -> y]], 1, 0] f[a_, b_] := f[b, a] Array[f, {m, m}] {{0, 1, 0, 0, 0}, {1, 0, 0, 0, 0}, {0, 0, 0, 1, 0}, {0, 0, 1, 0, 0}, {0, 0, 0, 0, 0}} Whenever you would want the {j,k} element of the adjacency matrix, just use f[j,k] instead. Each pair is computed only once. In this version, I've taken advantage of the symmetry of the problem, too, to cut the work in half. -----Original Message----- . . . . m n Where 1 to m represent actors and A to n represent events. My goal is to have an (m x m) matrix where cell i,j equals 1 if two actors are incident to the same event (in the sample above, 1 and 2 are both incident to B) and 0 otherwise (w/ zeros on the diagonal). I'm new to Mathmatica, and so I'm on the steep part of the learning curve ... All I've been able to figure out so far is how to get my incidence list into the program using Import[filename.txt]. But then what? How do I convert to the adjacency matrix? I've found the ToAdjacencyMatrix[] command in DiscreteMath`Combinatorica`, but I can't seem to get it to work ... ==== Well, I don't know how fast, but it is fairly simple, anyway. Suppose you have a series of values s for which you wish to obtain the edf. In[1]:= s = Table[Random[Integer, {0, 3}], {10}] Out[1]= {2,2,1,2,0,0,1,1,1,3} If no specification is made about their position on the x-axis, we assume that they correspond to the integers from 1 to 10. What we have then is the collection of pairs In[2]:= porig = Transpose[{Range[10], s}] Out[2]= {{1, 2}, {2, 2}, {3, 1}, {4, 2}, {5, 0}, {6, 0}, {7, 1}, {8, 1}, {9, 1}, {10, 3}} The edf gives, for each x, the proportion of points in s that are less than or equal to x, for all x. We obtain these proportions through the cumulative sums In[3]:= N[CumulativeSums[s]/Plus @@ s] Out[3]= {0.153846,0.307692,0.384615,0.538462,0.538462,0.538462,0.615385,0.692308,0. 769231,1.} so that for each of the pairs (x, y) below, y gives the proportion of points in s that are less than or equal to x: In[4]:= cumporig = Transpose[{Range[10], N[CumulativeSums[s]/Plus @@ s]}] Out[4]= {{1, 0.15384615384615385}, {2, 0.3076923076923077}, {3, 0.38461538461538464}, {4, 0.5384615384615384}, {5, 0.5384615384615384}, {6, 0.5384615384615384}, {7, 0.6153846153846154}, {8, 0.6923076923076923}, {9, 0.7692307692307693}, {10, 1.}} Now shift the x values one unit to the left, by dropping the last value and prepending 0 to them: In[5]:= ps = Transpose[{Prepend[Drop[Range[1, 10], -1], 0], CumulativeSums[s]/Plus @@ s}] Out[5]= {{0, 2/13}, {1, 4/13}, {2, 5/13}, {3, 7/13}, {4, 7/13}, {5, 7/13}, {6, 8/13}, {7, 9/13}, {8, 10/13}, {9, 1}} Then use Interpolation on this shifted set of points: In[6]:= ips=Interpolation[ps,InterpolationOrder[Rule]0] Out[6]= InterpolatingFunction[{{0,9}},<>] ips[x-1] is the edf you are looking for, as you may check by plotting it and displaying in the same graph together with the ListPlot of cumporig above. as Garza Mexico City ----- Original Message ----- > and f is continuous from the right. When the number of observations is large, the Which statement > evaluates fairly slowly (even if it has been Compiled). Since > InterpolationFunction evaluates so much faster in general, I've tried > to use Interpolation with InterpolationOrder -> 0. The problem is that > the resulting InterpolatingFunction doesn't behave the way (I think) > it ought to. For example, let g = Interpolation[{{1, 1/3}, {2, 2/3}, {3, 1}}, InterpolationOrder - 0] Then, g /@ {1, 2, 3} returns {2/3, 2/3, 1} instead of {1/3, 2/3, 1}. > In addition, g is continuous from the left rather than from the right. Obviously I am not aware of the considerations that went into > determining the behavior of InterpolationFunction when > InterpolationOrder -> 0. So I have two questions: (1) Does anyone have any opinions about how InterpolatingFunction > ought to behave with InterpolationOrder -> 0? (2) Does anyone have a faster way to evaluate an empirical CDF than a > compiled Which function? By the way, here's my current version: CompileEmpiricalCDF[list_?(VectorQ[#, NumericQ] &)] := > Block[{x}, Compile[{{x, _Real}}, Evaluate[ > Which @@ Flatten[ > Append[ > Transpose[{ > Thread[x < Sort[list]], > Range[0, 1 - 1/#, 1/#] & @ Length[list] > }], > {True, 1}]] > ]]] --Mark ==== I have no opinion on the behavior associated with InterpolationOrder->0, except that it should be DOCUMENTED, but isn't. Meanwhile, try this: lst = {{1, 1/3}, {2, 2/3}, {3, 1}}; ClearAll[empiricalCDF] empiricalCDF[{x_List, y_List}] := Compile[ {{z, _Real}}, Evaluate[ Which @@ Flatten[ Transpose[ {(z < #1 & ) /@ x, y}]]]] empiricalCDF[v:{{_, _}..}] := empiricalCDF[v] = empiricalCDF[ ({Join[#1[[1]], {Infinity}], Join[{0}, #1[[ 2]]]} & )[Transpose[ lst]]] empiricalCDF[lst] Plot[empiricalCDF[lst][x], {x, 1, 3}]; I split the work into two definitions for readability. If you evaluate: ?empiricalCDF after doing the plot above, you'll see that empiricalCDF[{{1, 1/3}, {2, 2/3}, {3, 1}}] has been compiled and saved for later use, and that it takes precedence over the SetDelayed rules listed after it. Hence, the compilation only occurs once for each list. As written, your CompileEmpiricalCDF would be compiled all over again every time it's used -- for each and every point you plot -- so of course it's slow. There are other problems, too. For instance, the pattern list_?(VectorQ[#, NumericQ] &) makes no sense at all. Maybe you meant list_?And[VectorQ[#],NumericQ[#]]& -----Original Message----- evaluates fairly slowly (even if it has been Compiled). Since InterpolationFunction evaluates so much faster in general, I've tried to use Interpolation with InterpolationOrder -> 0. The problem is that the resulting InterpolatingFunction doesn't behave the way (I think) it ought to. For example, let g = Interpolation[{{1, 1/3}, {2, 2/3}, {3, 1}}, InterpolationOrder -> 0] Then, g /@ {1, 2, 3} returns {2/3, 2/3, 1} instead of {1/3, 2/3, 1}. In addition, g is continuous from the left rather than from the right. Obviously I am not aware of the considerations that went into determining the behavior of InterpolationFunction when InterpolationOrder -> 0. So I have two questions: (1) Does anyone have any opinions about how InterpolatingFunction ought to behave with InterpolationOrder -> 0? (2) Does anyone have a faster way to evaluate an empirical CDF than a compiled Which function? By the way, here's my current version: CompileEmpiricalCDF[list_?(VectorQ[#, NumericQ] &)] := Block[{x}, Compile[{{x, _Real}}, Evaluate[ Which @@ Flatten[ Append[ Transpose[{ Thread[x < Sort[list]], Range[0, 1 - 1/#, 1/#] & @ Length[list] }], {True, 1}]] ]]] --Mark ==== Simplify: 5nth sqrt (3)/sqrt (6) 5nth Sqrt(3) ------------ Sqrt(6) 1st.. is this the correct notation for use on a computer 2nd.. how is the solution solved, step by step please. ==== I want to make a list of all symbols in the Global context, as in Names[Global`*] and compute a ByteCount for each symbol's OwnValues -- without evaluating the symbols. It seems possible in principle, but I haven't found a way. ==== with the following Java comands you can execute any Mathematica commands: ml.putFunction(EnterTextPacket, 1); ml.put(string); See JLinkUserGuide p. 164. ----- Original Message ----- > I've checked the documentation on Wolfram, but did not find the 'hook' > for ReadList from java. Any help would be appreciated. thanks. ==== I want to use the capability of mathematica from java. This is what I want to do: I want to read a data set and use ReadList to structure it into an array and then plot that data set with the output as a .gif. Any ideas? I've checked the documentation on Wolfram, but did not find the 'hook' for ReadList from java. Any help would be appreciated. thanks. Pete ==== Dear Mark, I suggest trying the following code before you put further energy in speeding up your functions. My code is very short and seems to be fast in my first test with a random array of 100000 integers. cdf[li_List]:= FoldList[#1+Length[#2]&, 0.0, Split[Sort[li]]]/Length[li] t=Table[Random[Integer,{1,1000}],{100000}]; Timing[cdf[t];] {0.44 Second,Null} Probably the speed could be increased further generating a compiled version. But this would require additional programming effort. Then you had to write a function which counts the number of equal data 'by hand' while scanning through the data list. Johannes > I'm trying to write a fast empirical cummulative distribution function > (CDF). Empirical CDFs are step functions that can be expressed in > terms of a Which statement. For example, given the list of > observations {1, 2, 3}, f = Which[# < 1, 0, # < 2, 1/3, # < 3, 2/3, True, 1]& is the empirical CDF. Note that f /@ {1, 2, 3} returns {1/3, 2/3, 1} > and f is continuous from the right. When the number of observations is large, the Which statement > evaluates fairly slowly (even if it has been Compiled). Since > InterpolationFunction evaluates so much faster in general, I've tried > to use Interpolation with InterpolationOrder -> 0. The problem is that > the resulting InterpolatingFunction doesn't behave the way (I think) > it ought to. For example, let g = Interpolation[{{1, 1/3}, {2, 2/3}, {3, 1}}, InterpolationOrder - 0] Then, g /@ {1, 2, 3} returns {2/3, 2/3, 1} instead of {1/3, 2/3, 1}. > In addition, g is continuous from the left rather than from the right. Obviously I am not aware of the considerations that went into > determining the behavior of InterpolationFunction when > InterpolationOrder -> 0. So I have two questions: (1) Does anyone have any opinions about how InterpolatingFunction > ought to behave with InterpolationOrder -> 0? (2) Does anyone have a faster way to evaluate an empirical CDF than a > compiled Which function? By the way, here's my current version: CompileEmpiricalCDF[list_?(VectorQ[#, NumericQ] &)] := > Block[{x}, Compile[{{x, _Real}}, Evaluate[ > Which @@ Flatten[ > Append[ > Transpose[{ > Thread[x < Sort[list]], > Range[0, 1 - 1/#, 1/#] & @ Length[list] > }], > {True, 1}]] > ]]] ====