title: Very long Refine/Solve batch run - is my code broken, or just complicated?

text: So I'm trying to run some Mathematica code in batch mode from my university cluster. Specifically, I'm trying to find the equilibria of a system of ordinary differential equations. Inspired by the answer to my original question on the Math site, I'm using a combination of Refine and Solve for this. Code for a simpler but related system, which runs swimmingly on either the cluster or my own machine is this: ![enter image description here](http://i.stack.imgur.com/huWKw.png) Since then, I've tried to run this adapted code for a considerably more complicated system: ![enter image description here](http://i.stack.imgur.com/vUVnV.png) This is taking _forever_. In the "longer than a day, devoured all the RAM on a 48 Gb cluster node and was well on its way to eating the RAM on a 96 Gb node before I killed it" sense of the word. Which is fine if the code is actually working - it's something of a hard problem for a computer to do, and if it just _takes_ a long time, so be it. My problem is there's also the distinct possibility I've done something wrong, created a runaway process, etc. and I'm not a good enough Mathematica programmer to spot it (I'm a rather poor Mathematica programmer in fact). Anyone able to spot anything clearly pathological about what I've done? Edit: Adding in editable/cut-and-pasteable code. Converting between the standard form and the input form has changed the code slightly, including swapping Simplify and Solve, but I have no idea if that matters.               Refine[(Simplify[#1, ϑ > 0 && ρp > 0 && σp >             0 && ρd > 0 && σd > 0 &&                  ρa > 0 && σa > 0 &&            1/ω > 0 && α > 0 && ψp > 0 && μp > 0 &&                  μa > 0 && θp > 0 && νup > 0 && θ >             0 && ζ > 0 && ψa > 0 &&                  θa > 0 && νua > 0 && γ > 0 && ϕ >             0 && κ > 0 && νcp > 0 &&                  τ > 0 && νca > 0] & )[          Solve[{ϑ*H - ρp*σp*            Cp*(Us/N) - ρd*σd*D*(Us/N)*ρa*σa*            Ca*(Us/N) == 0,               ρp*σp*Cp*(Us/N) + ρd*σd*            D*(Us/N)*ρa*σa*Ca*(Us/N) + ϑ*H == 0,               (1/ω)*Ua - α*Up - ρp*ψp*            Up*(H/N) - μp*σp*Up*(Cp/N) -                   μa*σa*Up*(Ca/N) - θp*            Up + νup*(θ*M + ζ*D) == 0,               α*Up - (1/ω)*Ua - ρa*ψa*            Ua*(H/N) - μp*σp*Ua*(Cp/N) -                   μa*σa*Ua*(Ca/N) - θa*            Ua + νua*(θ*M + ζ*D) == 0,               (1/ω)*Ca + γ*ϕ*D + ρp*ψp*            Up*(H/N) + μp*σp*Up*(Cp/N) +                   μa*σa*Up*(Ca/N) - α*Cp - κ*            Cp - θp*Cp + νcp*(θ*M + ζ*D) == 0,               α*Cp + γ*(1 - ϕ)*D + ρa*ψa*            Ua*(H/N) + μp*σp*Ua*(Cp/N) +                   μa*σa*Ua*(Ca/N) - (1/ω)*            Ca - κ*τ*Ca - θa*Ca +                   νca*(θ*M + ζ*D) ==           0, κ*Cp + κ*τ*Ca - γ*ϕ*            D - γ*(1 - ϕ)*D -                   ζ*D + νd*(θ*M + ζ*D) == 0,          Us + H + Up + Ua + Cp + Ca + D == 0,               Up + Ua + Cp + Ca + D == 0}, {Us, H, Up, Ua, Cp, Ca, D, N,          M}, Reals]], Null]      **Update** Per some suggestions below, am now trying the following:               Solve[{ϑ*H - ρp*σp*          Cp*(Us/N) - ρd*σd*D*(Us/N)*ρa*σa*          Ca*(Us/N) ==         0, ρp*σp*Cp*(Us/N) + ρd*σd*          D*(Us/N)*ρa*σa*Ca*(Us/N) + ϑ*H ==         0, (1/ω)*Ua - α*Up - ρp*ψp*          Up*(H/N) - μp*σp*Up*(Cp/N) - μa*σa*          Up*(Ca/N) - θp*Up + νup*(θ*M + ζ*D) ==         0, α*Up - (1/ω)*Ua - ρa*ψa*          Ua*(H/N) - μp*σp*Ua*(Cp/N) - μa*σa*          Ua*(Ca/N) - θa*Ua + νua*(θ*M + ζ*D) ==         0, (1/ω)*Ca + γ*ϕ*D + ρp*ψp*          Up*(H/N) + μp*σp*Up*(Cp/N) + μa*σa*          Up*(Ca/N) - α*Cp - κ*Cp - θp*          Cp + νcp*(θ*M + ζ*D) ==         0, α*Cp + γ*(1 - ϕ)*D + ρa*ψa*          Ua*(H/N) + μp*σp*Ua*(Cp/N) + μa*σa*          Ua*(Ca/N) - (1/ω)*Ca - κ*τ*Ca - θa*          Ca + νca*(θ*M + ζ*D) ==         0, κ*Cp + κ*τ*Ca - γ*ϕ*          D - γ*(1 - ϕ)*D - ζ*          D + νd*(θ*M + ζ*D) == 0,        N == Us + H + Up + Ua + Cp + Ca + D,        M == Up + Ua + Cp + Ca + D , ϑ > 0, ρp >         0, σp > 0, ρd > 0, σd > 0, ρa >         0, σa > 0, ω > 0, α > 0, ψp > 0, μp >         0, μa > 0, θp > 0, νup > 0, θ > 0, ζ >         0, ψa > 0, θa > 0, νua > 0, γ > 0, ϕ >         0, κ > 0, τ > 0, νca > 0, νd > 0}, {Us, H, Up,        Ua, Cp, Ca, D, N, M}, Reals]      It hasn't immediately solved the problem (let it run for ~hour on a local machine), but sending it to the cluster to see if it gets fixed.


title: Comparing Unequal Lists of Surnames to Establish Over or Under-Representation

text: I've been using _Mathematica_ to explore historical records. I have two long lists (one is 69k, the second 17k) of unequal length. They are like so (FWIW, the first list is the first 40 entries from the surnames of First World War attestation papers; the second list is the first 30 entries from the surnames of First World War war grave surnames - I have created percentage lists, to get a sense of how frequent their names were amongst all names):               atts={{"SMITH", 0.10031}, {"BROWN", 0.0555427}, {"WILSON",        0.0451438}, {"MCDONALD", 0.0411175}, {"JONES",        0.0375836}, {"TAYLOR", 0.0356139}, {"CAMPBELL",        0.033398}, {"JOHNSON", 0.0320656}, {"ANDERSON",        0.0303276}, {"WILLIAMS", 0.0298352}, {"THOMPSON",        0.0289952}, {"MARTIN", 0.0279669}, {"SCOTT", 0.0257944}, {"WHITE",        0.0250413}, {"STEWART", 0.0243606}, {"MILLER", 0.0240854}, {"CLARK",        0.0237088}, {"MACDONALD", 0.0224054}, {"WALKER",        0.0216233}, {"JOHNSTON", 0.0212322}, {"MOORE", 0.0212177}, {"YOUNG",        0.0211888}, {"ROBINSON", 0.0209426}, {"ROSS",        0.0196391}, {"WRIGHT", 0.018437}, {"MURRAY",        0.0182053}, {"ROBERTSON", 0.0181763}, {"MITCHELL",        0.0177273}, {"MCLEOD", 0.0176839}, {"JACKSON", 0.0175825}};              kia={{"SMITH", 0.0452787}, {"BROWN", 0.0260468}, {"WILSON",            0.0195784}, {"JONES", 0.0174415}, {"TAYLOR",            0.0169217}, {"CAMPBELL", 0.0166907}, {"McDONALD",            0.0146694}, {"ANDERSON", 0.0127057}, {"JOHNSON",            0.012417}, {"STEWART", 0.0123592}, {"SCOTT", 0.0120705}, {"WHITE",            0.0118394}, {"WILLIAMS", 0.0116662}, {"THOMPSON",            0.0115507}, {"MARTIN", 0.0107421}, {"WALKER", 0.0106844}, {"MOORE",            0.0103956}, {"CLARK", 0.0102224}, {"YOUNG", 0.0101646}, {"MILLER",            0.0101646}, {"ROSS", 0.00987583}, {"ROBERTSON",            0.00981808}, {"REID", 0.00958706}, {"WRIGHT", 0.00918279}, {"KING",            0.00912504}, {"MURRAY", 0.00895178}, {"MITCHELL",            0.008432}, {"JOHNSTON", 0.00825874}, {"HALL",            0.00820098}, {"MORRISON", 0.00814323}};      Does _Mathematica_ have a function that could help compare these two lists? I'm trying to think of a systematic way to establish that McDonald, for example, is under-represented in the war graves (4% of enlistees had the last name McDonald, whereas only 1.4% of dead did so). But since our lists are so long, would there be any way to do this systematically?


title: Problem with NDSolve

text: I'm solving a PDE that is partially elliptic. So I need to solve the equation iteratively. My problem is at the end of the code. I'm doing interactions manually in order to evaluate the code. At the second interaction, I don't know why NDSolve doesn't accept the new boundary condition. Please, find comments in the body of the code. Could some one help me?               vm = 5;     r0 = 127.25;     di = r0*2;     de = 273;     v = 0.000000464;     ρf = 958;     pr = 2.88;     kf = 0.656;     cf = 4186;     hi=35308.2;         αf = kf/(ρf*cf)     f=0.0174719     c1 = (-f*ρf*vm^2/(4*r0/1000))     tinf=298;     tme=373;     re = (vm*(2*r0/1000)/v)               ym[r_] := (((r0/1000) - r)*vm*(f/8)^(0.5))/v;     εm[r_] := v*(0.4*ym[r]/6)*(1 + (r/(r0/1000)))*(1 + 2*(r/(r0/1000))^2);            momento = NDSolve[{(1/r)*D[(v + εm[r])*r*D[u[r], r], r] == (1/ρf)*c1,u'[0.0000001] == 0, u[(r0/1000)] == 0},u[r], {r, 0.0000001, (r0/1000)}, MaxSteps -> 10000]; (*considerando difusividade de momento*);     mom[r_] = u[r] /. First[momento];          hi = (kf/(2*r0/1000))*((re - 1000)*pr*(f/8))/(1 + 12.7*(f/8)^0.5*(pr^(2/3) - 1))     umed=5.47              tmed[x_, t_] = 310;     ttf[x_, r_, t_] = tinf;     l = 160;     tempo = 80;          (*While[erro≥eps,*)          enf1 = NDSolve[{(ρf*cf*(D[tf[x, r, t], t] + mom[r]*D[tf[x, r, t], x])) == ((1/r)*        D[(r*(αf + εm[r])*D[tf[x, r, t], r]), r]),tf[x, r, 0] == (tme - tinf)*Exp[-1000*x] + tinf, tf[0, r, t] == tme, (D[tf[x, r, t], r] /. r -> 0.0000001) ==       0, -kf*(D[tf[x, r, t], r] /. r -> (r0/1000)) == (1 - Exp[-1000*x])*(1 - Exp[-1000*t])*hi*(tmed[x, t] - ttf[x, r0/1000, t])}, tf[x, r, t], {x, 0, l}, {r, 0.0000001, (r0/1000)}, {t, 0, tempo}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 13}}];      tff[x_, r_, t_] = tf[x, r, t] /. First[enf1];          Clear[med]          (*this calculates the average temperature*)     med[x_?NumericQ, t_?NumericQ] := (NIntegrate[(mom[r]*tff[x, r, t]*r), {r, 0.0000001,        r0/1000}])*(2/((r0/1000)^2*umed));          (*an evaluation of med[x,t] at x=0 shows that med is numeric at this position. this is valid for any time*)     med[0, tempo/2]     373.          (*so, it is created a new function that shall be numerically equal to med[x,t]*)     tmed1[x_?NumericQ, t_?NumericQ] := Evaluate[med[x, t]]      Up to here, everything is ok.               (*so, the code returns to NDSolve. The only thing that has changed here is tmed to tmed1. Why the mistake??*)          enf1 = NDSolve[{(ρf*cf*(D[tf[x, r, t], t] + mom[r]*D[tf[x, r, t], x])) == ((1/r)*        D[(r*(αf + εm[r])*D[tf[x, r, t], r]), r]), tf[x, r, 0] == (tme - tinf)*Exp[-1000*x] + tinf, tf[0, r, t] == tme, (D[tf[x, r, t], r] /. r -> 0.0000001) ==       0, -kf*(D[tf[x, r, t], r] /. r -> (r0/1000)) == (1 - Exp[-1000*x])*hi*(tmed1[x, t] - ttf[x, r0/1000, t])}, tf[x, r, t], {x, 0, l}, {r, 0.0000001, (r0/1000)}, {t, 0, tempo}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 13}}];     tff[x_, r_, t_] = tf[x, r, t] /. First[enf1];      > >     NDSolve::ndnum: Encountered non-numerical value for a derivative at x == > 0.`. >> >   PS: The original message on WOLFRAM community is here: http://community.wolfram.com/groups/-/m/t/317994?p_p_auth=NrBn8wS0


title: Optimise CPU-intensive list combining function

text: I have a list of numbers `list1` from which I want to create another list (`list2`) consisting of only those numbers that are equal to the sum of _any_ two members of `list1`, as long as the total is under 15. So for example if `list1 = {1,4,6,8}`, then I want `list2 = {2,5,7,8,9,10,12,14}`. I could do this by: `list2 = Union[Flatten[Table[Table[i + j, {i, list1}], {j, list1}]]]` which does give the right result, but then I run into two problems:   1. I have no idea how to filter the list so that it only displays numbers less than 15. I tried appending `/; x_ < 15` to the end, but this didn't do anything   2. This works fine for a list this small, but what I really want to do is scale this up to a list consisting of several thousand elements - for which it takes a very long time, and crashes my computer in the process. How could I make this function more efficient? Thank you


title: How to asynchronously download 100000 urls in an efficient way

text: `URLSaveAsynchronous` can download one url page, how to download 100000 urls in an efficient way? And better to with an whole progress bar. I've downloaded 90000 htmls a web-dictionary in mathematica this way, I think it's interesting.               progress= 0.;     progFunction[_, "progress", {dlnow_, dltotal_, _, _}]:= Quiet[progress = dlnow/dltotal]     Dynamic[ProgressIndicator[progress]]          cells = Table[Defer[URLSaveAsynchronous][urls[[i]], names[[i]], progFunction,      "Progress" -> True], {i, 10000}];      Here is the working Document.               CreateDocument[ExpressionCell[#, "Input"] & /@ cells]      Maybe can speed up it? Another aspect is how to do the whole job in Mathematica more comfortably. To build the download tasks took so much time in the working document which looks like I'm not downloading in an asynchronous way.


title: Gather lists according to another list

text: I have a first list of pairs like this one:               list1= {{1, 2}, {3, 4}, {5, 6}, {7, 8}, {9, 10}, {11, 12}, {13, 14}, {15, 16}, {17, 18}, {19, 20}, {21, 22}, {23, 24}, {25, 26}, {27, 28}, {29, 30}, {31, 32}, {33, 34}, {35, 36}, {37, 38}, {39, 40}, {41, 42}, {43, 44}, {45, 46}, {47, 48}, {49, 50}}      And a second one containing information about how the elements of the first list are related:               list2 = {{1, 46}, {5, 27}, {8, 42}, {21, 26, 29}, {31, 45}, {33, 36}, {35, 37}, {47, 49}}      How can I gather elements of the first list using list2 as a reference to get groups like these:               {{{1, 2}, {45, 46}, {31, 32}}, {{5, 6}, {27, 28}}, {{7, 8}, {41, 42}}, {{21, 22}, {25, 26}, {29, 30}}, {{31, 32}, {45, 46}}, {{33, 34}, {35, 36}, {37, 38}}, {{47, 48}, {49, 50}}}      So what it's happening here is that, element `{1, 46}` of list2 is telling that the pair containing element "1" (`{1, 2}`) of list1 is linked to element containing element "46" (`{45, 46}`) and so on. Note that when you use list 2 element `{35, 37}` to gather `{{35, 36}, {37, 38}}` you should consider that list 2 element `{33, 36}` already linked `{{33, 34}, {35, 36}}` so you will have a three element group `{{33, 34}, {35, 36}, {37, 38}}`. I tried to use graph related functions but couldn't make it work. Any ideas?


title: How to plot a decahedron with Mathematica?

text: I wrote               Graphics3D[{Opacity[.9], FaceForm[Yellow], PolyhedronData["Dipyramid", "Faces"]}]      and I get this picture: ![decahedron](http://i.stack.imgur.com/FQGu6.jpg) I would like to plot a decahedron without any defects on its surface. How can I do this? Thanks!


title: Draw an arbitrary convex polyhedron without excess diagonals drawn

text: I want to draw a set of convex polyhedrons whose vertices are defined by spherical coordinates on the surface of a unit sphere. Currently I followed the advice from here: http://mathematica.stackexchange.com/a/21842/651 But it draws diagonals on any tetragonal face. Each face is composed of triangles. Is it possible to accomplish the same but without tetragonal faces diagonaled? The current code:               Needs["TetGenLink`"]     mapping[{1,0,0}]={0,0,1}     mapping[{1,Pi,0}]={0,0,-1}     mapping:=CoordinateTransformData[{"Spherical"->"Cartesian"},"Mapping"]     verteces[n_]:=Flatten[Table [{1,Pi k/n, i 2Pi/Binomial[n,k]},{k,0,n},{i,0,Binomial[n,k]-1}],1]     Convex[n_]:=TetGenConvexHull[mapping/@ verteces[n]]     pts[n_]:=First[Convex[n]]     surface[n_] :=Last[Convex[n]]     b:=5     Graphics3D[{Yellow,Opacity[.9],GraphicsComplex[pts[b],Polygon[surface[b]]], Black, Line[{{0,0,-1.1},{0,0,1.1}}]}]


title: What is wrong with this Cobweb plot

text: Something is wrong with my cobweb plot code but I couldn't figure out where do the mistakes come from. So basically we have the iterative map $x_n == x_{n-1}^{2} - 3 \mu$. I would like to draw a cobweb diagram that shows in the range $-\frac{1}{12}<\mu<\frac{1}{4}$ there exists a fixed point $\frac{1}{2}(1-\sqrt{1+12\mu})$. Below are the code and the diagram: ![enter image description here](http://i.stack.imgur.com/vrHFf.png) ![enter image description here](http://i.stack.imgur.com/6mNGh.png) ![enter image description here](http://i.stack.imgur.com/MJM7P.png) My questions are:   1\. Why the diagonal line in the cobweb is not appearing and why is it replaced by a horizontal blue line?   2\. The label on top is wrong and I couldn't fix it.   3\. How can I fix the cobweb code so that I can show, say let $-\frac{1}{12}<\mu=0.2<\frac{1}{4}$ and initial value 2, there is a convergence towards a fixed point $\frac{1}{2}(1-\sqrt{1+12\mu})$. Many thanks. I really appreciate any helps.


title: "The function DiscretizeGraphics is not implemented for GraphicsComplex" confusion

text: Why does `DiscretizeGraphics` seems to work on one `GraphicsComplex` and not the other? Here is an example that works:               v = {{1, 0}, {0, 1}, {-1, 0}, {0, -1}};     p1 = Graphics[GraphicsComplex[v, Polygon[{1, 2, 3, 4}]]];     DiscretizeGraphics[p1]      But this does not               p2 = Graphics3D[First@ParametricPlot3D[{Cos[t],  Sin[u], c Sin[t]},         {u, 0, 2 Pi}, {t, 0, 2 Pi}]];     DiscretizeGraphics[p2]          (*The function DiscretizeGraphics is not implemented for \      GraphicsComplex[{{0.9999999999998993`,4.487989505128125`*^-7,*)      But `p2` is a GraphicsComplex? Looking at `FullForm[p2]` ![Mathematica graphics](http://i.stack.imgur.com/UmA56.png) Here is the `FullForm` for `p1` ![Mathematica graphics](http://i.stack.imgur.com/C8qcA.png) Are not `p1` and `p2` both `GraphicsComplex` ? `p1` is `2D` and `p2` is `3D`, but are they not both considered `GraphicsComplex`? ![Mathematica graphics](http://i.stack.imgur.com/mOlfR.png) It will good to know exactly what can and what can not be discretized. I tried to find this, but could not. All what I see are examples of usages so far. reference: http://www.wolfram.com/mathematica/new-in-10/data-and-mesh- regions/discretizing-graphics.html http://reference.wolfram.com/language/ref/DiscretizeGraphics.html?q=DiscretizeGraphics ![Mathematica graphics](http://i.stack.imgur.com/x61d6.png) I also looked at possible issues, and did not notice anything about this. Only one that came close is this `multiple volume primitives is not supported`. Is this the case here?


title: Expansion of $E(i c \mid m)$ at $c\to\infty$?

text: Currently, I am using a Windows machine with _Mathematica 8_. I noticed a difference in a series expansion of the function `EllipticE[]` in comparison with a result given by _Mathematica 9_ on Linux (which I was using previously). In _Mathematica 8_ on Windows the following input:               Series[EllipticE[I c, m], {c, Infinity, 0}] // PowerExpand // FullSimplify       produces a warning: > General::ivar: I c is not a valid variable. >> and the following output (slightly rearranged by me to better fit the browser): $$-\frac{i\left(6m+\frac{(1+m)}{\sinh^2(c)}\right)}{6\sinh(c)m^{3/2}}-\text{EllipticE}[m]+\frac{m \text{EllipticE}[\frac{1}{m}] + (1-m)\text{EllipticK}[\frac{1}{m}] }{\sqrt{m}}-$$ $$\frac{(1+m \cosh(2 c))\sqrt{\text{Limit}\big[-m\sinh^2(c)~,~i c \to 0\big]}}{\sinh^2(c)2m}$$ Same output as a code:               -((I (6 m Csch[c] + (1 + m) Csch[c]^3))/(6 m^(3/2))) -       EllipticE[m] +       (m EllipticE[1/m] - (-1 + m) EllipticK[1/m])/Sqrt[m] -       (((1 + m Cosh[2 c]) Csch[c]^2 Sqrt[Limit[-m Sinh[c]^2, I c -> 0]])/(2 m))      Now, in _Mathematica 9_ on Linux there was no such warning, and no `Limit` term appeared. I am confused about how to treat this `Limit` term, since it might just be a sign of something going terribly wrong in the guts of _Mathematica 8_. Does anyone have an advice on how to proceed? Maybe some of you can evaluate the same series expansion in a different version of _Mathematica_ so that we could compare results? ### EDIT Evaluating instead:               Series[EllipticE[c, m], {c, I Infinity, 0}] // PowerExpand // FullSimplify       worked as a charm without errors.


title: Randomly redrawing part of a random sample that matches a criterion

text: I would need some (a lot of) help to improve/rewrite some code. Since my real problem is a bit too complicated to explain, I try to show with a simplified scheme what I want to achieve. I create some random numbers on an interval between 0 and 1 with:               maximumRandomNumberFirstStep=10000;     randomNumberFirstStep=RandomReal[{0,1},maximumRandomNumberFirstStep];      Out of some reasons, I am only interested in values which are higher than a special limit value. So, I am searching for these values with:               limitValue=0.5;     helpListFirstStep=ParallelTable[If[randomNumberFirstStep[[i]]<limitValue,i,0],        {i,1,maximumRandomNumberFirstStep}];      For all values which are lower than this limit value I want to generate some new random numbers. At the moment I am doing this with:               maximumRandomNumberSecondStep=maximumRandomNumberFirstStep-Count[helpListFirstStep,0];     randomNumberSecondStep=RandomReal[{0,1},maximumRandomNumberSecondStep];      Now I want to add the new random numbers (randomNumberSecondStep) to all the old random numbers (randomNumberFirstStep) which are below my limit value. I am doing this with this very ugly code:               j=1;     randomNumberSecondStep2=Table[If[         helpListFirstStep[[i]]==0,         {randomNumberFirstStep[[i]],0},         {randomNumberFirstStep[[i]]+randomNumberSecondStep[[j]],j++}     ],{i,1,maximumRandomNumberFirstStep}];     randomNumberSecondStep2=randomNumberSecondStep2[[All,1]];      (By the way, this code does only work with `Table` but not with `ParallelTable` if I execute it more than once. Perhaps somebody has also a solution for that). If one creates a histogram `Histogram[randomNumberSecondStep2]`, it can be seen, that there are still values which are below my limit value. So, I have to copy my code from above, rename Second to Third and execute it again. I am doing this so often until none of my `randomNumberNthStep2`-values are below my limit value. This works, but I think it cannot be programed worse than I am doing it (copying, renaming etc. is really not nice, especially if your real problem is more complicated because of more variables). So, how can I write a better program in which all the copying and renaming steps can be avoided? I think there should be a solution in which all these steps are executed automatically. Unfortunately, my programing skills are too low to develop a nice solution. But I hope that some people here can give some help and hints. I would appreciate every support!


title: Alternatives to DateQ of Calendar Package

text: This answer notes that `Calendar` package has a function `DateQ` but the documentation suggests that the input is limited to a fixed format `YY,MM,DD` etc. Is there a more robust way to check any String for date-form ? (Self answer, below)


title: How to make binary image from black holes on white background

text: How to convert black holes on white background, so that I can see a binary image of holes on white image or bright dots on black image. ![enter image description here](http://i.stack.imgur.com/Eskg5.png) I can create a mask then also I am stuck with bright dot on black area, and then remaining white area.                ms = ColorNegate[image]      ![enter image description here](http://i.stack.imgur.com/TlMUh.png)


title: How can I run MUnit TestSuites outside WorkBench?

text: I would like to run MUnit tests inside TestSuite constructs from the command line or a notebook interface without loading the Workbench GUI. I tried               TestRun["suite.mt"]      where suite.mt contains:               TestSuite[{"file1.mt","file2.mt"}]      and each of the fileX.mt contain standard `Test[...]` constructs. However, this does not run the tests; I get the output that 0 tests were run. `TestRun["file1.mt"]` works. I cannot find documentation for `TestRun`, so maybe it does not support this, but in that case, how do I run testsuite files without loading the Workbench GUI?


title: Having Trouble Fitting Complicated Exponential Function with FindFit

text: I am working on high pressure shock experiments and am trying to figure out a the temperature of a sample under strong shocks. I have a measure value from 4 wavelengths of a pyrometer and am using a model to fit these points to a blackbody-like spectrum in order to back out a bulk true temperature, and a temperature of the hot spots/cracks that might also get viewed. The 4 temperature data points, namely:               {{1.8*10^-6, 1266.38}, {2.3*10^-6, 1137.87}, {3.5*10^-6, 901.019}, {4.8*10^-6, 751.279}}      with the second value of each pair corresponding to a real temperature. I am fitting them using the following hot spot model:               hotspotmodel = c2/(lam Log[(1/((1 - alpha)/(Exp[c2/(lam TT)] - 1) + alpha/(Exp[c2/(lam THS)] - 1))) + 1])      I am trying to use FindFit in te following way, but cannot get an answer that makes realistic sense:               pointsfit = FindFit[points, {hotspotmodel, TT > 0}, {{alpha, 0.05}, {TT, 500}, {THS, 2000}}, lam]          ""          function = Function[{lam}, Evaluate[hotspotmodel /. pointsfit]];     p1 = Plot[function[lam], {lam, 10^-7, 6 10^-6}, PlotStyle -> Thickness[.005]];     p2 = ListPlot[points2, AxesOrigin -> {0, 0},  PlotMarkers -> {Automatic, 12}];     Show[{p1, p2}, AxesOrigin -> {10^-7, 400}, PlotRange -> All,  AxesLabel -> {"\[Lambda (m)", "T (K)"}]      The problem is the function is not fitting properly, or returning sensible results. If I take away the constants, it returns a negative true temperature (TT), or a negative hot spot temperature (THS) and then has issues with trying to fit complex variables. My sense of the physical system says that true values should be within ~25% of those guesses. Any idea where I am going wrong? Thanks.


title: Notation for implicit derivative

text: Consider the following line of code:               D[x == y^3 + x y, x, NonConstants -> y]      The output would be:               1 == y + x D[y, x, NonConstants -> {y}] + 3 y^2 D[y, x, NonConstants -> {y}]      This is a confusing and cumbersome notation for the more natural:               1 == y + x y' + 3y^2 y'      I am trying to use the Notation package to help me replace the messy, default, output with the more natural one.   I have not succeeded. Of-course I have read this Q. The answer, if relevant, seems inaccessible to me, unfortunately.   Any help would be appreciated.


title: How to generate a real-time stream of data?

text: I'm experimenting with Inversion-of-Control (callback) programming patterns, and I would like a way to generate streams of data at specified real-time intervals. I came up with the illustrated way to use Animate (or Manipulate) to do it, and that's fine for experiments with the notebook front end, but I would like to have pure-code solutions I could put in a .m package. As always, I will be grateful for any hints, advice, and solutions! ![enter image description here](http://i.stack.imgur.com/dvVuZ.jpg)


title: Quick question about the FinancialBond function

text: If I want the fair price of a bond with coupon interval twice per year with interest rate $i$ p.a., I use "CouponInterval"->1/2. However, **the InterestRate should be the effective, or nominal interest rate?** Thank you for help.


title: How to generate a list from an iterator and some arbitrary function

text: Suppose I'd like to generate a table of integers:               Table[i, {i, 10}]      Great; now suppose I want only the integers that are even               Select[Table[i, {i, 10}], EvenQ]      This is fine for small things, but if I'd like to iterate over say $2^{n}$ things, only maybe $n$ of which will satisfy the predicate, it's a pretty terrible way to do it. Is there some good idiomatic way to do it? I could use a `Do` loop (I guess with `Append`ing?) but there must be a better way, right?


title: Eigenvector Anomaly

text: I'm trying to compute the eigenvectors for: $$ M = \left( \begin{array}{ccc} 1 & 4 \\\ 4 & 100 \end{array} \right) $$ Both myself and Mathematica report the **eigenvalues** as: $$ \lambda_1 = \frac{1}{2} \left(101+\sqrt{9865}\right) \approx 100.161 \\\ \lambda_2 = \frac{1}{2} \left(101-\sqrt{9865}\right) \approx 0.838647 $$ But when I ask for **eigenvectors** , the answer changes depending on whether I input the numbers as integers or floating-point numbers. When I ask for `Eigenvectors[{{1, 4}, {4, 100}}]` I get: $$ v_1 = \left(\frac{1}{8} \left(-99+\sqrt{9865}\right), 1\right) \approx (0.0403383, 1) \\\ v_2 = \left(\frac{1}{8} \left(-99-\sqrt{9865}\right), 1\right) \approx (-24.7903, 1) $$ When I ask for `Eigenvectors[{{1.0, 4.0}, {4.0, 100.0}}]` I get: $$v_1 \approx (0.0403055, 0.999187) \\\ v_2 \approx (-0.999187, 0.0403055)$$ When I calculate by hand, I get a solution which matches the first query. So, am I going crazy and overlooking some important maths (not unlikely, very tired...) or is this a bug? Possible related issues?:   * Difference In Eigenvectors   * eigenvector bug?


title: Generic Matrices

text: A multivariate Gaussian distribution for $k$ dimensions looks like this: $$\frac{1}{\sqrt{(2\pi)^k|\mathbf \Sigma|}}\mathbf e^{ -\frac 1 2 \mathbf x^T \mathbf \Sigma^{-1}\mathbf x }$$ $\mathbf x$ is a $k$-dimensional vector, and $\mathbf \Sigma$ a $k$ by $k$ symmetric positive definite matrix. Now I want Mathematica to integrate this expression over all $\mathbf x$, the result should be 1. How can I do this? Nothing else is known, in particular we don't know the number of dimensions.


title: Removing not-a-number entries from an array of numbers

text: I have a list of `{x, y, z}` coordinates from a wavefront sensor I'm working with. However, not all of the imported points have actual values, but they still get imported as `{x, y, "NaN"}`. I'm trying to find a way to get rid of these faulty points, but so far, I can't find a way to identify them. They don't affect some graphs, but I've run into trouble when performing operations on the set of points I've tried several things like `Delete`, `DeleteCases` and `StringTrim`, but to no avail. When I run, `DeleteCases[graph, "NaN"]`, for example, I get back the same set of points. Also, I am using _Mathematica_ 7\. Thanks.


title: Nonlinear System of ODEs

text: I want to solve this NonLinear System: $$\Theta_{\xi} = \frac{(1-\lambda)K_{\alpha}}{\beta}\Phi^{2}\Theta+\frac{2\beta-1}{6\beta}\xi\Theta$$ $$\Phi_{\xi} = \frac{1}{2\lambda K_{\alpha}\ \Theta^{3}\Phi+2\beta}\bigg[-3\lambda K_{\alpha}\ \Phi^{2}\Theta^{2}\Theta_{\xi}+\xi\Phi\bigg]$$ where I would like to use $\beta$ as a parameter, whereas $K_\alpha$ and $\lambda$ are known constants. The boundary conditions I have are: $$\Theta[0]=0,\ \Theta[10]=0$$ $$\Phi'[0]=0,\ \Phi'[10]=0$$ I tried to use NDSolve, but I get this error:                \[Lambda] = 0.1; K = 1;          eqn1 = ((1 - \[Lambda]) K)/\[Beta] \[CapitalPhi][x]^2 \[CapitalTheta][x] + (2 \[Beta] - 1)/6 x \[CapitalTheta][x];           eqn2 = \[CapitalPhi]'[x] == (-3 \[Lambda] K \[CapitalPhi][x]^2 \[CapitalTheta][x]^2 \[CapitalTheta]'[x] + x \[CapitalPhi][x] )/(2 \[Lambda] K \[CapitalPhi][x] \[CapitalTheta][x]^3 + 2 \[Beta]);               {xsoln, ysoln} = {\[CapitalTheta], \[CapitalPhi]} /. NDSolve[{\[CapitalPhi]'[x] == eqn2, \[CapitalTheta]'[x] == eqn1, \[CapitalTheta][0] == 0, \[CapitalPhi][0] =       0}, {\[CapitalTheta], \[CapitalPhi]}, {x, 0, 10}]          NDSolve::deqn: Equation or list of equations expected instead of 0 in the first argument{(\[CapitalPhi]^\[Prime])[x]==((\[CapitalPhi]^\[Prime])[x]==(x \[CapitalPhi][x]-0.3 <2>> <<1>>)/(2 \[Beta]+0.2 Power[<<2>>] \[CapitalPhi][<<1>>])),(\[CapitalTheta]^\[Prime])[x]==1/6 <<3>>+<<1>>,<<1>>,0}. >>      etc


title: Potential hill simulation

text: I attempted to simulate the motion of a particle approaching to a potential hill. The particle should experience a force of `-1/r^2`, if the equation for the potential hill is `1/r`. So the acceleration for the particle should be `-c/r^2`, where c is an arbitrary proportional constant. The particle is expected to bounce back when close enough to the hill, but it passed through the hill instead, despite of different values of c. Is there anything wrong with my code or the equations of motions. Secondly, when `PlotPoints->10` changed to `PlotPoints->30` in order to obtain an intact potential hill, the program becomes unbearable slow. How to get around this?               Clear[x, y, z, t, ux, uy, bx, by, cnt];     t = -2;     ux = 1; uy = 1;     bx = -1; by = -1;     obx = -1; oby = -1;     cnt = 0.2;     ax[x_, y_] := -cnt*x/( x^2 + y^2)^(3/2)  ;     ay[x_, y_] := -cnt*y/( x^2 + y^2)^(3/2)   ;          field = ContourPlot3D[        1/(x^2 + y^2)^(1/2) - z == 0, {x, -2, 2}, {y, -2, 2}, {z, 0, 50},         PlotPoints -> 10];          particle[t_] := (        bx = ux t + 0.5 ax[obx, oby] t^2;        by = uy t + 0.5 ay[obx, oby] t^2;        bz = 1/(bx^2 + by^2)^(1/2) ;        obx = bx; oby = by;        Graphics3D[{          Red, Sphere[{bx, by, bz}, 0.1]}]);          Dynamic[      Refresh[       Row[{Show[field, particle[t], ImageSize -> {300, 300}],          If[t < 2, t = t + 0.1, t = -2]}]       , UpdateInterval -> 1       ]        ]      ![enter image description here](http://i.stack.imgur.com/OceMq.gif)


title: Problem importing individual frames from QuickTime file

text: I am having difficulties with a code I am trying to run. Essentially I want to `Read` in a large `.mov` video file (60000+frames). Then take each frame individually, process it (e.g. crop and identify particles, I can do these steps) and `Export` a text file with processing data (can do), plus the cropped image. I am not sure which functions/the setup to use for reading in a video and exporting an image sequence. I have looked at `Import` with `.mov` files, and can identify the number of frames in the movie:               videof = Import[path <> file <> ".mov", "FrameCount"]      as well as view a single frame               frame1 = Import[path <> file <> ".mov", {"Frames", {n}}]      however there is something wrong in the second line of code; no matter the value of `n` I always get the first image. Advice is much appreciated.


title: Simultaneous equations

text: I am trying to solve two simultaneous equations that are claimed to have 145 solutions but I believe they have 143, I want to check I am right before sending an email to my professor. Mathematica is taking an incredible amount of time with this so I thought my code was wrong. Could you please check it for me?               Length[NSolve[{x == Sin[6 Pi y], y == Sin[6 Pi x]}, {x, y}]]      The equations are $y=sin(6 Pi x)$ and $x=sin(6 Pi y)$ apologies if my conjecture is wrong!


title: Generating an Array of Vectors

text: I'm working a lot with nested lists of $(x,y)$ coordinates (of a FEM-like mesh, but that's not important), like so:               nodes = { {{0,0}, {0,1}, {0,2}}, {{1,0}, {1,1}, {1,2}} }      This has the nice convenient structure that `nodes[[i,j]]` gives me the $(x,y)$ coordinates of my "$(i,j)$th" node, so that `nodes[[i,j]][[1]]` is my $x$ coordinate and `nodes[[i,j]][[2]]` is my $y$ coordinate. This works nicely for my purposes. How would I construct a _variable_ size array of this kind using `array[]` et. al? If I use for instance               nodes = array[array[X,2],{4,4}]      I get a 4-by-4 array where `nodes[[i,j]]` is of the form `{X[1],X[2]}[1,1]`. This is close to what I want, and I could change a bunch of my code to work with this, but I'd rather not - I'd love to have something where `nodes[[i,j]]` looks like this:               {X[i,j][1], X[i,j][2]}      The reason I'd like to do this is I have defined some functions that take an array of nodes as input and computes some quantity, i.e. `F[nodes_]:= (*something*)`, and I'd like to compute derivatives of this function with respect to each of its variables, i.e. $F_{x_i}$ or $F_{y_j}$. I've tried various things with no success, and I've also tried looking for a similar construct in the documentation and here on SE, but I didn't see anything. Any ideas? Edit: I had too many square braces in some of my code before. Fixed.


This question follows on from my earlier one, Seeking strategies to deploy a function securely without a front end. As it pursues a different track to solve the same problem it seemed distinct enough and perhaps interesting enough to others to merit its own question. I need to deploy a single function into a Windows based technology environment so that the users can call the function from some of their process workflows. The users don't need a front end or user interface of any kind. In fact a front end would cause problems. I have some additional code security concerns, but approaches to that have begun to emerge (see the earlier question referenced above). Running without a front end would also keep code from prying eyes. **Question:** Can one launch a Player Pro kernel independently of the front end via something like a Mathematica Script to run the contents of the computable document without displaying the document itself? References or examples of how to do this appreciated.


title: Solve non linear system of equations

text: I am trying to solve a system of two equations with two unknowns. In these equations I have, a part from constants:   * Unknown nr 1, $$D_{\perp}$$   * Unknown nr 2, $$\omega_C$$   * Known function of r: $$\mu(r)$$ The full system looks like: equation1: $$ D_{||}= \frac 1 {2 \omega_0}(-\alpha-1)\sqrt{(-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2}+ \frac{\alpha}{2\omega_0}\sqrt{(-2\omega_0-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2}+ \frac{1}{2\omega_0}\sqrt{(2\omega_0-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2} $$ and equation2: $$ 1= \mu(r) \frac{\alpha}{2\omega_0} \ln{\left[\frac{-2\omega_0-\omega_C+\mu(r)D_{||}+\sqrt{(-2\omega_0-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2}}{-\omega_C+\mu(r)D_{||}+\sqrt{(-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2}}\right]}+ \frac{\mu(r)}{2\omega_0} \ln{\left[\frac{2\omega_0-\omega_C+\mu(r)D_{||}+\sqrt{(2\omega_0-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2}} {-\omega_C+\mu(r)D_{||}+\sqrt{(-\omega_C+\mu(r)D_{||})^2+(\mu(r)D_{\perp})^2}}\right]} $$ So the solution of the system of equations will be $$D_{\perp}(r),\;\omega_C(r)$$. What I've tried to do is simply                   Solve[{equation1, equation2},{Dorthogonal, omegaC}]      but Mathematica keeps on running forever without any output. I have also tried:                   DorthFun[r_]:=Solve[{equation1, equation2},{Dorthogonal, omegaC}][[1,1]]         omegaCFun[r_]:=Solve[{equation1, equation2},{Dorthogonal, omegaC}][[1,2]]      and it just keeps on running... It doesn't return any errors. Just...eternal running. Forrest Gump Syndrome... I have also tried to solve the system putting $$\mu(r)=1$$ without any change. I have given Mathematica about 20 minutes. Should I give it more time or does this mean that Mathematica cannot solve this? Or is there something I could do differently? Thank you for your help! My code looks like this:                   NSolve[{         Dparallel==1/(2 omega0) Sqrt[(-omegaC + Dparallel)^2 + (Dorth)^2] (-alpha-1)+alpha/(2omega0)Sqrt[(-2omega0-omegaC+Dparallel)^2+(Dorth)^2]+1/(2omega0)Sqrt[(2omega0-omegaC+Dparallel)^2+(Dorth)^2],         1 == alpha/(2omega0)Log[(-2omega0-omegaC+Dparallel+Sqrt[(-2omega0-omegaC+Dparallel)^2+(Dorth)^2])/(-omegaC+Dparallel+Sqrt[(-omegaC+Dparallel)^2+(Dorth)^2])]+1/(2omega0)Log[(2omega0-omegaC+Dparallel+Sqrt[(2omega0-omegaC+Dparallel)^2+(Dorth)^2])/(-omegaC+Dparallel+Sqrt[(-omegaC+Dparallel)^2+Sqrt[(-omegaC+Dparallel)^2+(Dorth)^2]])]         }, {Dorth, omegaC}]


title: Multiplying SI Prefixes

text: I'm doing a few calculations and would like to keep my units straight. I figured that using the built-in Units functionality could help with this. However, using SI prefixes doesn't seem to give the desired output. I try executing this code:                In[15]:= A = 4 Milli Meter*Pi*3^2 Micro Meter // N      And instead of getting, say, 113E-9 Meter^2, I get:                Out[15]= 113.097 Meter^2 Micro Milli      Milli, Meter, and Micro are all globally defined names. What am I doing wrong?


title: Numerical method for solving a polynomial equation

text: I am looking for a numerical solution of a equation which contains, in general, one polynomial equation with unknown variable x. I have tried `Reduce`, `Solve`, `NSolve` and `FindRoot`, but I'm getting nothing after running my code for a few days. Probably there is solution in some of numerical iterative or other methods, unknown for me.                a1 = Sqrt[2*10^11*3.055*10^-5];      a2 = 0.01;      a3 = 7849;      a4 = 10;      a5 = 50;      a6 = 11;      a7 = 700;           k1 = Rationalize[      Solve[(k (k^3 a1 - I a2 - k (a4 + a5^2 a3)) + x - 2 k a5 a3 a6 -            a3 a6^2) == 0, k][[4]][[1]][[2]], 0];      k2 = Rationalize[      Solve[(k (k^3 a1 - I a2 - k (a4 + a5^2 a3)) + x - 2 k a5 a3 a6 -            a3 a6^2) == 0, k][[1]][[1]][[2]], 0];      k3 = Rationalize[      Solve[(k (k^3 a1 - I a2 - k (a4 + a5^2 a3)) + x - 2 k a5 a3 a6 -            a3 a6^2) == 0, k][[3]][[1]][[2]], 0];      k4 = Rationalize[      Solve[(k (k^3 a1 - I a2 - k (a4 + a5^2 a3)) + x - 2 k a5 a3 a6 -            a3 a6^2) == 0, k][[2]][[1]][[2]], 0];           expr = (I  a1 (k1 - k3) (k2 - k3) (k1 - k4) (k2 - k4))/(      k1 + k2 - k3 - k4) - a7 a6^2;           NSolve[expr == 0, x]


title: Problem with CurrentValue["MouseOver"] when graphics are wrapped with Deploy

text: Please take a look at following toy example:               Graphics[{      Dynamic[{       If[CurrentValue["MouseOver"], Red, Blue], Thickness@.05, Line[{{0, 0}, {.5, 0}}]}],      Dynamic[{       If[CurrentValue["MouseOver"], Red, Blue], Disk[{.5, 0}, .2]}]       }, PlotRange -> {{0, .6}, {-.3, .3}}]      > ![enter image description here](http://i.stack.imgur.com/U6aIC.png) You can put cursor over the `Disk` or `Line`. As you can see the `Disk` covers the `Line`, so you can only get: > ![enter image description here](http://i.stack.imgur.com/1jioK.png) or > ![enter image description here](http://i.stack.imgur.com/wUYbr.png) But as soon as we add `Deploy` (`Deploy@Graphics[...`), the `Disk` becomes transparent to the `"MouseOver"` event, and we can get the following result (with cursor on the right place): > ![enter image description here](http://i.stack.imgur.com/Zfhfl.png) This behaviour is so unexpected, I'm tempted to call it a bug. This change in behaviour is huge! In this context, `Deploy` as a quick and nice function seems to be unreliable. :/ V9.0.1 Win7 * * * The walkaround is similar like there: Delay of events detection dependent of “MouseDown” position. But I doubt the reason is the same.


title: Plot sets in the complex plane

text: I just want to know if there's a way to plot sets in the complex plane. For example $$A=\\{z\in \mathbb{C},z+e^{z}=0\\},\\\A=\\{z\in \mathbb{C},\Re(z)+e^{z}\geq0\\}.$$


title: Creating a package of special plots

text: In my work as a math teacher, I create lots of graphs to illustrate ideas. I have lots of math books and will often try to duplicate a particular style of graph. Below, one very simple example, created with Mathematica. (but as a "one off") ![enter image description here](http://i.stack.imgur.com/22GIY.png) I wish my Mathematica skills were better, I spend too much time grading and lesson planning, and I know I could be more effective with lesson writing if my toolkit was bigger and more robust. So I decided to try to write some functions to create these kinds of plots, and thought it would be best to put them into a package so I could use them easily. I have a program that created the plot, and it works fine, but I'm embarrassed to share it, it's so clunky...I don't care about speed and concise code, I wanted to be able to understand all the steps and pieces...               specialA[xmin_, xmax_, xscale_, ymin_, ymax_, yscale_, xtickmin_,        xtickmax_, xtickinc_, ytickmin_ , ytickmax_, ytickinc_, funclist_,        stylelist_, ratio_] :=       Module [{xvals, lightx, darkx, gridx1, gridx2, xgrid, yvals, lighty,         darky, gridy1, gridy2, ygrid, xticks, yticks},            xvals = Range[xmin, xmax, xscale];       lightx = Select[xvals, Not[IntegerQ[#/5]] &];       darkx = Select[xvals, IntegerQ[#/5] &];       gridx1 = {#, GrayLevel[0.8]} & /@ lightx;       gridx2 = {#, Black} & /@ darkx;       xgrid = Flatten[{gridx1, gridx2}, 1];       yvals = Range[ymin, ymax, yscale];       lighty = Select[yvals, Not[IntegerQ[#/5]] &];       darky = Select[yvals, IntegerQ[#/5] &];       gridy1 = {#, GrayLevel[0.8]} & /@ lighty;       gridy2 = {#, Black} & /@ darky;       ygrid = Flatten[{gridy1, gridy2}, 1];       xticks = Range[xtickmin, xtickmax, xtickinc];       yticks =  Range[ytickmin, ytickmax, ytickinc];            Plot[funclist, {x, xmin, xmax},         PlotStyle -> stylelist,        GridLines -> {xgrid, ygrid},         PlotRange -> {{xmin, xmax}, {ymin, ymax}},         AspectRatio -> ratio,         Ticks -> {xticks, yticks},         ImageSize -> 430,         Epilog -> Text[0, {0, 0}, {1, 1}],        BaseStyle -> {FontSize -> 18, FontFamily -> "Times"}]       ]      So, for example....                   f[x_] := 2 x^2 - x - 8;     g[x_] := -x^2 - 2 x + 2     specialA [-7, 7, 1, -32, 32, 2, -5, 5, 5, -30, 30, 10, {f[x], g[x],        f[x] - g[x], g[x] - f[x]}, {Red, Blue, Green, Magenta}, 2]      I'd welcome any comments on that and welcome any nice plotting routines for "high school math" graphs that others might be using. But here is my real question. I'm no expert but have created a few simple Mathematica packages. I put this function into a package, wrote a usage statement, defined the function. When I use it, I get a grid, but the functions don't display. I am sure this has something to do with contexts , but I just have no clue. How can I improve my function definition, and how would I include it in a working package of similar functions to draw other graphs? (I save my package online in case that would be helpful... https://dl.dropbox.com/u/24411154/grapher.nb) I know this is a huge question, I do appreciate the feedback and help from this group. Tom


title: How to integrate ParametricNDSolve solution with respect to a parameter?

text: I have just upgraded to the new version of Mathematica because of its new built-in `ParametricNDSolve` function. I need to solve a first-order non- linear ordinary differential equation that has two parameters. I have no problem obtaining appropriate solutions to this equation, but once I have the solution I need to integrate out the dependence on the parameters. For instance, I would like to write something like               sol=ParametricNDSolveValue[{-y*a'[x]==(a[x]^2-1)f[x]+a[x],a[-10]==0.6},a,{x,-10,10},{y}]     Integrate[sol[y][x],{y,0,1}]      If I have an `InterpolationFunction` which has two arguments I can very easily integrate over one of the two arguments, but this does not seem possible to do with a `ParametricFunction`. The best idea I have had so far is to form an `InterpolationFunction` from the `ParametricFunction` and then integrate from there, but that seems inefficient and I would like something a little faster and more straightforward. Any help would be greatly appreciated.


title: Equivalent of Python's "all" function in Mathematica

text: The Python function                def isPrime(n):         return all(n % i for i in xrange(2, n))      checks if a number is a prime number by using `all`. How can I write a function similar to `all` in _Mathematica_?


title: How get Singular.m package to interact with Singular itself?

text: I've installed the package `Singular.m` by Kauers and Levandovskyy from http://www.risc.uni-linz.ac.at/research/combinat/software/Singular/ as well as the Singular CAS itself. How, exactly does one get the package to communicate with the Singular system -- specifically to get results back from Singular to Mathematica? In a terminal, I just type `singular` and that starts the system. So at the end of the package `.m` file I modified the final line so it now reads `SingularCommand = "singular"`. Now after I load the package and evaluate one of the interface functions, such as one in the `demo.nb` notebook, I get an error: `Singular::error : Unable to read Singular's output.`


title: How to save a web page in mathematica?

text: For a given web page http://www.nytimes.com/, is it possible to save it as : > **first** : only HTML(within css) > > **second** : more elements(include pictures and so on) just like what a browser can do. I have tried this, but it only generates a large picture with plain text.               Export["F:\\nytimes.html", Import["http://www.nytimes.com"], "HTML"]


Is Mathematica a Turing-complete language? If so, how can that be proved? If not, why?


title: Define a function that returns a function with a variable number of arguments

text: Previously, I've defined a function with a variable number of arguments with the following method.               x = {x1, x2, x3};     f1 = Function[       Evaluate[x],       x + 1 //       Evaluate     ];     f1 @@ x     Output: {1 + x1, 1 + x2, 1 + x3}      Also, I've defined a function that returned a function as follows.               g[x_] := Function[       {z},       x + z //       Evaluate     ];     f2 = g @@ {z};     f2 @@ {x1}     Output: x1 + z      But when I attempt to put these ideas together, I have an issue. The hang-up seems to be in passing the argument of the first function to the argument of the second. Note that in our second example, above, we weren't passing `x` as an argument to the inner function. However, this is exactly what we want to do:               h[x_] := Function[       Evaluate[x],       x + 1 //       Evaluate     ];     f3 = h @@ x;     f3 @@ x     Output: h[x1, x2, x3][x1, x2, x3]      The desired output is `{1 + x1, 1 + x2, 1 + x3}`. My motivation for writing such a function is to bury it in a package.


title: I've some data points in Mathematica. I want to plot these in MATLAB. How to do that?

text: I want to export these data points into matlab and get a plot. Like I've to plot ft vs t. Through the following command- Table[{[Tau], ft[2, 5, [Tau]]}, {[Tau], 0, 15, .01}], i generated some data set. Now how do I export it to MATLAB to plot?


title: How to convert a sparsearray to a graph?

text: s = Import["hello.mat"];     t = GraphEmbedding[Graph[s[[1]]], "GridEmbedding"];      I get the following error: > GraphEmbedding::graph: A graph object is expected at position 1 in >   >   >     GraphEmbedding[Graph[SparseArray[<1232>, {223, 223}]], GridEmbedding]. >   I would prefer to do SpringElectricalEmbedding and then GridEmbedding. What should I do? Please help.


title: How can I get GridEmbedding for my graph?

text: I evaluate the following               s = Import["hello.mat"];     t = GraphEmbedding[s[[1]], "GridEmbedding"];      and I get the these error messages > GraphPlot::mthd: The value of option Method -> GridEmbedding should be > Automatic, "SpringElectricalEmbedding", "SpringEmbedding", "LayeredDrawing", > "LayeredDigraphDrawing", "RadialDrawing", "HighDimensionalEmbedding", > "CircularEmbedding", "SpiralEmbedding", "LinearEmbedding", or > "RandomEmbedding". > > GraphPlot::grph: SparseArray[<1252>, {240, 240}] is not a valid graph. How can I do `GridEmbedding` for my sparse array `s`?


title: How can I plot a circular region and lines?

text: I am having difficulty plotting something along these lines in _Mathematica_. ![Imgur](http://i.imgur.com/l7CTZIy.png) I will explain what is going on here. Basically, I am solving an interval equation $ax=b$, in the first example $a = [1, 20]$ and $b =[3, 4]$ and therefore $b/a = [3/20, 4]$. As you can see from the graph I want a circle representing the projective line. A straight line will come off the two end points of the interval which represents the solutions and those lines will go through the center of the circle and the region bounded by the lines should be shaded in. Similarly the other example below is an unbounded one where $a = [-1,2]$ and $b = [5,100]$ and hence $b/a = (-\infty, -5] \cup [5/2, \infty]$. Same type of plot though, the lines come off $-5$ and $5/2$ on the x axis and go through the center of the circle and the bounded regions are shaded. Could someone please help me plot this?


title: Using Mathematica to make a graphics array of circles

text: I need to make .svg file that contains an array of black circles that all have radii of 5 units ad are evenly spaced apart from each other 32 units. The difficulty in creating this design is the sheer number of circles necessary. I will need a 3175 by 3175 array. At first I tried doing this manually in Adobe Illustrator, but then it occurred to me that I could write a code to do it for me in Mathematica.               n = 50;     mat = Flatten[Table[{32 i, 32 j}, {i, 0, n - 1, 1}, {j, 0, n - 1, 1}], 1];     g = Disk[#, 5] & /@ mat;     Export["circle_array.svg", Graphics[g], "SVG"];      Above is the code I wrote, and it for a small number of rows/columns, it works just fine. After it's created, I open it up in Adobe Illustrator, and it looks just fine. However, when I try to increase the number of rows much greater than 100 X 100, after it's exported and I open it up in Adobe Illustrator. I notice that not all the circles are evenly spaced away from each other, and I don't understand why. If anyone could help, it would be greatly appreciated.


title: Function generating a rule with a Patterns: e.g. toRule[ lhs, rhs] into lhs[x_] -> rhs

text: I want to create a rule as a function of `z` from a left and right hand side. My (non-working) code and usage is:               toRule[l_, r_] := l[x_] -> r          toRule[f, a+b]      I want the output to be:               (* => f[x_] -> a + b*)      However: when I do this with the code above I get:               (* => f[x$_] -> a + b*)      Where is the $ coming from, and how can I get rid of it?


title: Restarting Mathematica causes an animation to fail to run

text: I've tried animating two different sets of 3D plots. Interestingly, one plot will always work, regardless of whether I restart Mathematica, yet another plot will only work in its current Mathematica session. For instance, I have:               h = 1;     α = 1;     H = (1/2)*(p^2 + x^2);     rho[n_] = (1/Pi)*(-1)^n*Exp[-2*H/h]*LaguerreL[n, 4 H/h];     list = ParallelTable[Plot3D[rho[n], {x, -10, 10}, {p, -10, 10}, PlotRange -> All,     PlotLabel -> Style[n, FontSize -> 40], ImageSize -> 800,      AxesLabel -> {Style[x, FontSize -> 30],        Style[p, FontSize -> 30]}], {n, 0, 10}];     ListAnimate[list, AnimationRate -> 1, AnimationRunning -> False]      This yields: ![enter image description here](http://i.stack.imgur.com/99THp.png): If I save my file and restart Mathematica, while none of my variables are initialize, the animation is still able to run. However, I have another piece of code to make a different animation:               Clear[n]     α = 0.15;     h = 1;     η = 1;     a = -n - Sqrt[1 + 16 r];     r = η^2/(4*α^2);     qrn = (1/2)*(n + 1/2) + (1/4)*Sqrt[1 + 16 r];     c = 1 - 2*qrn;          Clear[x, p]     ψ[p_] = (1/(1 + α*p^2)^qrn)*  Hypergeometric2F1[a, -n, c, (1/2) +    (I/2)*Sqrt[α]*p]          For[i = 0, i <= 2, i++, Clear[list2]]     For[n = 0, n <= 2, n++,      list2[n] = Append[list2[n],      ParallelTable[{x = RandomReal[{-5, 5}], p = RandomReal[{-10, 10}],        NIntegrate[        Exp[I*u*x/h]*         Conjugate[ψ[p + u/2]]*ψ[p - u/2], {u, -10000,          10000}, MaxRecursion -> 12]}, {40}]]] // AbsoluteTiming     Animate[ListPlot3D[list2[n][[2]], PlotRange -> All, ImageSize -> 400], {n, 0, 2, 1}, AnimationRunning -> False]      This gives: ![enter image description here](http://i.stack.imgur.com/IWGNy.png) After saving, and restarting Mathematica, I now get: ![enter image description here](http://i.stack.imgur.com/fN0Ge.png) Creating my animation with 40 points only takes a couple seconds, so it isn't terrible that I'd have to recreate it. However my plan is to increase the number of points to more than 12000 and increase `n` to 10, which will take much more time than a couple seconds. Why does one animation run after re-opening Mathematica, yet another animation will give me an error message? As a related question, is there a way for Mathematica to always remember the values of my variables, even after restarting Mathematica? Thanks.


title: How do I find the value of specific printing options?

text: I'd like to programmatically save (and then temporarily change and restore) specific printing options. I get as far as               Options[Notebook, PrintingOptions]      but don't know what to do next to extract specific options, for example `PrintingMargins`, from the (long) resulting list.


title: Dataset and string conversion, using Internal`StringToDouble

text: i'm trying to get familiar with the new Dataset-stuff in Mathematica 10. As it will (hopefully) provide some benefit for me, I have started to convert my app to move from list and list of rules to Associations and Datasets. Starting point for me is an XML object, which is converted to a Dataset now. I was hoping to reduce the programming code by e.g. applying functions to multiple columns at once. This is an example Dataset:               testset = Dataset[{<|"Element" -> "A", "a" -> "529251", "b" -> "520358"|>, <|     "Element" -> "B", "a" -> "51"|>, <|"Element" -> "C",      "a" -> "177"|>, <|"Element" -> "L", "a" -> "125"|>, <|     "Element" -> "S", "a" -> "1343"|>}]      This gives: ![dataset](http://i.stack.imgur.com/QOSPH.png) As the XML sometimes contains integers and reals in scientific form, I'm using               Internal`StringToDouble      to convert the stuff:               fnStr2Real[x_,r_:0]:=If[Head[x]===Missing,r,Internal`StringToDouble[x]]     SetAttributes[fnStr2Real, {Listable}];     testset[fnStr2Real, {"a", "b"}]     (* Trace yields *)     (* Query[fnStr2Real, {"a", "b"}][testset] *)      The issue: there is some kind of rounding in the resulting recordset: ![result](http://i.stack.imgur.com/2Wvkh.png) If I apply my function directly to my list or association, everything is fine:               Function[x, fnStr2Real[#[x] & /@ Normal[testset]]] /@ {"a", "b"}     (* {{529251., 51., 177., 125., 1343.}, {520358., 0, 0, 0, 0}} *)      I know that Internal`StringToDouble is an internal function, but I found several posts in this forum using it, and I really like it! Is there an issue with M10, Dataset and Internal`StringToDouble? Are there any other robust string coversion functions availabe which can handle integers, reals and reals in scientific form **inside a dataset**? Thanks!


title: How to combine the sum of two rational polynomials into a single rational polynomial

text: I have browsed the web and found no information. I have two rational polynomials $f1(x,y)$ and $f2(x,y)$. When I perform `f1 + f2`, I always get two separate functions. However, I would like to combine them into one rational expression. How do I do that?


title: Save state of buttons, sliders, listpicker in enterpriseCDF?

text: How do you make enterpriseCDFs save the state of buttons sliders, listpickers etc? I manipulated a CDF and saved it on exit but all the buttons etc were reset when I started it again. I've made a simplified example: If you execute the last line it produces a toggler bar and a button in a single output cell. I export that cell as a stand alone enterprise CDF. When I open it in cdf player, toggle a few buttons, save, close and reopen, the state is reset to the original first three buttons pressed. Is it possible to get the state of the buttons to be remembered across saves?               data[n_] := Range[n];     series = {1, 2, 3};     cdf =        DynamicModule[{},          Grid[{           {TogglerBar[Dynamic[series], Range[6]]},           {Dynamic[              PopupWindow[Tooltip[Button["Show data"], data[series]], data[series],                WindowSize -> All,                WindowTitle -> "data",                 WindowElements -> {"VerticalScrollBar", "HorizontalScrollBar","StatusArea"},                ScrollingOptions ->                   {"HorizontalScrollRange" -> Automatic, "VerticalScrollRange" -> Automatic}],               TrackedSymbols :> Full]}},           Alignment -> {Top, Top}],           SaveDefinitions -> True];     cdf      I also tried exporting with the command:               CDFDeploy[NotebookDirectory[] <> "testCDF.cdf", cdf,        "Method" -> "Standalone", "Target" -> "PlayerPro"];      with same result (not sure if target needs to be playerpro though; the aim was for it to be an enterpriseCDF state-savable with the standard player.


title: Using function with multiple definitions in Manipulate

text: I'm trying to use `Manipulate` to visually try out different values of lambda in a Box-Cox transformation. I've created a `boxcox` function with two definitions to deal with both the normal case and the case when lambda is 0:               boxcox[data_, 0] := Log[data]     boxcox[data_, l_] := (data^l - 1)/l      Then I use this function inside `Manipulate` but I keep getting tons of errors. It looks like Manipulate is only using the general definition and starts complaining about dividing by zero.               Manipulate[      pdata = Partition[boxcox[data, u], 12];      ranges = Max[#] - Min[#] & /@ pdata;      means = Mean[#] & /@ pdata;      mrdata = Transpose[{means, ranges}];      mrlm = LinearModelFit[mrdata, x, x];      Show[       ListPlot[mrdata, Axes -> False, Frame -> True,         AxesOrigin -> {Automatic, 0}],       Plot[mrlm[x], {x, Min[means], Max[means]}]       ],      {u, 0.00, 1.00}      ]      Here is the data I'm using in case it matters:               data = {154., 96., 73., 49., 36., 59., 95., 169., 210., 278., 298., 245., \     200., 118., 90., 79., 78., 91., 167., 169., 289., 347., 375., 203., \     223., 104., 107., 85., 75., 99., 135., 211., 335., 460., 488., 326., \     346., 261., 224., 141., 148., 145., 223., 272., 445., 560., 612., \     467., 518., 404., 300., 210., 196., 186., 247., 343., 464., 680., \     711., 610., 613., 392., 273., 322., 189., 257., 324., 404., 677., \     858., 895., 664., 628., 308., 324., 248., 272.}


title: Text replacement rules in $FrontEnd?

text: I want to set up a part of my Mathematica package so that typing `<``-``Space` is replaced with the left arrow character (`←` or `\[LeftArrow]`) just like it does with `-``>``Space`. Does anyone know how to do this? I looked briefly at event handling but didn't see anything immediately useful. I already know, of course, that `Esc``<``-``Esc` will insert `←`, and that doing so is not that much additional effort, but I'm more interested in the concept involved here. Any ideas?


title: Some way to identify the source of a notebook file?

text: Is there a way to tell where a Mathematica file came from, meaning some type of machine ID (yes, I know about `$MachineID`). I want to be able to identify students who are sharing files, i.e. one student creating the file and giving it to others. Anti-cheating is what I'm trying to accomplish.


title: Detect each polygons in an image

text: Considering this image: ![](http://i.imgur.com/Y87duSz.jpg) Is there a way to colorize each tiles independently? My problem is that               img = Import["http://i.imgur.com/Y87duSz.jpg"];     Colorize@img      gives ![Mathematica graphics](http://i.stack.imgur.com/ENWxM.png) where polygons with very light colours are merged. Note that `EdgeDetect[img, 1, .05]` produces: ![Mathematica graphics](http://i.stack.imgur.com/v4q0C.png) if it can help. * * * **Ideally** , the idea would be to colour the tiles independently (`Colorize` colours them in terms of their "blackness"). That would imply that one detects the polygons and recreate them, which is harder than "only" colouring them in terms of their pre-existing colours. Please note I have absolutely no knowledge in image-processing.


title: Can the Notation function support complex structures on its left hand side?

text: In my investigation of the possibilities offered by the Notation package, I came over a limitation in the complexity of the left hand side symbol. For instance, although I am able to enter                Notation[        ParsedBoxWrapper[          OverBar[SuperStar[af]]]        \[DoubleLongLeftRightArrow]        ParsedBoxWrapper["afstarbar"]]      without Mathematica complaining, **evaluating** `afstarbar` **returns an incorrectly formatted symbol** (the whole symbol is replaced by a red error square instead of displaying $\overline{af^*}$). The output is:                OverBar[SuperStar[$CellContext`af]]           An unknown box name (OverBar) was sent as the BoxForm for the expression.      Check the format rules for the expression.      I am unfamiliar with the concept of "boxing", but I tried the following:               Notation[       ParsedBoxWrapper[         MakeBoxes[           OverBar[SuperStar[af]]]]       \[DoubleLongLeftRightArrow]       ParsedBoxWrapper["afstarbar"]]      And the display is slightly better, as the $af^*$ is here, but the bar is replaced by the red error square and the output is:               \!\(\*OverscriptBox[     SuperscriptBox[\(af\), \(*\)],      Notation`Private`singleBlank[]]\)          An unknown box name (Notation`Private`singleBlank) was sent as the BoxForm      for the expression. Check the format rules for the expression.      Any trick on how to get the symbol to display correctly with the Notation package?


title: Dynamically updating interactive objects in Mathematica

text: Is it possible to create interactive code in Mathematica (via EventHandler[] or Manipulate[] perhaps) that continually updates some calculation based on external data in the time when no user input is received? As a trivial example, code that shows the current time, with seconds, continually updating, and which allows the user to click on any of the numbers to change their color. Thank you.


title: Winkel Tripel and Robinson Projections?

text: How can I do WorldPlot with the WorldProjection -> Winkel Tripel or Robinson projections. Mathematica has Mollweide and several others defined, but not these popular projections adapted by National Geographic.


title: How do I emulate MATLAB's comet plot?

text: I have tried the following code:               ParametricPlot[{Sin[12 u] Cos[u], Sin[12 u] Sin[u]}, {u, 0, 2 Pi}]      and it produced the following image: ![Flower](http://i.stack.imgur.com/0ZgdU.png) MATLAB has a `comet` command that will draw this parametric plot in real time, allowing the viewer to watch as the particle traces out the plot. Is there a way to do this in _Mathematica_?


title: How to find derivative of a numerical solution, where precision is ambiguous?

text: I am trying to take the derivative of a numerical solution. I am concerned that the way I'm doing this may be problematic due to numerical error; I think there must be a better way but I'm not very experienced in Mathematica. The relevant function is called GammaFcn and is defined as follows:               K = 10; c = 1; a = 1; b = 3; d = .01;  tradeoff = 1; mig = .002;     N1[x_] := K*(1 - c*x + b*x^a - mig);     f[x_, y_, t_] := 1/(1 + N1[x]*E^(c*t*(y - x)));     Gstar[x_, y_, t_] := s*(y*f[x, y, t] + x*(1 - f[x, y, t]));     sStar[x_, y_, t_] :=      s /. NSolve[(1 - s/K) + b*(Gstar[x, y, t]/s)^a -        x*c*(1 - f[x, y, t]) - y*c*f[x, y, t] - mig == 0, s][[1]]     N2[x_, y_, t_] := f[x, y, t]*sStar[x, y, t];     Beta1[x_] := mig*N1[x];     Estar[x_] := (d)/(Beta1[x]);     Istar[x_] := 1 - Estar[x];     lambda[x_] := Beta1[x]*Istar[x];     I1star[SigmaI_, x_] :=     Estar[x]*(1 - Estar[x])*mig*N1[x]/(d + SigmaI*lambda[x]);     M1[SigmaI_, p_, q_, t_, z_] := Sum[z^(n)*(d + SigmaI*Istar[p]*mig*N1[p])* E^(-t*(d + SigmaI*Istar[p]*mig*N1[p]))* E^(-N1[q]*mig*t)*(N1[q]*mig)^n*t^n/n!, {n, 0, Infinity}]     M2[p_, q_, t_, z_] := Sum[z^(n)*((d)*E^(-t*(d))* E^(-mig*(sStar[q, p, t] - N2[q, p, t])*t)*(mig*(sStar[q, p, t] - N2[q, p, t]))^n*t^n/n!), {n, 0, Infinity}];     M3[p_, q_, t_, z_] := Sum[z^(n)*((d)*E^(-t*(d))* E^(-mig*N2[p, q, t]*t)*(mig*N2[p, q, t])^n*t^n/n!), {n, 0, Infinity}];     Part1[SigmaI_, p_, q_, z_] := NIntegrate[M1[SigmaI, p, q, t, z], {t, 0, Infinity}]     Part2[p_, q_, z_] := NIntegrate[M2[p, q, t, z], {t, 0, Infinity}]     Part3[p_, q_, z_] := NIntegrate[M3[p, q, t, z], {t, 0, Infinity}]     Const1[SigmaI_, Res_, New_] := (d/(d + SigmaI*Istar[Res]*mig*N1[Res]));     Correct[SigmaI_, Res_, New_, Z_] := Const1[SigmaI, Res, New]*Part1[SigmaI, Res, New, Z] + (1 - Const1[SigmaI, Res, New])*(Part1[SigmaI, Res, New, Z]*      Part2[Res, New, Z])          GammaFcn[SigmaI_, p_, q_, z_] := (1 - Estar[p] - I1star[SigmaI, p]) + Estar[p]*Correct[SigmaI, p, q, z] + I1star[SigmaI, p]*(SigmaI*Part3[p, q, z] + (1 - SigmaI));      I am interested in solutions to GammaFcn[SigmaI, p, q, z]=z, where 0>z>1; suppose these solutions are given by a function GammaFcnFixedPt[SigmaI, p, q]. In particular, I want to calculate the first derivative of GammaFcnFixedPt with respect to p at the point p=q, with some fixed SigmaI and q. I do this with the following function, for some small epsilon (because when p=q, GammaFcn=1):               GammaFcnFixedPtDeriv[SigmaI_, q_, epsilon_] := (-1 + z /. FindRoot[GammaFcn[SigmaI, q + epsilon, q, z] - z, {z, .99}])/  epsilon      My problem is that when I do this, when epsilon gets sufficiently small, the results diverge. So with decreasing epsilons, the function will appear to be converging onto the derivative, but then (presumably due to precision errors), it diverges. Here are pairings of epsilon and the approximated derivative (from above function), for SigmaI=.1, q=.3. : {1/10, -0.200505}, {1/10^2, -1.44905}, {1/10^3, -0.901373}, {1/10^4, -0.886587}, {1/10^5, -0.886699}, {1/10^6, -0.886358}, {1/10^7, -0.847299}, {1/10^8, -0.297589}, {1/10^9, 0.969691} So in the above pairings, before epsilon goes below 1/10^7, there it looks as though the sequence is approaching somewhere around .88, but thereafter diverges quickly. Again I assume this is due to some machine error; but I'm not positive it is even due to this error. But if this is the case, how do I know what epsilon to pick for the most accurate derivative? If I add "PrecisionGoal -> 10" to the integrals in the code in the first block above:               Part2[p_, q_, z_] := NIntegrate[M2[p, q, t, z], {t, 0, Infinity}, PrecisionGoal -> 10]     Part3[p_, q_, z_] := NIntegrate[M3[p, q, t, z], {t, 0, Infinity}, PrecisionGoal -> 10]      {1/10^1, -0.200505}, {1/10^2, -1.44905}, {1/10^3, -0.901373}, {1/10^4, -0.886587}, {1/10^5, -0.886699}, {1/10^6, -0.886431}, {1/10^7, -0.855601}, {1/10^8, -0.0572474}, {1/10^9, -13.9658} And doing the same thing with PrecisionGoal->12 gives: {1/10^1, -0.200505}, {1/10^2, -1.44905}, {1/10^3, -0.901373}, {1/10^4, -0.886587}, {1/10^5, -0.886699}, {1/10^6, -0.886431}, {1/10^7, -0.855601}, {1/10^8, -1.83026}, {1/10^9, -13.9658} When I set it much higher than that it takes a long time. So how can I know that the derivative is converging to (in this case) approximately -0.887, rather than convergence not occurring (limit not existing)? And is there a better way to find the derivative than I have done? Thanks very much.


title: Running .NET/Link on Mono (On MacOS or any UNIX or Linux Distro)

text: Has anyone have any luck running .NET/Link on Mono? When I run .NET/Link to call Mathematica from within a Mono program, this is the error message I receive.               Unhandled Exception:     System.TypeInitializationException: An exception was thrown by the type initializer for Wolfram.NETLink.Internal.NativeLink ---> System.DllNotFoundException: mathlink       at (wrapper managed-to-native) Wolfram.NETLink.Internal.MacMathLinkAPIProvider:MLBegin (intptr)       at Wolfram.NETLink.Internal.MacMathLinkAPIProvider.extMLBegin (IntPtr zero) [0x00000] in <filename unknown>:0        at Wolfram.NETLink.Internal.NativeLink..cctor () [0x00000] in <filename unknown>:0        --- End of inner exception stack trace ---       at Wolfram.NETLink.MathLinkFactory.createMathLink0 (System.String cmdLine, System.String[] argv) [0x00000] in <filename unknown>:0        at Wolfram.NETLink.MathLinkFactory.CreateMathLink (System.String cmdLine) [0x00000] in <filename unknown>:0        at Wolfram.NETLink.MathLinkFactory.createKernelLink0 (System.String cmdLine, System.String[] argv) [0x00000] in <filename unknown>:0        at Wolfram.NETLink.MathLinkFactory.CreateKernelLink (System.String cmdLine) [0x00000] in <filename unknown>:0       Now the C# code for calling Mathematica is:                   _kl = MathLinkFactory.CreateKernelLink("-linkmode launch -linkname '\"/Applications/Mathematica.app/Contents/MacOS/MathKernel\" -mathlink'");      The same argument worked for J/Link, but for some reason does not work on NETLink.


title: How do I create a matrix of functions?

text: I'm trying to create a matrix in which the elements are functions of two variables, but I can't figure out how to do it. Is it possible? The only way I could figure out is to define the matrix in terms of two variables, then replace those variables with values later on. Here's my code:               row1 = {};     row2 = {};     For[ii = 1, ii <= Length[modes], ii++,       {l, m} = modes[[ii]];       G1lm = sglm[-2, l, m, \[Theta]] + (-1)^l sglm[-2, l, -m, \[Theta]];       G2lm = sglm[-2, l, m, \[Theta]] - (-1)^l sglm[-2, l, -m, \[Theta]];       AppendTo[row1, G1lm Cos[m \[Phi]]];       AppendTo[row1, G1lm Sin[m \[Phi]]];       AppendTo[row2, G2lm Sin[m \[Phi]]];       AppendTo[row2, -G2lm Cos[m \[Phi]]];       ];     proj = {row1, row2};      Theta and Phi are my variables. Later on, when I want proj for a specific Theta and Phi I use replace:               pN = proj /. {\[Theta] -> a, \[Phi] -> b};      where a and b are some numbers. The problem with this solution is that if I use Theta and Phi for something else later on, say I set Theta = 0.5, then proj will permanently have Theta = 0.5, regardless of the replace I use. This happened to me: I was calling proj in a loop, but I kept getting the same results. After ripping my hair out for a few hours, I realized I had Theta set to a specific value in a block of code further up in my notebook. The simple solution is to simply not use Theta and Phi anywhere else in the notebook, but that doesn't seem very stable. The root problem here is I'm using global variables to do something that should really be local to proj, but I don't know how to fix that. The ideal would be if I could somehow make proj a function of Theta and Phi so that later on I can simply do `pN = proj[a,b]`. I know how to make the Gnlm's functions of theta, but then I don't know how to make row1, row2, and by extension proj, functions of Theta (and Phi). Any help would be greatly appreciated.


title: How to use Compile with NDSolve function for large no of differential equations?

text: Matrix `M[t]` and and Hamiltonian of same order `L^n x L^n`.               L=3; n=3;tmax = 2;     Rho[t_] = Table[M[i,j][t], {i,L^n}, {j,L^n}]; //Timing          Rho0=ParallelTable[M[i,j][0],{i,L^n},{j,L^n}];     Hop[t_, i_] =     Transpose[       Table[If[i == j, 1, 1+Cos[t]], {i, {1}}, {j, L}]].Table[If[i == j, 1, 1+Sin[t]^2], {i, {2}}, {j, L}];          I1[i_Integer] := I1[i]=IdentityMatrix[L^(i-1)];     I2[i_Integer] := I2[i]=IdentityMatrix[L^(n-i)];          H[t_] = Sum[FixedPoint[ArrayFlatten, I1[i] \[TensorProduct] Hop[t,i]      \[TensorProduct]I2[i]],{i,n}];      I have made a code for large no. of time dependent differential equations with initial values. I want to compile this for using NDSolve. Here is my code...               Rho1[t_] =     Table[         D[Rho[t],t][[i,j]] == H[t][[i,j]],         {i,L^n}, {j,L^n}]; // Timing          Fun[t_]= -I (#1.#2 - #2.#1)&[H[t],Rho[t]];     Eq0[i_,j_]:=If[i == 1 && j == 1,#==1.,#==0. ]&[Rho0[[i,j]]];     Eq0a=ParallelTable[Eq0[i,j],{i,L^n},{j,L^n}];     Eqsopwo[t_]=Join[Thread /@ Thread[D[Rho[t], t] == Fun[t]],Eq0a]//Flatten;     \[Rho]Termsop[t_]=Rho[t]//Flatten;          sol= NDSolve[Eqsopwo[t],Diagonal[Rho[t]],{t,0,tmax},DependentVariables->\[Rho]Termsop[t],MaxSteps->10^5,MaxStepFraction-> 1/50, Method-> {"EquationSimplification"->"Solve"}]//Chop;     Plot[Evaluate[{Re[Diagonal[Rho[t]]]}/.sol],{t,0,tmax},PlotRange->All]


title: No Symbols Matching $InterfaceEnvironment

text: I'm running a webMathematica server (running Mathematica 9) and am getting some weird stuff in the log - specifically, I keep getting the following message:               General pool request background 61 > [      ] Input MSP`Utility`DestroyRequest[] generated message:     Remove::rmnsm: There are no symbols matching "Global`$InterfaceEnvironment".      I can't find any documentation for this, so it's possible that this is a user created symbol that I'm not seeing, but I think this is unlikely. Is this error something I should worry about?


title: Proper treatment of roots and powers in Series?

text: I have the following problem in _Mathematica_ 9 on Linux. I let _Mathematica_ compute the `Series` expansion:               Block[{$Assumptions = {r, p}\[Element] Reals},            Series[Log[r^2 - Sqrt[r^4 + p^2]], {p, 0, 0}]     ]      But the output is a rather disappointing naked infinity:               Log[r^2 - Sqrt[r^4]] + O[p]^1      Similarly, for:               Block[{$Assumptions = {p, r} \[Element] Reals},            Series[1/(r^2 - Sqrt[r^4 + p^2]), {p, 0, 0}]     ]      I get:               1/(r^2 - Sqrt[r^4]) + O[p]^1      Why is _Mathematica_ handing me naked infinities like that, even though I specified that the computation is to be carried out in the real numbers? What do I do wrong? How can I avoid this? EDIT: To make my point more clear, I would expect the Series to give results like $$\log(c_1 p^2)+O(p^1)$$ and $$\frac{1}{c_2 p^2}+O(p^1)$$ as a regulated infinity, instead of $\log(0)$ and $\frac{1}{0}$.


title: A Case of DeleteCases for WeatherData

text: I am trying to use `ListVectorPlot` in combination with `WeatherData` and `CountryData` to show wind velocity vectors over a given country. The problem I am having stems from "NotAvailable" and "NotApplicable" data from a function in which the wind's velocity vector is constructed. The function is               WindFoo[x_, y_] := -WeatherData[{y, x}, "WindSpeed"] {1/Cos[y °],         1} Through[{Sin, Cos}[WeatherData[{y, x}, "WindDirection"] °]]      Note: I am not able to display the degree symbol, so it is written out in the above code. How can the code above be modified to exclude the cases of "NotAvailable" and "NotApplicable" data to work in manner similar to the wind velocity vectors shown over Australia in the following tutorial: http://blog.wolfram.com/2009/01/13/visualizing-weather-patterns-in- mathematica-7/


title: How to control the order of meshing？

text: I want the first pic，but i got the second. What can I do? ![enter image description here](http://i.stack.imgur.com/HS45U.png) ![enter image description here](http://i.stack.imgur.com/NJKfa.png) Below if the code to produce the second pic:               data = {{14.808253637272422`,          5.992415735398138`, -10}, {14.937895181494618`,          5.835230243884557`, -10}, {15.104476099263962`,          6.190504320078969`, -10}, {15.549905577158725`,          6.242337037739479`, -10}, {16.214292480395272`,          5.946484831244263`, -10}, {17.172819193670374`,          6.044620763590807`, -10}, {14.808281657120743`,          5.992346489742819`, -5}, {14.937944420202394`,          5.83506166050832`, -5}, {15.104458590180693`,          6.1905995697624006`, -5}, {15.549909194291196`,          6.242289007293565`, -5}, {16.21428558314019`,          5.946304785847829`, -5}, {17.172828128632702`,          6.044675503473975`, -5}, {14.808370059168467`, 5.992128021118139`,          0}, {14.93796801813446`, 5.834980865961923`,          0}, {15.104521660293425`, 6.190256467374961`,          0}, {15.549917959537805`, 6.242172617112528`,          0}, {16.214298214579777`, 5.946634515939104`,          0}, {17.172827172286173`, 6.044669644434473`,          0}, {14.808196500526428`, 5.992556937854838`,          5}, {14.937952476296239`, 5.835034078072448`,          5}, {15.104459789437353`, 6.190593045787643`,          5}, {15.549890996009355`, 6.242530654951938`,          5}, {16.21428524490538`, 5.946295956593386`,          5}, {17.17285709312181`, 6.044852953886377`,          5}, {14.808354331354467`, 5.992166889379481`,          10}, {14.937875586457652`, 5.835297333327767`,          10}, {15.1044365297471`, 6.190719578861889`,          10}, {15.549897447086584`, 6.242444993694446`,          10}, {16.214291019501033`, 5.946446696177589`,          10}, {17.172815114360944`, 6.044595771776242`, 10}};          pica = Graphics3D[Sphere[data, 0.1]]          picb = ListPlot3D[data, InterpolationOrder -> 1,        MeshFunctions -> {#2 &, #3 &}, Mesh -> All, PlotRange -> All]          Show[pica, picb]


title: Using Collect to gather explicitly imaginary terms

text: Collect[Expand[(x + I y)^6], I]      yields               x^6 + 6 I x^5 y - 15 x^4 y^2 - 20 I x^3 y^3 + 15 x^2 y^4 + 6 I x y^5 - y^6      How can I obtain               x^6 - 15 x^4 y^2 + 15 x^2 y^4 - y^6 + (6 x^5 y - 20 x^3 y^3 + 6 x y^5) I


title: Creation of one column matrix

text: I want to create an one column matrix with the 1st element to be -1, and all the other elements to be zero. However, I want to be able to make a matrix of this kind with arbitrary dimension, that is, an arbitrary number of rows.


title: How to find a max point from discrete data of function?

text: I have a lot of points like $\\{\\{x_1,y_1,z_1\\},\\{x_2,y_2,z_2\\},...,\\{x_n,y_n,z_n\\}\\}$ which obtained by a function $z=f[x,y]$, and I want to find the max point of $z$, say $z_i$, and return the point coordinents $\\{x_i,y_i,z_i\\}$. How do I do? I have checked the help, but find there is no function can return the point coordinents from a list.


title: How to drop the identical rows and columns of a matrix?

text: I'm currently trying to find a way to eliminate from a square matrix of generic dimension the duplicate (up to sign) rows and columns, obtaining a smaller dimensional matrix. I searched around here but I couldn't find an answer, neither I managed to combine methods for list manipulation in order to arrive at my result. I'll make an example in order to clarify. Suppose I have a matrix of this form: $$ \left(\begin{array}{cccccc} x_1 & w_2 & -x_1 & -w_1 & -w_2 & w_1\\\ z_2 & x_2 & -z_2 & z_1 & -x_2 & -z_1 \\\ -x_1 & -w_2 & x_1 & w_1 & w_2 & -w_1\\\ -y_2 & y_1 & y_2 & x_3 & -y_1 & -x_3\\\ z_2 & x_2 & -z_2 & z_1 & -x_2 & -z_1\\\ y_2 & -y_1 & -y_2 & -x_3 & y_1 & x_3\\\ \end{array} \right) $$ I want to eliminate rows and columns 3,5,6, (copies of 1,2,4) ending up with the matrix $$ \left(\begin{array}{ccc} x_1 & w_2 & -w_1\\\ z_2 & x_2 & z_1\\\ -y_2 & y_1 & x_3\\\ \end{array} \right) $$ I would like to implement this operation so that it can act on matrices of arbitrary dimension. Thank you very much! Stefano


title: Can anyone help to solve this set of coupled ODE?

text: After some time Mathematica just give me the same thing while I was trying to solve: DSolve[{3 E^(-[Lambda]2 [Phi]0) r a[r] - 1/2 [Lambda]1 [Phi][r] + 3/8 E^(-[Lambda]2 [Phi]0) [Lambda]1 [Phi][r] - 1/2 [Lambda]2 [Phi][r] + 3/2 E^(-[Lambda]2 [Phi]0) r^2 Derivative[1][a][r] - 2 Derivative[1][b][r] + 3/2 E^(-[Lambda]2 [Phi]0) Derivative[1][b][r] - r [Lambda]2 Derivative[1][[Phi]][r] + r^2 !(*SuperscriptBox[(b), * TagBox[ RowBox[{"(", "3", ")"}], Derivative], MultilineFunction->None])[r] == 0, -((6 b[r])/r^2) - (3 [Lambda]1 [Phi][r])/(8 r) + ( E^([Lambda]2 [Phi]0) [Lambda]1 [Phi][r])/(2 r) + ( E^([Lambda]2 [Phi]0) [Lambda]2 [Phi][r])/(2 r) - 21/2 r Derivative[1][a][r] + (9 Derivative[1][b][r])/(2 r) + ( 2 E^([Lambda]2 [Phi]0) Derivative[1][b][r])/r + E^([Lambda]2 [Phi]0) [Lambda]2 Derivative[1][[Phi]][r] - 9/2 r^2 (a^[Prime][Prime])[r] - 3 (b^[Prime][Prime])[r] - E^([Lambda]2 [Phi]0) r !(*SuperscriptBox[(b), * TagBox[ RowBox[{"(", "3", ")"}], Derivative], MultilineFunction->None])[r] == 0, -((2 [Lambda]1 a[r])/r) + ( 8 E^([Lambda]2 [Phi]0) [Lambda]1 a[r])/(3 r) - ( 8 E^([Lambda]2 [Phi]0) [Lambda]2 a[r])/(3 r) - ( 4 [Lambda]1 b[r])/r^3 + ( 16 E^([Lambda]2 [Phi]0) [Lambda]1 b[r])/(3 r^3) - ( 16 E^([Lambda]2 [Phi]0) [Lambda]2 b[r])/( 3 r^3) - ([Lambda]1^2 [Phi][r])/r^2 + ( 4 E^([Lambda]2 [Phi]0) [Lambda]1^2 [Phi][r])/(3 r^2) + ( 4 E^([Lambda]2 [Phi]0) [Lambda]2^2 [Phi][r])/(3 r^2) + ( 16 E^([Lambda]2 [Phi]0) [Lambda]2 Derivative[1][b][r])/( 3 r^2) - ( 8 E^([Lambda]2 [Phi]0) [Lambda]2 (b^[Prime][Prime])[r])/( 3 r) + ([Phi]^[Prime][Prime])[r] == 0}, {a[r], b[r], [Phi][r]}, r] or you can use: $DSolve[{3 E^(-\\[Lambda]2 \\[Phi]0) r a[r] - 1/2 \\[Lambda]1 \\[Phi][r] + 3/8 E^(-\\[Lambda]2 \\[Phi]0) \\[Lambda]1 \\[Phi][r] - 1/2 \\[Lambda]2 \\[Phi][r] + 3/2 E^(-\\[Lambda]2 \\[Phi]0) r^2 Derivative[1][a][r] - 2 Derivative[1][b][r] + 3/2 E^(-\\[Lambda]2 \\[Phi]0) Derivative[1][b][r] - r \\[Lambda]2 Derivative[1][\\[Phi]][r] + r^2 \\!\\(\\*SuperscriptBox[\\(b\\), \\* TagBox[ RowBox[{"(", "3", ")"}], Derivative], MultilineFunction->None]\\)[r] == 0, -((6 b[r])/r^2) - (3 \\[Lambda]1 \\[Phi][r])/(8 r) + ( E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]1 \\[Phi][r])/(2 r) + ( E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2 \\[Phi][r])/(2 r) - 21/2 r Derivative[1][a][r] + (9 Derivative[1][b][r])/(2 r) + ( 2 E^(\\[Lambda]2 \\[Phi]0) Derivative[1][b][r])/r + E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2 Derivative[1][\\[Phi]][r] - 9/2 r^2 (a^\\[Prime]\\[Prime])[r] - 3 (b^\\[Prime]\\[Prime])[r] - E^(\\[Lambda]2 \\[Phi]0) r \\!\\(\\*SuperscriptBox[\\(b\\), \\* TagBox[ RowBox[{"(", "3", ")"}], Derivative], MultilineFunction->None]\\)[r] == 0, -((2 \\[Lambda]1 a[r])/r) + ( 8 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]1 a[r])/(3 r) - ( 8 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2 a[r])/(3 r) - ( 4 \\[Lambda]1 b[r])/r^3 + ( 16 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]1 b[r])/(3 r^3) - ( 16 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2 b[r])/( 3 r^3) - (\\[Lambda]1^2 \\[Phi][r])/r^2 + ( 4 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]1^2 \\[Phi][r])/(3 r^2) + ( 4 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2^2 \\[Phi][r])/(3 r^2) + ( 16 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2 Derivative[1][b][r])/( 3 r^2) - ( 8 E^(\\[Lambda]2 \\[Phi]0) \\[Lambda]2 (b^\\[Prime]\\[Prime])[r])/( 3 r) + (\\[Phi]^\\[Prime]\\[Prime])[r] == 0}, {a[r], b[r], \\[Phi][r]}, r]$


title: how to define roots of an equation as a function?

text: I need to solve an equation depending on 3variables in only one variables and define the roots as a function of the remaining variables? I used some thing like that:               root[x_,y_]:=Solve[x^2+y^2+z^2==0,z][1];      but it doesn't ?


title: Testing if a specific expression is a solution for ordinary differential equation

text: I would like to know what is the full command details of wanting to test if a specific expression is a solution for ordinary differential equation (without using DSolve). Thanking you in advance for your assistance.


title: how can I make a control for my vector equation in my simulation?

text: I'm trying to build a control selector or menu for equations to put into my curvature simulations. However all the code I've tried hasn't worked.Does anyone know how I could make a control selector or menu for a position vector r[t] in my simulation?               r[t_] := {3 Cos[t], Sin[t]}     uT[t_] := Simplify[r'[t]/Norm[r'[t]], t \[Element] Reals]     vN[t_] := Simplify[uT'[t]/Norm[uT'[t]], t \[Element] Reals]     x[t_] := r[t][[1]]     y[t_] := r[t][[2]]          y2[t_] = r[t] + ((1/(uT[t]/r'[t])) vN[t])          OlcusionPoint[t_] := y2[t][[1]]     (*-((((x'[t])^2+(y'[t])^2)*y'[t])/((x'[t]*y''[t])-(x''[t]*y'[t])))*)     OclusionPoint2[t_] := y2[t][[2]]          (*+((((x'[t])^2+(y'[t])^2)*x'[t])/((x'[t]*y''[t])-(x''[t]*y'[t])))*)     Manipulate[      Show[Graphics[{RGBColor[0.6, 0.6, 0.8],          Disk[{OlcusionPoint[t], OclusionPoint2[t]},           Abs[1/(uT[t]/r'[t])]]}],        ParametricPlot[{x[t], y[t]}, {t, 0, 2 Pi}],        Graphics[{PointSize[.025], Point[{x[t], y[t]}]}]], {t, 0.0001,        2 Pi}]


title: Nested function problem

text: I have trouble defining the following function in _Mathematica_               sp[M_, u_] := 1/Sum[u^k M!/(M - k)!, {k, 0, M}]      After defining it _Mathematica_ displays an error. And in fact it prints a pair of values, but my application requires only one, i.e 2 argument function. As the result I need to plot this function in 3D with `u` changing from `0.01` to `1`, and `M` changing `1` to `100`. Any suggestions?


title: marking the points of x-intercept on a plot

text: OK, I have a nice plot of some oscillating motion. I can read some stuff off the graph, but I thought it would be cool to show a dot where the curve hits the x-axis and a little label showing the coordinates, and thus get a direct measure of the period of oscillation. (I know there are equations for that but that's not quite the point of the exercise). Anyhow, here's the code I am using               Clear[x1, x2, v1, v2];Clear[x1, x2, v1, v2];     ti = 0; tf = 4.; delta = (1/40);     x1[ti] = -0.04; x2[ti] = -0.04; v1[ti] = 0.0; v2[ti] = 0.0;     m1 = 0.1890; m2 = 0.1898; k = 3.65; kc = 3.5;     Do[ F1[t] = -k*x1[t] - kc*(x1[t] - x2[t]);       F2[t] = -k*x2[t] - kc*(x2[t] - x1[t]);      v1[t + delta] = v1[t] + (F1[t]/m1)*delta;       v2[t + delta] = v2[t] + (F2[t]/m2)*delta;     x1[t + delta] = x1[t] + v1[t + delta]*delta;          x2[t + delta] = x2[t] + v2[t + delta]*delta,      {t, ti, tf, delta}]      x1data = Table[{t, x1[t]}, {t, ti, tf, delta}];     x1plot = ListPlot[x1data, AxesLabel -> {"t", "x1"}, Joined -> True]      and with it I get a nice graph. I tried using NSolve and didn't get any output at all. So I think NSolve must not be geared to this. Understand I am a rank beginner with Mathematica. So it's ok to assume I am stupid :-) And I know that there are likely more efficient ways to do the curve and such; this is a lab exercise tho. We're testing a simulation against the real data. thanks!


title: Real and Imaginary parts of solutions to a complex linear ODE system

text: Consider a complex linear ODE system $x'=Ax$, where $$A=\left( \begin{matrix} 0&1\\\ -2&-i \end{matrix}\right). $$ One can first find the eigenvalues and eigenvectors using `JordanDecomposition` and get the general solutions with `MatrixExp`. Here is my **question** :   Is there a quick way in _Mathematica_ that I can get the solution in the form such that the real and imaginary parts are separated: $$ v(t)=c_1v_1(t)+c_2iv_2(t) $$ where $v_1(t),v_2(t)\in{\mathbb R}$?


title: Programming paradigm change

text: I'm new to Mathematica, I've been learning by myself, and, being a Java developer for some time, it makes it hard to adapt. I've tried to solve a Reddit challenge, and I came up with this solution:               str = "[can {and it(it (mix) up)} look silly]"     total = ""     str = StringTake[str, {2, -2}]     inix[strn_] :=        Min[StringPosition[strn, "["], StringPosition[strn, "{"],        StringPosition[strn, "("]];     finix[strn_] :=        Max[StringPosition[strn, "]"], StringPosition[strn, "}"],        StringPosition[strn, ")"]];         howmany[strn_] :=        Length[StringPosition[strn, "["]] +        Length[StringPosition[strn, "{"]] +        Length[StringPosition[strn, "("]];     brackets[x_, y_] :=        If[x != \[Infinity], (total =           StringTake[str, {1, x - 1}] <> " " <>           StringTake[str, {y + 1, -1}] <> " " <> total;           str = StringTake[str, {x + 1, y - 1}]) ,           (total = str <> " " <> total; str = "")];          Do[brackets[inix[str], finix[str]], {i, howmany[str] + 1}];     StringSplit[total, Whitespace]     StringJoin[Riffle[%, " "]]      Is there a better solution? From seeing some Mathematica code, this looks to be terribly big code. I would want just some pointers on how to improve.


title: How to set the NDSolve method to LSODA

text: * I notice that off all the `Method` options available for `NDSolve[...]`, `LSODA` is invoked quite automatically for stiff equations when `Method -> Automatic`. Is there a way to set the `Method -> LSODA`? I ask because as `LSODA` switches between Adams' and Gear's method, it should have sub-options to it that I could set (like order etc.)?   * Under `tutorial/NDSolveStiffnessSwitching` I encountered `Needs["DifferentialEquations`NDSolveProblems"]`. What other example problems are available and where can I find documentation on them? * * * Edit (03/24/2011): In light of comments made to this post, I'd also refer readers to page 409 of "The Mathematica Guidebook for Numerics". Springer has this book online and can be found here.


title: Positioning Vertices in a Graph

text: I made a graph that has a given number at each vertex: ![Lattice Graph for D_16](http://i.stack.imgur.com/NWURD.png) I want to place vertex 1 above vertex 2 if value of vertex 1 is larger than of vertex 2. And also, vertices with same value should be located on the same horizontal line. Above graph satisfies the condition. (I think because of luck.) But, another graph constructed by the same function does not satisfy it: ![Lattice Graph for D_10](http://i.stack.imgur.com/ZeKeU.png) Here 5 is not above 2. How can I implement this condition? (The height of the gap between 2 and 5 does not need to be exactly 3 times the height gap between 1 and 2.) ==== Editted ==== Here is the code for function makes the graph. (I simplified it.)               LatticeGraph[group_] :=         Module[{subgroup = (*omitted*), nodes, edges, list},             nodes = Table[Property[Labeled[i, Placed[subgroup[[i]], Center]],(*VertexStyle Omitted*)], {i, 1, Length[subgroup]}];             list = {(*Omitted and it's used for generate edges*)}             edges = {(*Omitted and it was just the list of v1->v2*)};                  Graph[nodes, edges, VertexSize -> Large,EdgeShapeFunction -> GraphElementData["Line"], ImagePadding -> 20]]


title: Density matrix - Mathematica

text: I would to write a density matrix "matr" related to a quantum system in a mixed state. "Ham" is the Hamiltonian of N-spin system. The code I wrote is:               {evl, eve} = Quiet[Eigensystem[Ham]];     ord = Ordering[evl];     evl = evl[[ord]];     eve = eve[[ord]];     matr = ConstantArray[0, {2^(Nspin/2), 2^(Nspin/2)}];     Monitor[Do[re = matrix[eve[[1]], a2, a1, Nspin];       matr[[a2, a1]] = re; matr[[a1, a2]] = re, {a1, 1, 2^(Nspin/2)}, {a2,         1, a1}], {a1, a2}]      The time needed to solve the problem for a large number of Nspin it's quite long. I woul like to know if is possible to modify the code in order to make the process faster.


title: Parallelization of REvaluate[]?

text: Is there a way to run several RLink computations via REvaluate[] in parallel, e.g. with ParallelTable[]?


title: Create a Matrix out of two matrices

text: I have two matrices `A = {{a1,a2},{b1,b2}}` and `B = {{x1,x2},{y1,y2}}`. I would like to create out of these matrices a new matrix `C={{{a1,x1},{a2,x2}},{{b1,y1},{b2,y2}}}` which I would like at the end transform to this form `{{a1,x1},{a2,x2},{b1,y1},{b2,y2}}`. Of course, I am looking for a solution to the general case. Thanks.


title: Plot hermite spline

text: I'm new to Mathematica and cannot figure out a way to plot a 2d hermite spline. My points are: $p_0=\\{143.05;195.606\\};$ $p_1=\\{159.373;181.845\\};$ $p_2=\\{160.891;179.26\\};$ $p_3=\\{176.423;171.878\\};$ I need to plot the spline from $p_0$ to $p_1$ I calculate Catmull-Rom tangents in the following way: $\left(t_1=0.5 \left(p_2-p_0\right);\right)$ $\left(t_2=0.5 \left(p_3-p_1\right);\right)$ Calculating the weights: $$ \left(h_1(\text{t$\\_$})\text{:=}2 t^3-3 t^2+1;\right) \left(h_2(\text{t$\\_$})\text{:=}3 t^2-2 t^3;\right) \left(h_3(\text{t$\\_$})\text{:=}t^3-2 t^2+t;\right) \left(h_4(\text{t$\\_$})\text{:=}t^3-t^2;\right) $$ And finaly the spline function: $$ P(\text{t$\\_$})\text{:=}\left\\{h_1 p_1.x+h_2 p_2.x+h_3 t_1.x+4 h_4 t_2.x;h_1 p_1.y+h_2 p_2.y+h_3 t_1.y+4 h_4 t_2.y\right\\}; $$ And trying to plot it: $$\text{ParametricPlot}[\\{P(t)\\},\\{t,0,1\\}]$$ But I'm getting some strange errors:               SetDelayed::write: Tag Plus in (1-3 t^2+2 t^3)[t_] is Protected. >>      What does it mean? Please provide me with a complete solution to plotting a hermite spline. Thanks! **UPDATE:** My mathematica code is:               Subscript[p, 0] = {143.0497; 195.6056};     Subscript[p, 1] = {159.3734; 181.8449};     Subscript[p, 2] = {160.8912; 179.2604};     Subscript[p, 3] = {176.423; 171.8781};          Subscript[t, 1] = (Subscript[p, 2] - Subscript[p, 0]) * 0.5;     Subscript[t, 2] = (Subscript[p, 3] - Subscript[p, 1]) * 0.5;          Subscript[h, 1][t_] := 2 * t^3 - 3*t^2 + 1;     Subscript[h, 2][t_] := -2*t^3 + 3*t^2;     Subscript[h, 3][t_] := t^3 - 2*t^2 + t;     Subscript[h, 4][t_] := t^3 - t^2;          P[t_] := {Subscript[h, 1][t] * Subscript[p, 1].x +          Subscript[h, 2][t] * Subscript[p, 2].x +          Subscript[h, 3][t] * Subscript[t, 1].x +          Subscript[h, 4] [t]* Subscript[t, 2].x;         Subscript[h, 1][t] * Subscript[p, 1].y +          Subscript[h, 2][t] * Subscript[p, 2].y +          Subscript[h, 3][t] * Subscript[t, 1].y +          Subscript[h, 4][t] * Subscript[t, 2].y};          ParametricPlot[{P[t]}, {t, 0, 1}]


title: How to express an expression with only ArcTan and ArcTanh?

text: I have an expression which is simply               (j/k) x^(j/k) LerchPhi[x,1,j/k)]       where `0 < j < k`. Manually I have been able (tedious work) to obtain only `ArcTan`'s and `ArcTanh`'s. How could I ask _Mathematica_ to do this automatically in the most compact form ?


title: Defining and evaluating new continuous distributions based on Mathematica's distributions

text: A number of functions operate with Mathematica's set of distributions e.g.               Through[{Mean, StandardDeviation, Variance}[UniformDistribution[{umin,umax}]]]      which returns the correct result                   {(umax + umin)/2, (umax - umin)/(2 Sqrt[3]), 1/12 (umax - umin)^2}      However, when I add two continuous distributions, thus creating a kind of annular distribution, it will not be evaluated.               Through[{Mean, StandardDeviation, Variance}[       UniformDistribution[{umin, umax}] -         UniformDistribution[{\[Epsilon] umin, \[Epsilon] umax}]]]      I think I had a very good look at the Mathematica documentation but I simply could not find any indication of how to define "my own distribution". Please note, that I know how to handle such functions, i.e., how to calculate the `Mean` ... The issue of this question is, are there ways in Mathematica, I have missed, that allow me to use Mathematica's functionality on distributions with my distributions I have defined using Mathematica's distributions.


title: Using a predicate to control the application of a rule in a replacement

text: How can I make a replacement rule that will only be applied to those parts of an expression that will not break a condition placed on the whole expression? For example, suppose I have $\frac{xy}{(x-y)(x+y)}$ and I want to change $y$ to $x$ getting the result $\frac{x^2}{(x-y)2x}$; i.e., without turning original expression into infinity. Further, I don't want the replacement to raise warnings or errors. At the least, I need a way to do a replacement just once, which seems like it should be very easy, but I have no idea how to do it. I came up with               ReplaceOnce[expr_, rule_Rule] := MapAt[# /. rule &, expr, Position[expr, rule[[1]]][[1]]]      How can I get to my goal from here?