title: Custom formatting for DateListPlot

text: I'd like to control more aspects of a `DateListPlot`, for example: shading for weekend days, and/or indicators for daytime/nighttime areas. By way of illustration, here's a simple example of a set of time data points (recent questions on mathematica.stackexchange):               questions =        First[Rest[         Import["http://api.stackexchange.com/2.0/questions?page=1&\     pagesize=100&order=desc&sort=creation&site=mathematica", "JSON"]]];          questionTimes =        Cases[questions, HoldPattern["creation_date" -> value_] :>  value,         Infinity];          hours = (Mod[#, 86400] /86400) 24  & /@ questionTimes ;          DateListPlot[Transpose[{questionTimes, hours}], Filling -> Bottom,       GridLines -> False, Frame -> {True, True, False, False}]      and the plot is like this: ![stack exchange](http://i.stack.imgur.com/xwQlR.png) I can't see how to show the necessary information along the x-axis, nor how to shade different areas of the graph to show day/night. _Edit_ : I now realise that the date/times returned by the SE API are in Unix Epoch (1970), and I hadn't noticed because I wasn't able to see the day numbers or years on my first attempt at a plot...


title: How does the speed of WriteString compare to other languages?

text: I'm doing something where the most sensible approach seems to be to open a file and write to it as I process data from some large tables. Old school. While testing the program, I remembered the low level write statements I use to use to do this kind of thing in FORTRAN and how fast they seemed to be, and I wondered how `WriteString` compared, so I did a little test. It took 10.875 seconds to write 1,000,000 records with record number, "Hello World." That seems pretty fast, but how does that compare to other languages? I use only Mathematica these days so I can't easily do any comparisons.               funit = OpenWrite[NotebookDirectory[] <> "writetest.csv"];     Timing[Do[WriteString[funit, i, ",Hello World\n"];, {i, 1, 1000000}]]          (*  {10.875000, Null}  *)


title: Problems using hand-made transformation functions

text: I experience difficulties when using custom transformation rules. Below I give two concrete problems. I am interested in fixing each of them and in understanding what is wrong with my attitude to the issue in general. 1) Suppose I try to make symbolic computations involving gamma function $\Gamma(z)$ more effective. Particularly, I want $\it{Mathematica}$ to make use of the property $\frac{\Gamma(z+1)}{\Gamma(z)}=z$ in simplification, which it doesn't do on its own. Evaluation of               Simplify[Gamma[z + 1]/Gamma[z]]      returns               Gamma[1 + z]/Gamma[z]      In order to handle this task let me define the first transformation function               tf1[e_] := e /. Gamma[x_]/Gamma[y_] /; Simplify[x - y] == 1 :> y;      Now evaluation of                Simplify[Gamma[z + 1]/Gamma[z],       TransformationFunctions -> {Automatic, tf1}]      indeed gives               z      However, execution of               Simplify[Gamma[z + 1]^2/Gamma[z]^2,      TransformationFunctions -> {Automatic, tf1}]      results in no real simplification and gives               Gamma[1 + z]^2/Gamma[z]^2      Well, one unwitty way to deal with this problem is to define a separate transformation function for the ratio of squares. But the ratio of the third powers will then require a third transformation function and so on. I would like to handle all these possibilities universally. I've tried something like               tfUniversal[e_] :=      e /. (Gamma[x_]/Gamma[y_])^n_ /; Simplify[x - y] == 1 :> y^n;      but this didn't work. Is there a way to solve this problem? 2) The second chapter of my adventures includes teaching "Mathematica" identity $\Gamma(z)\Gamma(1-z)=\frac{\pi}{\sin{\pi z}}$. Here problem comes right at the start. I define transformation function               tf2[e_] :=      e /. Gamma[x_] Gamma[y_] /; Simplify[x + y] == 1 :> Pi/Sin[Pi x];      and then evaluate               Simplify[{Gamma[-z] Gamma[1 + z], Gamma[z] Gamma[1 - z]},      TransformationFunctions -> {Automatic, tf2}]      which gives               {-\[Pi] Csc[\[Pi] z], Gamma[1 - z] Gamma[z]}      so that simplification occurs in the first case but not in the second. I wonder what is the problem here. Overall, I am very often troubled with such issues. Is it possible to find some references in the literature on using custom transformation rules? Any help is highly appreciated, thanks :)


title: Plotting TransformedDistributions of a Poly-Aeppli distribution

text: I do not know how to plot a `TransformedDistribution` of the named Polya- Aeppli distribution and do not know if the transformation can be handled after reading: Just what kind of transformations can TransformedDistribution handle?. I'm having trouble understanding the format of the plotting transformed distributions but hope that I can keep `TransformedDistribution` as part of the command.


title: How to use Interpolation to fill in missing data

text: I have a bit of data that looks like this:               data = {79., 79.1, 78.3, 78.6, 79., 77.2, 79.3, 79.3, "", "", 79.3, 76.6, \     79.3, 79.8, 78.4, 79., 78.4, 78.4, 79.4, 79.5, 79.4, 79.3, 79.2, \     79.2, 79.3, 79.6, "", 79.2, 79.3, 79., 79.1, 79.2}      I'm looking to do Linear Interpolation to fill in the missing data. The `Interpolation` functions with `InterpolationOrder->1` seems to be just what I am looking for. But I can't figure out how to feed the data to the `Interpolation` function in order to get a a sensible result.


title: How to define the Poisson bracket in Mathematica

text: If $f,g$ are functions of the independent variables $\\{q_1, q_2, ..., q_N, p_1, p_2, ..., p_N\\}$, then the Poisson bracket is defined as: \begin{equation*} [f, g] = \sum_{i=1}^N \left(\frac{\partial f}{\partial q_i}\frac{\partial g}{\partial p_i} - \frac{\partial f}{\partial p_i}\frac{\partial g}{\partial q_i}\right) \end{equation*} I'd like to define Poisson bracket in _Mathematica_. Just so to prove I'm not lazy, I wrote the following snippet (I doubt it is correct, let alone slick).               ClearAll["Global`*"];     pbra[f_, g_, q_, p_] :=         With[{c = Join[q, p]},          Sum[           D[f@@c, q[[k]]] D[g@@c, p[[k]]] - D[f@@c, p[[k]]] D[g@@c, q[[k]]],          {k, 1, Length@q}]]      And call it like (dummy example):               f[x_, y_, px_, py_] := x+y+px*py;     h[x_, y_, px_, py_] := x-y+px/py;     p[f, g, {x, y}, {px, py}]      How would _you_ define Poisson bracket ? For example, I am wondering whether it would be more sound to define it in a way that the user calls it like:               f[x_, y_, px_, py_] := x+y+px*py;     h[x_, y_, px_, py_] := x-y+px/py;     p[f[x, y, px, py], g[x, y, px, py]]


title: Solving a constrained vector system with NDSolve

text: I'm trying to replicate the pendulum DAE example in the Documentation Center (example/ModelConstrainedSystemsAsDAEs) using Mathematica 9's vector support with NDSolve. So far I have:               deqns = {m X''[t] == \[Lambda][t]-{0, m 9.81}};     constraints = {Norm[X[t]] == 1};     initialConditions = {X[0] == {1, 0}, X'[0] == {0, 1}};     params = {g -> 9.81, m -> 1, l -> 1};          pendulumSol =       First[NDSolve[{deqns, constraints, initialConditions} /.         params, {X, \[Lambda]}, {t, 0, 15},           Method -> {"IndexReduction" -> Automatic}]];      but am getting:               NDSolve::ndincd: "Initial conditions for derivatives of the                        function X[t] do not have consistent dimensions."      Is this an issue with NDSolve properly identifying the vector nature of the equations, or have I overlooked something? Any help would be appreciated.


title: Need help in building a pattern

text: Suppose I have the following list:               list = {w, r, {f, first, t, y}, {a, y}, last, {k, r, first, e}, m, e, last, r, t}      I am expecting this               Cases[list, {__, {__, first, x__}, y__, last, __} :> {{first, x}, y,         last}]      to result in this               (*{{{first, t, y}, {a, y}, last}, {{first, e}, m, e, last}}*)      and this               Cases[Flatten@list, {__, first, x__, last, __} :> {first, x, last}]      to result in this               (*{{first, t, y, a, y, last}, {first, e, m, e, last}}*)      but the results for both are `{}`. How can I build a suitable pattern to get the desired results?


title: Sporadic bugs in multidimensional SparseArray objects

text: **Udate** This `SparseArray[]` bug is (apparently) fixed on version 10.0               $System//Print;     $Version//Print;          Mac OS X x86 (64-bit)     10.0 for Mac OS X x86 (64-bit) (June 29, 2014)      Thanks, Wolfram! * * * **Resolution**   * WRI affirms that `SparseArray[__]` is broken (in 9.0 and 9.01)   * WRI's suggested workaround is `SparseArray//Normal` My experience with this particular `SparseArray` bug has been that reliable coding requires that the `Normal` workaround be applied to _every_ instance of `SparseArray`, because the failures of `SparseArray` are sporadic, occult, and associated to unknown preconditions. Ouch! My thanks are extended to all who commented. * * *               Wolfram Technical Support      email: support@wolfram.com     Re: [CASE:646986]          Sporadic bugs (?) in multidimensional `SparseArray` objects          Hi John,          Thank you for getting back to us.          I understand that for now your temporary workaround      of `aTest = aTest // Normal` is holding the fort.          Our developers are working on this issue and      we anticipate that the resolution will appear      in the future release of Mathematica.          Sincerely,          <name redacted>          Wolfram Technology Group     Wolfram Research, Inc.       * * * What appear to be sporadic bugs (?) in multidimensional `SparseArray` objects (?) are plaguing our code development efforts.               $Version//Print;     9.0 for Mac OS X x86 (64-bit) (January 24, 2013)       The following code elicits wrong answers and/or kernel crashes:               Do[         (* aTest is a SparseArray[...] that crashes the kernel and/or            yields incorrect results unless wrapped in Normal[...] *)              aTest = RandomReal[NormalDistribution[0,1], {2,2,2,2}] +             I*RandomReal[NormalDistribution[0,1], {2,2,2,2}]//           SparseArray[               {                   {i_,j_,k_,1 ,i_} :> #[[k,1,i,j]]*#[[k,2,1,j]],                   {i_,j_,k_,2 ,j_} :> #[[k,1,i,j]]*#[[k,2,2,i]],                   {i_,j_,k_,n_,j_} :> #[[k,1,i,j]]/;Equal[n-2,i]               },               {2,2,2,4,2},               Complex[0.0,0.0]           ]&//Flatten[#,1]&//             Transpose[#,{4,1,2,3}]&//               Flatten[#,2]&//                 Transpose;              (* enable the following line to "fix" aTest *)         (* aTest = aTest//Normal; *)              (* the following all should be 0.0 (to machine precision) *)              (aTest//Normal)-aTest//           Flatten//Abs//Max//Print["      (aTest//Normal)-aTest =?= 0 = ",#]&;              id = SparseArray[{i_,i_}\[Rule]1,{16,16},0];         (aTest//Normal).id-aTest.id//           Flatten//Abs//Max//Print["(aTest//Normal).id-aTest.id =?= 0 =",#]&;     ,{10}];      **Results** (when the kernel doesn't crash immediately)                     (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =7.43925           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =3.14861           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =1.27373           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =1.8666           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =2.03269           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =3.26388           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =1.30013           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =3.14138           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =0.953447           (aTest//Normal)-aTest =?= 0 = 0.     (aTest//Normal).id-aTest.id =?= 0 =0.841252      **Remarks** Variants of the above code crash/fail for Complex-valued and Integer-valued SparseArrays; apparently the multidimensionality is the essential trait. Help, advice, and workarounds all are appreciated!


title: String pattern matching in function definition

text: I have a function f and its downvalues. I want to code an extra downvalue for arguments of the form "part1:part2" I don't understand why               f[a__~~":"~~b__] := {a,b}      does not work. I can only think of               f[string_String] /; StringMatchQ[string,__~~":"~~__] := StringCases[string, a__~~":"~~b__ :> {a,b}][[1]]      but this doesn't seem very smart as the string structure is examined two times.


title: Most influential and/or disrupting features introduced with new versions of Mathematica

text: This is something that has been bugging me for ages. With every upgrade, _Mathematica_ has added new functionalities and amendend some of its bugs. In the process, though, some older functionality has inevitably been disrupted or made obsolete. This is to be expected of course, and more often than not these modifications and additions are only limited to a change in names, or very limited changes in behaviour (with some completely transparent to users). Some changes, though, have a deeper influence on the way we code and use the front-end with _Mathematica_. Two major game changers (IMBO) have been the introduction of the two-dimensional notation in all versions of Mathematica 3 (Traditional Form notation), and the graphics revolution in Mathematica 6 (graphics are no longer side effects, but real objects). Other important changes that come to mind are the introduction of the packed-array facility in Mathematica 4, the continuing evolution of numerical evaluation through all versions, and the availability of curated data from the Web. The incorporation of most external packages in the core is another game changer, since it contributed to hiding the code and made it harder for users to understand what is going on under the hood. Also, the introduction of Manipulate and of the CDF format have changed the way people use _Mathematica_ to convey ideas. I am well aware that there are sections of the documentation dedicated to _all_ these changes, namely the " **What's new in Mathematica X.x** " parts, that can be found here http://www.wolfram.com/mathematica/quick-revision- history.html and more specifically (for integer values of version numbers) here: > http://www.wolfram.com/mathematica/new-in-9/ > > http://www.wolfram.com/mathematica/new-in-8/ > > https://www.wolfram.com/products/mathematica/newin7/ > > http://www.wolfram.com/products/mathematica/newin5/summary.html > > _?? version 4 what's new ??_ > > http://reference.wolfram.com/legacy/v3/index105.html > > _older versions seems to lack a dedicated part of the documentation_ The purpose of this 'question' is to create a **_selected_** listing of the most disrupting and/or influential changes introduced with every new version. The utility of such a listing is multifold:   1. Understanding older code (need to adapt some old program? reading an old book on Mathematica?)    2. Getting rid of older coding habits (What, you're still using $DisplayFunction?)    3. Getting to know new, more efficient ways, to do stuff (Maybe Total is better than Plus@@ with Packed Arrays?)    4. Having a list _ordered_ by perceived importance of the most relevant changes (Yeah, let's see the really important changes of paradigm that I have overlooked...) The template for the answers, in my opinion, should carry a title like this > **From X.x to Y.y : name of the feature of behaviour** for example (these are just examples of the titles of three separate answers)                   From 2.2 to 3.0 : Frontend with two dimensional notation         From 3 to 4 : Packed Arrays         From 5.1 to 5.2 : 64-bit-enhanced arbitrary-precision numerics      And then elaborate what is that has changed, why is that important and how it has changed the way you coded or used Mathematica. **ONE** answer, **ONE** feature. My idea is to have answers to this post pop up naturally by voting, and then from time to time update this very post with a list by importance for each version change. As many other threads of this kind, this is a Collaborative Wiki post.


title: Saving a graph of a great number of nodes

text: I'm developing a code where I get graphs with over 10K nodes. These graphs have a particular relationship between clustering and degree that I need to save for future analysis. The problem is when I use `EdgeList` and then save the list in a _Mathematica_ file (.nb), when I recall it, the size of some communities are different to the ones showed by the graph. The following sizes are obtained after change the graph by the `EdgeList`               Community = FindGraphCommunities[g];     aux = Length[#] & /@ FindGraphCommunities[g]      > >     {1845, 1656, 1536, 1419, 1206, 441, 349, 342, 210, 205, 168, 126, >               113, 103, 73, 58, 57, 45, 40, 8} >   Then I use `EdgeList[g]`, click on Show full Output and copy the entire cell and save it into the file.nb. When I call again the graph `g3` from the file.nb it gives me the following sizes               g = g3;     Community = FindGraphCommunities[g];     aux = Length[#] & /@ FindGraphCommunities[g]      > >      {1884, 1766, 1720, 1364, 643, 420, 420, 394, 266, 218, 211, 159, 117, >             101, 93, 76, 29, 28, 24, 21, 20, 16, 10} >   I don't know what happen there, could be the command `EdgeList` or something behind...Then I try to save the graph by Exporting it as a `gml` file and after 10 minutes it was still running... So, do you have any suggestion for this, how can I export efficiently this information or do I have to wait for the exportation in `gml` format?


title: Cropping a Plot without rasterization

text: Suppose that I make a plot of $x^2$ vs. $x$:               myPlot = Plot[x^2, {x, -1, 1}, Frame -> True]      ![myPlot1](http://i.stack.imgur.com/PkklF.png) If I click on the plot, I see the orange handles with which I can resize the plot: ![myPlot2](http://i.stack.imgur.com/xfyTX.png) Now, suppose that I want to actually crop the plot -- I would like to effectively shrink the orange box so that some of the plot is actually cropped. ( _Why in the world would I want to do this?_ It's a long story...) I have found, from answers to a question that I previously asked, that this is possible using `ImagePad` with negative arguments. For example, suppose that I want to crop 40 pixels off the right side of the plot. I can do the following:               ImagePad[myPlot, {{0, -45}, {0, 0}}]      ![myPlot3](http://i.stack.imgur.com/Fh4Vy.png) It looks good so far. **_However,_** when I resize the resulting cropped plot by dragging the orange handles, it looks like the image has been rasterized: ![myPlot4](http://i.stack.imgur.com/BhwSu.png) **But, in contrast, I need the image to remain a vector image (non-rasterized) after cropping.** Is that possible? I am running both version 8 and 9.


title: Manipulate strings within a list

text: Let's say I have               Solve[x^4 + 3 == 0, x]      with output               {{x -> -(-3)^(1/4)}, {x -> -i (-3)^(1/4)}, {x ->i (-3)^(1/4)}, {x -> (-3)^(1/4)}}      how do i lose the "x->" part of each string within this list and get something like               {-(-3)^(1/4), -i (-3)^(1/4), i (-3)^(1/4), (-3)^(1/4)}      Thanks for the answers/links. Judging by the answers I realise now that I have formulated my question poorly. I wanted to know how I lose symbols that are in front of a specific symbol. Let's have a look at two more examples:               list={a: horse, b: chicken, c: fish}      how do I lose "a: ","b: ","c: " or               list2={section 1, section 2, section 3}      how do I lose "section"


title: SSH known hosts in Windows (Error code = 124)

text: I had to re-install Raspbian OS and now it has new SSH key. SSH clients usually let ignore this fact and save the new key, but not now. Mathematica gives no way to go if the Host Key was changed. I get an error               SSH received a different Host Key while launching kernel 'Raspberry Pi'. Error code = 124.      when trying to start remote kernel and I have found no any information about where is the Mathematica's SSH client on Windows... I just need to clear known hosts... The easy way is to change Raspberry's IP, but it's far less interesting, moreover not always possible


title: VertexLabelStyle not found in Mathematica 8.0.0.0

text: I tried to use `VertexLabelStyle`, but seemed not recognized in my system with version number 8.0.0.0. Is this expected or I have missed something? From the online documentation page, I don't see this option was added after version 8.0.0.0; it just says "New in version 8".


title: Manipulate[] inside DialogInput[]

text: I would like to get input from an interactive dialog that uses `Manipulate[]`. But I find that my `Manipulate[]` breaks when I put it inside a `DialogInput[]`. So, when I try `DialogInput[Manipulate[Plot[Sin[e x], {x, 0, 2 \[Pi]}], {e, 1, 10}]]` the plot inside the dialog is not connected to the `Manipulate[]` slider. What am I missing here?


title: Ragged Diamonds

text: diamond[n_Integer] :=       Table[1, {#}] & /@ (Range[n]~ Join~Range[n - 1, 1, -1])          diamond[5] // MatrixForm      ![enter image description here](http://i.stack.imgur.com/e0sFG.jpg) How could a more functional, "tablefree", solution look like? **EDIT** Thanks for the nice comments and answers so far :) But I'm looking for a `1-2-3-4-3-2-1` or `1-2-3-4-5-4-3-2-1` sequence which `DiamondMatrix` doesn't seem to offer.


title: Why round to even integers?

text: According to the Mathematica help: > Round rounds numbers of the form x.5 toward the nearest even integer. For example:               Round[{0.5, 1.5, 2.5, 3.5, 4.5}]      gives               {0, 2, 2, 4, 4}      What's the rationale behind this? Why not the usual x.5 always rounds up?


title: How can I reduce computation time while still obtaining a good approximation for my function?

text: I am new to any CAS (and Mathematica, for that matter) and new to StackExchange too, so forgive me and correct me on any mistakes. I have this function: $J_p=\sum_{m,n=1}^{\infty} \epsilon_{mn}f_{mn}\sum_{k=-\infty}^{\infty}\frac{J_k^2(\beta)(m\Omega+k\omega)}{1+(m\Omega+k\omega)^2}$ where $\epsilon_{mn}=-\frac{m n}{4\pi^2}\int_0^{2\pi}\epsilon(p_x,p_y)\exp(-i(m p_x+n p_y))\,dp_x dp_y$ and $f_{mn}=-\frac{m n}{4\pi^2}\int_0^{2\pi}\frac{\exp(-i(m p_x+n p_y))}{1+\exp(-\epsilon(p_x,p_y))}\,dp_x dp_y$ where again $\epsilon(p_x,p_y)=\sqrt{1+4\cos\left(\frac{p_y}{2}\right)\cos\left(\frac{p_x\sqrt{3}}{2}\right)+4\cos^2\left(\frac{p_y}{2}\right)}$. Here is my Mathematica code to evaluate this:               Off[NIntegrate::ncvi];     epsilonCoeffsMMA[cl_] := Module[{reComp, imComp},        reComp[m_, n_] := (-m n)/(4 \[Pi]^2)           NIntegrate[           Re[(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/2)              Exp[-I (m px + n py)]], {px, 0, 2 \[Pi]}, {py, 0, 2 \[Pi]},            Method -> "Trapezoidal", MaxRecursion -> 100];        imComp[m_, n_] := (-m n)/(4 \[Pi]^2)           NIntegrate[           Im[(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/2)              Exp[-I (m px + n py)]], {px, 0, 2 \[Pi]}, {py, 0, 2 \[Pi]},            Method -> "Trapezoidal", MaxRecursion -> 100];        emnMatrix = Table[0, {m, 1, cl}, {n, 1, cl}];        Do[emnMatrix[[m, n]] = reComp[m, n] + I imComp[m, n], {m, 1,           cl}, {n, 1, cl}];        ];     boltzECoeffsMMA[cl_] := Module[{reComp, imComp},     reComp[m_, n_] := (-m n)/(4 \[Pi]^2)           NIntegrate[           Re[Exp[-I (m px + n py)]/(1 +                Exp[-(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] +                    4 Cos[py/2]^2)^(1/2)])], {px, 0, 2 \[Pi]}, {py, 0,             2 \[Pi]}, Method -> "Trapezoidal", MaxRecursion -> 100];        imComp[m_, n_] := (-m n)/(4 \[Pi]^2)           NIntegrate[           Im[(1 + 4 Cos[py /2] Cos[(px Sqrt[3])/2] + 4 Cos[py/2]^2)^(1/2)              Exp[-I (m px + n py)]], {px, 0, 2 \[Pi]}, {py, 0, 2 \[Pi]},            Method -> "Trapezoidal", MaxRecursion -> 100];        fmnMatrix = Table[0, {m, 1, cl}, {n, 1, cl}];        Do[fmnMatrix[[m, n]] = reComp[m, n] + I imComp[m, n], {m, 1,           cl}, {n, 1, cl}];        ];     jPMMA[coeffLim_, kernLim_] := Module[       {cl = coeffLim, kl = kernLim, px, py},       epsilonCoeffsMMA[cl]; boltzECoeffsMMA[cl];       coeffMatrix = emnMatrix fmnMatrix;       sumMatrix =         Table[Sum[(          BesselJ[k, \[Beta]]^2 (m \[CapitalOmega] + k \[Omega]))/(          1 + (m \[CapitalOmega] + k \[Omega])^2), {k, -kl, kl}], {m, 1,           cl}, {n, 1, cl}];       jParaMMA = Total[coeffMatrix sumMatrix, 2];       ];      This generates a function `jParaMMA` which I can `Plot` after I have made the call `jPMMA[a,b]` for some integers; `a` and `b`. For example               jPMMA[10, 10];     Plot[Evaluate@     Re[jParaMMA /. {\[Beta] ->2, \[Omega] -> {0, 2, 4, 6, 8}}], {\[CapitalOmega], 0, 20}, PlotRange -> Full]      ![Plot of the function jParaMMA for \[Beta\]->2 and \[Omega\] -> {0, 2, 4, 6, 8}](http://i.stack.imgur.com/ZOwiT.jpg) for which               First[Timing[jPMMA[10, 10]]]      gives               115.437500      My question is: How can I obtain similar results, possibly with more terms (i.e. from running `jPMMA[50, 60]`, say.) in a shorter time,? Thank you. PS: I used the `Off[NIntegrate::ncvi]` because I do not know how to eliminate it from my numerical integration and I'd be glad to obtain some help for that too. Also, I used the `Trapezoidal` method because I noticed it gave a faster approximation even when coupled with `MaxRecursion -> 100`. I have tried with the Cuba library implementation in both Mathematica and Maple, which I was led to by this post, and the approximations are appreciably close.


title: Publishing a Mathematica notebook

text: I am considering publishing a Mathematica notebook, used for analysis of some data, as a “supporting information” file alongside a research paper. What are the common practices for doing so, and what should I consider in making this choice? In particular, I am wondering:   * Is the `.nb` format suitable, or should I transform it into a CDF?   * Or is it better to include a PDF of the Mathematica notebook (thus readable by anyone), attaching the notebook itself inside the PDF (thus allowing Mathematica owners/users to use it directly)?   * Do I need to audit/edit the content of the file to remove any metadata that could leak?


title: Joining and interpolating data points

text: I wonder what is the best practice for interpolating curves? Usually I'm using `BSplineCurve` and adjusting `SplineWeights` so it would fit better (assigning more weight around sharp edges to drag the curve closer to it). Or if I can guess what formula describing the points I use `FindFit`. But often I can't guess the formula and adjusting weights is very tedious so it's easier just manually draw the curve along the points. So what is the best way to join points in _Mathematica_? Consider these five sets of points describing five resonance curves for example:                data = {{{100.434, 0.}, {102.461, 0.0909091}, {104.392, 0.318182}, {105.321, 0.545455},               {106.226, 1.}, {107.108, 0.545455}, {107.965, 0.318182}, {109.608, 0.136364},               {111.154, 0.0909091}},              {{100.434, 0.}, {102.461, 0.06}, {104.392, 0.22}, {105.321, 0.46}, {106.226, 1.},               {107.108, 0.4}, {107.965, 0.12}, {111.154, 0.02}, {113.958, 0.}},              {{100.434, 0.030303}, {102.461, 0.0505051}, {104.392, 0.0909091},               {105.321, 0.272727}, {105.867, 0.494949}, {106.226, 1.}, {106.582, 0.636364},               {107.108, 0.212121}, {107.965, 0.0505051}, {111.154, 0.}},              {{100.434, 0.}, {102.461, 0.0555556}, {104.392, 0.333333}, {105.321, 0.611111},               {105.867, 1.}, {106.226, 0.944444}, {106.405, 0.583333}, {106.582, 0.777778},               {106.933, 1.}, {107.965, 0.444444}, {109.608, 0.166667}, {111.154, 0.0555556}},              {{100.434, 0.0188679}, {102.461, 0.0566038}, {104.392, 0.}, {105.321, 0.54717},               {105.867, 0.849057}, {106.226, 1.}, {106.405, 0.433962}, {106.582, 0.886792},               {106.933, 0.924528}, {107.281, 0.660377}, {107.965, 0.320755},               {111.154, 0.0566038}}}


title: Nonlinear FEM and FindRoot

text: I'm trying to develop a kind of nonlinear FEM application using mathematica to solve a bvp like the following: $$ \gamma(u') ~u^{iv} + 2 \gamma'(u') u''' u''+ u''^3 = f(x) $$ where $u = \tilde{u}(x)$, $\gamma = \tilde{\gamma}(u')$ and $'$ denotes derivative with respect to the function variable, thus: $$ u' = \frac{\text{d} u(x)}{\text{d}x} \quad \quad \gamma' = \frac{\text{d} \gamma(u')}{\text{d} u'} $$ The next step is getting the weak form: $$ \int_{\Omega} \left( \gamma ~u^{iv} + 2 \gamma' u''' u''+ u''^3 \right) \psi = \int_{\Omega} f(x)\psi $$ where $\Omega$ is the domain and $\psi = \tilde{\psi}(x)$ is the test function. In my code, I'm using a _Cubic Hermite Shape Functions_ , then the approximate function reads as follows: $$ \hat{u} = \sum_i \alpha_i \psi_i^{(a)} + \beta_i \psi_i^{(b)} $$ where $\psi_i^{(a)}$ and $\psi_i^{(b)}$ are piecewise functions, and $\alpha_i$ and $\beta_i$ are coefficients related to the values of function and its first derivative, respectively. Let bear in mind the chain rule: $$ \frac{\text{d} \gamma(u')}{\text{d} x} = \frac{\text{d} \gamma(u')}{\text{d} u'} \frac{\text{d} u'}{\text{d} x} = \gamma' u'' $$ After some calculation, the weak form assumes the following form: $$ \int_{\Omega}\gamma~u^{iv} = \left[\left(\gamma u''' + \frac{1}{2}\gamma' u''^2 \right)\psi-\left(\gamma \psi \right)' u'' \right]_{\partial \Omega} + \int_{\Omega} \frac{1}{2} \gamma'' u''^3 \psi + \frac{3}{2} \gamma' u''^2 \psi' + \gamma u'' \psi '' $$ $$ \int_{\Omega} 2 \gamma' u''' u'' \psi = \left[\gamma' u''^2 \psi \right]_{\partial \Omega} - \int_{\Omega} \gamma'' \psi u''^3 \psi + \gamma' u''^2 \psi' $$ As starting point, however, I'm using a linear equation, setting $\gamma(u') = 1 $ (but at the end, the procedure is the same). Then, the simple equation is: $$ u^{iv} = f(x) \Longrightarrow \int_{\Omega} u^{iv} \psi = \int_{\Omega} f(x) \psi $$ then $$ \left[u''' \psi - u'' \psi' \right]_{\partial \Omega} + \int_{\Omega} u'' \psi'' - \int_{\Omega} f(x) \psi = 0 $$ If the domain is discretised in _n_ points ( _n-1_ elements), this equation represents a system of _2n_ (non)linear equations. Let suppose to apply BCs on the function and its first derivative at both edges. The system reads as follows: $$ R(1) = \alpha_1 - u(x_0) \\\\\\\ R(2) = \beta_1 - u'(x_0) \\\\\\\ ... \\\\\\\ R(i) = \int_{\Omega_i} u'' \psi'' - \int_{\Omega_i} f(x) \psi \\\\\\\ ... \\\\\\\ R(2n-1) = \alpha_n - u(x_1) \\\\\\\ R(2n) = \beta_n - u'(x_1) \\\\\\\ $$ where the domain is $\Omega = [x_0, x_1]$. Now the goal is solving this system (in general nonlinear). I'm trying to do that using FindRoot. Thus I defined the following functions:               (* integration of forcing term *)     funFF[X_, f_] := Module[       {F, n},       n = Length[X];       F = ConstantArray[0, 2 n];       Do[         jr = 2 j - 1;         F[[jr ;; jr + 3]] =          F[[jr ;; jr + 3]] +          Flatten[NIntegrate[             f[x] NN0[x, X[[j]], X[[j + 1]]], {x, X[[j]], X[[j + 1]]}]],         {j, 1, n - 1}         ];       F       ];            (* assembly of system *)       SYS1[X_, U_?VectorQ, RHS_, fun_] := Module[         {EQ, x1, x2, f1, p1, f2, p2, F, P, R, n},         n = Length[X];          (* value of function *)         F = U[[1 ;; Length[U] ;; 2]];         (* value of first derivative *)         P = U[[2 ;; Length[U] ;; 2]];         (* initialitazion *)         EQ = ConstantArray[0, 2 n];         (* starting loop *)         Do[           jr = 2 j - 1;           (* values for the current element *)           x1 = X[[j]];           x2 = X[[j + 1]];           If[Length[F] == 0,             f1 = F;             f2 = F,             (* if non-constant *)             f1 = F[[j]];             f2 = F[[j + 1]]           ];           If[Length[P] == 0,             p1 = P;             p2 = P,             (* if non-constant *)             p1 = P[[j]];             p2 = P[[j + 1]]           ];           (* j-th equation *)           EQ[[jr ;; jr + 3]] = EQ[[jr ;; jr + 3]] +                   Flatten[NIntegrate[fun[x, x1, x2, f1, p1, f2, p2], {x, x1, x2}]],           {j, 1, n - 1}         ];         (* assembly of system *)         R = EQ - RHS;         (* boundary conditions *)         R[[1]] = U[[1]] - wf0;         R[[2]] = U[[2]] - wp0;         R[[2 n - 1]] = U[[2 n - 1]] - wf1;         R[[2 n]] = U[[2 n]] - wp1;         (*output*)         R     ];      where I assumed that $w(x)$ is the exact solution, thus `wf0 = w(x0)`, `wp0 = w'(x0)`, `wf1 = w(x1)$` and `wp1 = w'(x1)`. Then:               (* forcing term *)     f1[x_] = w''''[x];          (* integration of forcing term *)     F1 = funFF[XX, f1];          (* function to be integrated*)     fun1[x_, x1_, x2_, f1_, p1_, f2_, p2_] := Ne2[x, x1, x2, f1, p1, f2, p2] NN2[x, x1, x2];          (* solution *)     U1 = UU /. Flatten[FindRoot[SYS1[XX, UU, F1, fun1] == 0, {UU, U0}]      where:   * `XX` is the grid   * `U0` = {\alpha_1, \beta_1, \alpha_2, \beta_2, ... } is the guess   * `Ne2` contains the interpolation inside an element   * `NN2` contains the four shape functions   * `fun` is the integrand for each element, i.e $u'' \psi ''$ This simple example works, but it is really slow even if it's a linear case. How can I improve/optimize the code? I'm new in this kind of programming, so I'm sure I'm missing some speeding-up tricks. Any suggestion is very appreciated. Best, Petrus


title: Speeding Up NIntegrate and ?NumericQ

text: I have attached below a piece of code that has been running successfully but takes extremely long for mathematica to compute. I believe it is the last "rms" calculation and the fact that I am using very small numbers that is really pushing mathematica. Any suggestions or help would be greatly appreciated.               c = 3*10^8;     e = 1.6*10^-19;     m = 9.1*10^-31;     rc = 2.8*10^-15;     gamma = 1600/0.511;     n = (10^-10)/e;     sbeg = 9.6*10^-5;     send = 1.1*10^-5;     lm = 0.5;     theta = 0.04;     r = lm/theta;               sL[phi_?NumericQ] := (r*phi^3)/24          lambda[z_, phi_] := n/(((send - sbeg)/theta phi + sbeg)*Sqrt[2 Pi])        Exp[-(z^2/(2*((send - sbeg)/theta phi + sbeg)^2))]               de[s_, phi_?NumericQ] :=  With[{sL = sL[phi]},      (-((2*n*rc*m*c^2)/(3^(1/3) r^(2/3))))      *((lambda[s - sL, phi] - lambda[s - 4 sL, phi])/sL^(1/3)       + NIntegrate[(1/(s - sprime)^(1/3)*D[lambda[sprime, phi], sprime]),            {sprime, s - sL,s}])]          total[s_?NumericQ] := NIntegrate[de[s, phi], {phi, 0, theta}]          rms = Sqrt[       NIntegrate[         1/(sbeg*Sqrt[2 Pi])           Exp[-(s^2/(2*(sbeg)^2))]*(total[s])^2, {s, -10 sbeg,           10 sbeg}] - (NIntegrate[          1/(sbeg*Sqrt[2 Pi])            Exp[-(s^2/(2*(sbeg)^2))]*(total[s]), {s, -10 sbeg,            10 sbeg}])^2]


title: How do I plot x[t] vs. x'[t] (where x[t] and x'[t] are solutions to NDSolve)?

text: I have a differential equation which I solved using NDSolve. I can easily plot x[t] vs. t, x'[t] vs. t, but.... how do I plot x[t] vs. x'[t]? I tried using the Evaluate function to simplify things, but I still have no luck. Here's what I mean:               x1 = Evaluate[x'[t] /. sol];     x2 = Evaluate[x[t] /. sol];          Plot[x2, {x1, 0, 50}, PlotRange -> Automatic, AxesLabel -> {x[t], x'[t]}]      How can I plot these? This approach also did not work:               Plot[x[t]/.sol, {x'[t]/.sol, 0, 50}, PlotRange -> Automatic,      AxesLabel -> {x[t], x'[t]}]      Help!


title: About generating power series

text: For an arbitrary function $f(x,y)$ I am defining functions `LogMT1` and `LogMT2` as follows,                Nn = 5;            LogMT1 = Sum [ f[x^n, y^n]/(n*(1 - x^(2*n))), {n, 1 , Nn}];             LogMT2 =Sum[Log[(1 + x^n)/((1 - x^(n/2)*y^(n/2))*(1 - x^(n/2)*y^(-n/2)))], {n, 1 , Nn}];      Now I want to know the power series expansion of both these functions LogMT1 and LogMT2 as a power series in x - and I expect them to come as $\sqrt{x}$, $x$, $x\sqrt{x}$, $x^2$ and so on and each of them should be multiplied by a function of $y$ as a coefficient.   * I want to know how this can be done?  (..of course my eventual goal is to be able to determine $f(x,y)$ such that `LogMT1 = LogMT2` for arbitrarily large values of $Nn$ and it would be great if someone can suggest a Mathematica way of being able to do that...) * * * Here is a function $f$ which seems to solve the above equation for arbitrarily large values of $Nn$ to arbitrarily large values of powers of $x$ as far as one can see this way, (..the point is that I don't know how this function $f$ can be derived..)                  $Assumptions = y > 0;             f[x_, y_] =  Sqrt[x] (Sqrt[y] + 1/Sqrt[y]) + x (1 + y + 1/y) + x^(3/2) (y^(3/2) + 1/y^(3/2)) + x^2 (y^2 + 1/y^2) + ((x y)^(5/2) (1 - 1/y^2))/(     1 - Sqrt[x y]) + (x/y)^(5/2)/(1 - Sqrt[x/y]) (1 - y^2) // Simplify;           Nn = 30;(*you can increase this but it takes longer time*)          LogMT1 = Sum[f[x^n, y^n]/(n (1 - x^(2 n))), {n, 1, Nn}];          LogMT2 = Sum[ Log[(1 + x^n)/((1 - x^(n/2) y^(n/2)) (1 - x^(n/2) y^(-(n/2))))], {n, 1, Nn}];          Series[LogMT1 - LogMT2, {x, 0, Nn/2}] // Simplify      O(x^{31}) is the output showing that the equation is satisfied to that order.


title: How to 'center' a Histogram

text: How can a generate a histogram so the peak of the distribution is at x=0? I want to center these histograms for comparison               (*k=2*)     Histogram@Table[Total@RandomInteger[1, {2}], {1000}]          (*k=4*)     Histogram@Table[Total@RandomInteger[1, {4}], {1000}]          (*k=10*)     Histogram@Table[Total@RandomInteger[1, {10}], {1000}]          (*k=50*)     Histogram@Table[Total@RandomInteger[1, {50}], {1000}]          (*k=100*)     Histogram@Table[Total@RandomInteger[1, {100}], {1000}]          (*k=10 000*)     Histogram@Table[Total@RandomInteger[1, {10000}], {1000}]          (*k=100 000*)     Histogram@Table[Total@RandomInteger[1, {100000}], {1000}]      edit: range was incorrect edit 2: Got it! Histogram[RandomVariate[NormalDistribution[0, 1],k]]


title: Listable compiled functions used in NMinimize

text: This is an example from MMA:               cListable = Compile[{x}, (x-1)^2, RuntimeAttributes -> {Listable}]      In case of single entry:               cListable[ 10.1]      Returns:                82.81      and if we give it a list:                cListable[ {-5, 0, 5}]      Returns:                {36., 1., 16.}      Is it possible to use cListable, or any other compiled function, in `NMinimize` so that the function I want to minimize take a list as its argument and as a result use all of my CPU cores to minimize the function? The point of doing this is using all the CPU cores available and reduce the computation time.


title: Export ContourPlot as DXF

text: I'm having some trouble exporting a ContourPlot as a .DXF (to be used in AutoCAD). My code is:               LG3 = ContourPlot[1/2 - Sum[Sinc[(n \[Pi])/2] Cos[n (1 x + 4 ArcTan[y, x])], {n, 1, 200}] == 0.5,       {x, -20, 20},       {y, -20, 20},       PlotPoints -> 100]     SetDirectory[NotebookDirectory[]];     Export["LG3.dxf", LG3];      The relevant output is: > Export::nodta: Graphics3D contains no data that can be exported to the DXF > format. >> Any help would be appreciated. Thanks.


title: Performance tuning in Mathematica?

text: What performance tuning tricks do you use to make a _Mathematica_ application faster? MATLAB has an amazing profiler, but from what I can tell, _Mathematica_ has no similar functionality.


title: Faster evaluation of an expression for brute force optimization

text: I'm trying to bruteforce an optimization problem and I want to generate roughly a million function evaluations. I was wondering what would be the best for me to do in terms of having the expressions evaluated faster. Curently the expressions are in plain symbolic form and I tried using the ReplaceAll[] or /. operator and it runs really slowly. My expressions are quite large, when I output them into a text file that amasses to a few thopusand lines of equations, so yeah, it's kinda slow because of that. I did optimize it a little in terms of that when in bruteforce, only a single variable of an expression evaluates for a For loop at certain depth and the entire expression doesn't really evaluate a gazillion times, but it's still kinda slow. I know there are a trillion ways to tell mathematica the same thing (turn the expression into a function, etc..) and I was wondering what is the type of definition that would make the evaluation run faster. tl;dr So basically, a simple, somewhat noobish solution to the problem of faster expression evaluation is what I'm looking for


title: Automatic Cuda/CPU code; Fast programme run

text: I have developed my thesis codes by means of Mathematica software in a .nb file format. Despite using the available Parallel code capabilities and doing all the numerical integrations by means of parallel programming; the duration of programme run is so so much high; something like 2 to 3 weeks; which makes the debugging and verification part of the thesis almost unavailable. My question goes to two branches:   1. How can I make my programme more efficient and run faster.   2. I have heard that by applying the Cude or GPU programming; the duration of run will significantly be reduced. How is it possible to automatically change and turn my present codes to a Cude GPU code? As far as I am facing lack of time for my thesis; answers to my questions and helping me decrease the run duration will be so appreciated.


title: How to avoid conflicts between local variable names and symbolic arguments in Block constructs?

text: I am using `Block` inside the definition of some recursive functions. If the Block definitions use short variable names, like `{s = First[param], r = Rest[param]}`, an infinite recursion occurs when `param` contains one of the symbols `s` or `r`. Using tortured local symbols names minimizes this risk, but   * it makes code unreadable,   * it does not prevent the collision risk, but only minimizes it. **How to deal with local variable assignment without any collision risk?** Substitution rules ~~and Modules~~ have their own problems with recursive functions, which is why I use Blocks, but it's maybe the wrong tool. **Edit** : There is no good reason not to use Module, except if one wants to locally modify the recursion limit, it seems. For my part, I was just confused, and **I should just use Module**.


title: How can I find a particular character in an image that has many similar Chinese characters?

text: The Chinese characters "竟" and "竞" are very similar in their appearance. I am trying to find "竟" in the following picture. I tried use `TextRecognize` but it didn't work. Can you tell me any other way? ![enter image description here](http://i.stack.imgur.com/GR1Ml.jpg) ![enter image description here](http://i.stack.imgur.com/xlscr.jpg)


title: Mathematica function/package for shuffle permutations

text: Does anyone knows a way to compute a list of all (p,q)-shuffles in mathematica? For a definition of the shuffle permutations see for example http://ncatlab.org/nlab/show/shuffle I'm dreaming of a function with p and q as arguments that gives me: 1.) The number of all (p,q)-shuffles. 2.) A list of the actual "shuffle" (p+q)-permutations 3.) The sign of each such permutation


title: Visualizing several long lists of numerical information to see relative frequency

text: _Mathematica_ has very powerful visualization techniques. However, I'm at a loss at how to best make the following chart readable.               BarChart3D[{viewerCount1, viewerCount2, viewerCount3},       ChartLabels -> {Placed[{word1, word2}, None],         Placed[filelist, None]}, ChartLayout -> "Grid",       BarSpacing -> {0.5, 0},       LabelingFunction -> (Tooltip[Row[Flatten[{#3, #1}], " - "]] &)]      ![enter image description here](http://i.stack.imgur.com/STnJH.png) Each row (i.e. each viewerCount) is the number of times a given word appears across a large corpus of files (so if `viewerCount1` was the count for the word 'coffee,' the first row would show what files are most relevant for that word at a glance). The goal is that very quickly, an archivist could see what files are best without having to parse textual data, etc. The preprocessing of all this data has taken place in _Mathematica_ so I would like to keep visualizations in there too (as opposed to learning Processing or another language). There are 385 files being searched here, so each of those is a long list of numbers. I have terrible handwriting, I'm afraid, but here is what I have been trying to create: ![enter image description here](http://i.stack.imgur.com/nVdOB.png) What's different in my 'dream' rather than what I have?   * Bars are more distinctive, not as flimsy   * If the value is zero, still see a little flag rectangle that can draw user to tooltip it   * spacing between the different rows, just to visually distinguish them I suspect `BarChart3D` isn't the best command for this, but I've struggled to get other visualization techniques working. Alternatively, perhaps the problem with this visualization lies before the visualization stage itself... **An Appendix for convenience:**               viewerCount1={0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 1, 0, 0, 0, 0, 0, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 10, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 1, 8, 0, 2, 0, 0, 0, 4, 0, 1, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 9, 0, 0, 2, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10,      0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 3, 2, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 29, 0, 3, 0, 0, 0, 0, 0, 0,      0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,      0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10, 0, 0, 1, 0, 0, 0, 0,      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0}


title: plotting the solutions of system of differential equations

text: The following code include the system of differential equations:                   Derivative[1][x][t] == -10^12 x[t] - y[t]+2*10^-10,          Derivative[1][y][t] == x[t] - 10^12 y[t] + z[t]-2,          Derivative[1][z][t] == - 2*y[t],          x[0] == y[0] == 0, z[0] == 1      I want to plot z[t] vs t. But, unfortunately Mathematics cannot plot it. How can I plot it?


title: DSolve to solve a system of differential equations

text: I need to solve a system of differential equations as follows:               Clear[t]          system =       { x'[t] ==  2.2758*y[t] - x[t]/200,          y'[t] == - 2*2.2758*E^(-(0.7071^2))*E^(-((1.77*t)^2/13.8^2))*Cos[-0.5*t]*z[t]                  + 2.2758*x[t]  - y[t]/200,          z'[t] == 2*2.2758*E^(-(0.7071^2))*E^(-((1.77*t)^2/13.8^2))*Cos[-0.5*t]*y[t]};     initialvalues = {x[0] ==  0, y[0] ==  0, z[0] ==  -1};          sol = DSolve[ Join[ system, initialvalues], {x, y, z}, t]      Unfortunately, _Mathematica_ is not able to yield a result. Is there a way to do it? Also I have another question. I want to make the values of `x` , `y` and `z` be numbers (numerical values) at `t == 0`. How I can do that ?


title: Why are my FrameTicks not working on ArrayPlot?

text: For some reason the ticks don't work on this:               ArrayPlot[Thread@{Range[0, 1, 0.01]}, AspectRatio -> 4,       ColorFunction -> "Rainbow",      PlotRange -> {{0.`, 1.`}, {0.`, 1.`}},      FrameTicks -> {{None, {{0, 75}, {1/5, 100}, {2/5, 125}, {3/5,            150}, {4/5, 175}, {1, 200}}}, {None, None}},      Frame -> True,      PlotRangePadding -> 0,      BaseStyle ->        Directive[Opacity[1], FontFamily -> "Bitstream Charter", 10],      ImageSize -> {Automatic, 100}]


title: How can I set the axis for ScalingFunction in a BarChart?

text: Why is the scaling function being applied to the wrong axis? Here, the log ticks are on the Y axis, but they should be on the X axis:               BarChart[{0.1, 1, 10, 100}, ChartElementFunction -> "GlassRectangle",       ChartStyle -> 45, ScalingFunctions -> "Log", BarOrigin -> Left,       Frame -> True]


title: Simplify expression after substitution

text: I tried to substitue variable x of one equation with another variable `z` (`z=nx/(nx+m)`). See as follow:               Solve[z == n*x/(n*x + m), x]      > >     x -> -((m z)/(n (-1 + z))) >   Then:               PDF[FRatioDistribution[n, m], x] /. x -> -((m z)/(n (-1 + z)))      > >      (m^(m/2) n^(n/2) (-((m z)/(n (-1 + z))))^(-1 + n/2) (m - (m z)/(-1 + > z))^(1/2 >       (-m - n)))/Beta[n/2, m/2] >   How can I collect the term `z` and `1-z` together? Factually, if collecting the results by ``z and`1-z`, it can be very simple (that is Beta distribution). Any suggestion will be greatly appreciated.


title: "Masking" an InputField

text: In building an assignment for my class, I would like set up an InputField, where they will enter a function. I was wondering if there was a way to "continuously" check if their function is being defined correctly (the thing is, they often forget the underscore after the variables as in f[x_], they use ( ) instead of [ ], etc.). I am imagining some red "X" that becomes a green checkmark when the beginning of the definition is correctly written (similar to the way some websites show you if you've re-entered your password correctly when signing up for an account). If this seems too complicated, I would even be interested in seeing how to mask something simpler, like perhaps a phone number. Is this possible in Mathematica? Thanks!


title: Plotting results of numerical differential equation

text: I'm numerically solving the ODE equation $a = -kv$, i.e. $ x'' = -kx'$.               sol = NDSolve[{y''[x] + 0.1 * y'[x] == 0, y[0] == 0, y'[0] == 10}, y, {x, 0, 100}]      I need to do plots of `x` versus time, `v` versus time, and `v` versus `x`. I figured out how to do the first two plots, but I'm having problems with plotting `v` versus `x` (i.e. `y'` versus `y`)


title: Context "Unique to This Notebook" makes variables black even if they're not defined?

text: I'm trying to use this technique to keep variables separated between different notebooks. It seems to be working fine, except that **variables turn black as soon as they are evaluated** — any ideas why this would be happening? Here are some examples: ![](http://i.stack.imgur.com/tVLX1.png) And with a variable that hasn't been used at all: ![enter image description here](http://i.stack.imgur.com/M5ltH.png) Suddenly it's black, even though it has no value: ![](http://i.stack.imgur.com/3EyDk.png)


title: World plot without borders between countries

text: I would like to plot the world map without the countries' borders and just the coastlines (in _Mathematica_ 8). I need to overlap a temperature contour plot on the world map. When I plot country data with lines, it is so crowded because of the lines related to the borders. Plotting with polygons does not help, since I cannot overlap the second plot on top of it. For example in the following command how can I add an option to exclude the borders?                   Graphics[CountryData[#, "Polygon"] & /@ CountryData["Continents"],                   ImageSize -> 400] /. Polygon -> Line      In other words, I need something nice like this but not with polygons and only the coastlines:                   map = Graphics[{GrayLevel[0.5], AbsoluteThickness[1.6],         CountryData[#, {"FullPolygon", {"Equirectangular", {0,             180}}}]} & /@ {"Antarctica", "World"} /. Polygon -> Line,         Frame -> True, FrameStyle -> Directive[Thick],         PlotLabel -> "Precipitation",         FrameLabel -> {"Longitude (°E)", "Latitude (°N)"},        Axes -> True, AxesStyle -> Directive[Gray],        LabelStyle -> Directive[Black, 18, Bold], PlotRangePadding -> None,       ImageSize -> 600,        GridLinesStyle -> Directive[AbsoluteThickness[2], GrayLevel[0.25]],       GridLines -> {Automatic, Automatic},        Method -> {"GridLinesInFront" -> True}];


title: Dealing with a huge dataset

text: I have two ingredients here:   * a big dataset contained in a list, with ~ 20M values.    * a function that takes as each element of the list as input and yields True or False I want to save somewhere the elements of the list that yielded True. Usually I would do something like that:               list = Range[10];     fun = PrimeQ;     Reap[       If[fun[#]&,Sow[#]]& /@ list     ]      This works perfectly, except the fact that we have to wait the end of the computation in order to be able to see all the results. When dealing with such huge lists sometimes waiting for the computation to finish is not an option. This is what I am doing now to split the computation is more chunks,               SaveResult[list_, partitions_, fun_] := Table[       Print["doing iteration ", i];       If[             #[[2]] =!= {},             #[[2]] >>> NotebookDirectory[] <> "/Data/" <> ToString[i]            ] &@           Reap[(If[fun[#], Sow[#]]) & /@ Partition[list, partitions][[i]]],        {i, 1, 9}       ];      My question is: is there a better way to saving partial results or dealing with huge datasets?


title: Manipulating formal vectors

text: I want to work with a formal vector space spanned by integer triplets: $$ V=\text{span}\\{(a,b,c):a,b,c\in\mathbb N^3\\}, $$ such that Mathematica can do additions for me. For example,if I input {1,0,1} + {2,0,1} + 3{1,0,1} - {4,4,7} - 3{2,0,1} Instead of getting {-4, -4, -5} where it is calculated entrywise, I would like Mathematica to return 4{1,0,1} - 2{2,0,1} - {4,4,7} What is the simplest way to implement this?


title: discrete Fourier transform

text: As an answer to my code I will get a periodic discrete time dependent function called "data" which I want to get a Discrete Fourier transform of it using just one period of it,but I think some thing is wrong with the way I do Fourier transform.here is my code:               a = 0.05; L = 15; T = 20 \[Pi]; hdc = 2.1;     sol = NDSolve[{a*D[u[t, x], t] == D[u[t, x], x, x] - Sin[u[t, x]],      u[0, x] == 0, Derivative[0, 1][u][t, 0] == Tanh[t/0.01]*hdc,      Derivative[0, 1][u][t, L] == 0}, u, {t, 0, T}, {x, 0, L},      MaxStepSize -> 0.005, MaxSteps -> 10^6];     q := NIntegrate[(      Evaluate[First[Derivative[1, 0][u][tp, x] /. sol]])^2, {x, 0, L},     Method -> "LocalAdaptive", MinRecursion -> 50,      MaxRecursion -> 100];     data = Parallelize[Table[q, {tp, 25.8, T, 0.5}]];     ListLinePlot[data, PlotRange -> All]     ListLinePlot[Abs[Fourier[data]]]      ![enter image description here](http://i.stack.imgur.com/DWCoI.png) ![enter image description here](http://i.stack.imgur.com/GzxPP.png)


title: How can I make a "self-referential" replacement operation self-contained?

text: Fairly often I find use for replacement rules that call themselves on the right-hand side of the rule, e.g.:               rule = {p___, a, b, c, q___} :> Join[{p, "abc"}, {q} /. rule];      Unfortunately this makes the rule dependent on the specific definition of the Symbol `rule`; it can no longer be copied and modified as an independent expression. Ideally I would like to use the `#0` syntax for `Function` to call the replacement but the most direct application does not work:               # /. {p___, a, b, c, q___} :> Join[{p, "abc"}, #0 @ {q}] & @ "x"      > RuleDelayed::rhs: Pattern p___ appears on the right-hand side of rule > {p___,a,b,c,q___}:>Join[{p,abc},(#1/. > {p___,a,b,c,q___}:>Join[{p,abc},#0[{<<1>>}]]&)[{q}]]. >> Due to the specific mechanism of `#0` its substitution cannot be delayed using e.g. `Slot @@ {0}`; the substitution will never occur. The only pseudo-solution I have found is to move `#0` out of the RHS of the rule, then substitute it using another replacement:               # /. {p___, a, b, c, q___} :> "Join"[{p, "abc"}, "#0" @ {q}] /. {"Join" -> Join, "#0" -> #0} & @      {1, 2, a, b, c, 3, {4, a, b, c}, 5}      > >      {1, 2, "abc", 3, {4, "abc"}, 5} >   However this destroys the purpose of my construct which was to avoid repeatedly scanning from the beginning of the expression. (There are other uses for self-referential replacement but I believe they are all seriously degraded by this "solution.") It is also fragile to use replacement of `"Join"` and `"#0"` so additional complexity or dependency would be required. Is there some other way I can shield the replacement rule (substituted for `#0`) from itself without changing the replacement traversal of the expression? Is there an approach not involving `#0` to accomplish my goal?


title: How do I numerically solve a custom function?

text: Whenever I use functions like `FindRoot` or `NDSolve`, it sends `x` through the function and deals with the result. That would be fine if I was sending a simple math function through, but I have something more like a short program. Is there any way to make it solve it by putting specific values through the function and looking at the results?


title: The number of real values assumed by combinations of complex roots of a decic

text: Define _complex conjugate_ as the pair of complex numbers $a+bi,a-bi$. Assume you have 5 such pairs and they are the 10 roots $x_i$ of a 10th-deg equation with integer coefficients. _In general_ , how many combinations, $$F= x_1 x_2+x_3 x_4+x_5 x_6+x_7 x_8+x_9 x_{10}$$ _**will be real_** , and how do we tell Mathematica to enumerate them all?


title: Using Goto and Label in compiled function?

text: Below I post some rudimentary code to show that there is something wrong using Goto and Label commands in compiled functions (or I believe there is something wrong). It seems that the command Label is not found. Does anybody know what's going on?               toy=Compile[{{x,_Real}},       Label["t"];       y=RandomReal[{-1.,1.}];       If[x*y<-0.,Print[y*x];Goto["t"],Print[y*x," ","done"]]     ];          toy[1.]     (*Goto::nolabel: "Label ("t") not found.*)


title: Write equation as a multinomial in the variables x and Sin[y] and extract the coeficient of x*Sin[y]^2

text: I need help solving this problem. Thank you. (a*x^2 + b*x*Sin[y] + c*Sin[y])^2 + (a*Sin[y]^2 + b*x)^3 Write the above as a multinomial in the variables x and Sin[y] and extract the coeficient of x*Sin[y]^2. Ans : 2bc


title: About calculating Integrals

text: I am trying to Integrate the following Integral (all of the variables are reals) $$\int Exp(-(\cos^{-1}\left[\frac{\text{n}_0 (\text{v}_0-\text{x}_0)+\text{n}_1 (\text{v}_1-\text{x}_1)+\text{n}_2 (u \text{v}+\text{v}_2-\text{x}_2)}{\sqrt{\text{n}_0^2+\text{n}_1^2+\text{n}_2^2} \sqrt{(\text{v}_0-\text{x}_0)^2+(\text{v}_1-\text{x}_1)^2+(u \text{v}+\text{v}_2-\text{x}_2)^2}}\right]/m)^2) \, du$$ The formula, simplified, in Mathematica input:               Integrate[Exp[-(ArcCos[(a + b*u)/(c*Sqrt[d + (e + b*u)^2])]/m)^2], u]      If I simply enter that into Mathematica, it instantly returns the same expression. If I enter it as an Integral with set limits (0 and t), it calculates on and on and on.... How do I go about this? Is there any tricks that can be applied to   1. Find out if there is a closed Form for this (obviously including erf terms)   2. Get an Approximation for the Integral (For example using series)   3. Make Mathematica take more time trying to find a closed form for the indefinite Integral/antiderivative?


title: An apparently "simple" limit?

text: Let $c$ and $h$ be real values. I was using _Mathematica_ to compute the limit $(h \rightarrow \infty)$ of the following expression: $$ \frac{\left(h^2 +c^2 h^2 + \sqrt{4 h^2+\left(h^2+c^2 h^2\right)^2}\right)^2}{4\left(h^2 + \frac{1}{4}\left(h^2 +c^2 h^2 + \sqrt{4 h^2+\left(h^2+c^2 h^2\right)^2}\right)^2\right)}\text{//FullSimplify} $$ which is equal to $$=\frac{2}{4+(1+c^2)^2 h^2 - (1+c^2)\sqrt{4h^2+(1+c^2)^2 h^4}}$$ Just as the rough estimate anyone can make by counting powers of $h$ in the numerator and denominator predicts, the _Mathematica_ function `Limit[f[h], h -> ∞]` gives $1$ in the nonsimplified case and gives $0$ in the simplified case. Now, I am a little confused about whether the correct limit is $1$ or $0$ ? And is the function `Limit` supposed to work that way as to not being able to distinguish the correct result from the wrong one ? Here is the _Mathematica_ code               a = (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2/(        4 (h^2 +            1/4 (h^2 + c^2 h^2 + Sqrt[4 h^2 + (h^2 + c^2 h^2)^2])^2)) //         FullSimplify;          Limit[a, h -> Infinity]      as you know, you can add or remove the `//FullSimplify` at the end of the definition for $a$ to switch between the two cases. **UPDATE:** As mentioned in the comments below, I have reported this as a bug to Wolfram support. They wrote back, were very polite and thankful, mentioned that this issue has been forwarded to their developer team and recommended to use `Assumptions` in order to be sure that everything goes well and additionally even speed up the computation.


title: Vertical RadioButtonBar in Manipulate

text: I want to make a vertical `RadioButtonBar` in `Manipulate`. So far I only have a workaround that looks like this:               Manipulate[       Which[       sel == 0,       Text["one"],            sel == 1,       Text["two"],            sel == 2,       Text["three"]       ],       Control[{{sel, 0, ""}, {0 -> "one"}, ControlType -> RadioButton}],       Control[{{sel, 1, ""}, {1 -> "two"}, ControlType -> RadioButton}],       Control[{{sel, 2, ""}, {2 -> "three"},        ControlType -> RadioButton}],            ControlPlacement -> Left       ]      Is there any way not to have a line for each of the radiobuttons? Thanks!


title: Using DynamicModule variables outside the DynamicModule

text: Take the following example for hiding/showing a notebook through a `Checkbox`:               x = True;     CreateDocument[{"hello"}, Visible -> Dynamic[x]];     Checkbox[Dynamic[x]]      It works correctly until the kernel is restarted, but not after. The documentation explains that the variables in a `DynamicModule` are owned by the front end, not the kernel. So to make this piece of GUI resistant to kernel quits/crashes, I tried:               DynamicModule[      {x = True},      CreateDocument[{"hello"}, Visible -> Dynamic[x]];      Checkbox[Dynamic[x]]     ]      This one does not work though. Can someone explain why, and whether there is to use this `DynamicModule` variable outside the `DynamicModule`? Is this behaviour somehow related to the `DynamicModuleParent` option? * * * **Update:** As Brett Champion has pointed out in the comments below, there's a similar example in the documentation in the DynamicModule Wormholes section. This shows how to make it possible to access a particular `DynamicModule`'s variable from another `DynamicModule` using the `InheritScope -> True` setting. However, my example I don't have a second `DynamicModule` in the newly created notebook, only a single `Dynamic`, so I am not sure how to connect the variables.


title: Outputting a fraction on one line with an oblique division sign

text: How can I make _Mathematica_ output a fraction as `x/2` as opposed to the `x` on top with a horizontal bar and `2` on the bottom? Perhaps I need to work more with `Riffle` but so far I cannot get it. I have tried concatenating strings, but _Mathematica_ reverses the order of the two elements. Specifically, I am trying to put the symbol `theta` with a subscript `i` and then a division by `2`. It is going onto a figure, and is very difficult to read with the vertical division. My code is:               Do[       alabel[i] =         Text[StyleForm[Subscript["\[Theta]", i]/2, FontSize -> 10,                        FontColor -> Black, FontFamily -> "Times"],              {Re[ww[tt[[3]]]], Im[ww[tt[[3]]]]} + .2 {Cos[aloc3[[i]]], Sin[aloc3[[i]]]},              ioffset[i]]       , {i, 1, 6}];          th3 = Graphics[Table[alabel[j], {j, 1, 6}]];      (The coordinate `ww[[tt[[33]]]]` is just some point in the plane).


title: Will Version 9 functions all work in Version 10?

text: Having developed a substantial code base in _Mathematica_ 9, leads me to ask these important questions prior to downloading and installing _Mathematica_ 10. Will any _Mathematica_ functions from the previous version either not work in _Mathematica_ 10 or work differently? If so, does anyone know of a reference to such changes? I've looked for one but only find the new features.


title: Why the inequality does not take into account the domain?

text: I have this inequality:               Reduce[(4000-1000k)/(k-4) < 0]      and the answer is               k ∈ Reals      I would expect `k != 4`.


title: Finding differences between Pi with varying number of decimals

text: I have the following code               In[32]:= N[Pi, 2]          Out[32]= 3.1          In[33]:= N[Pi, 1]          Out[33]= 3.          In[34]:= N[Pi, 2] - N[Pi, 1]          Out[34]= 0.*10^-1      Why can't Mathematica find the difference between the two numbers?


title: How to find a specific error message?

text: Is there a way to find an existing error message in the `System`` context, that is used for a certain case? I want to include a message in a custom function that compares all the arguments of the function, and checks if they have the same length. I'm pretty sure that there is one such, but I can't seem to find it. And I don't want to insert new symbols or messages to the namespace if I can use existing ones. Now I am looking for a method that is general; i.e., it can be applied to other cases, not just this one with argument lengths. For example I am looking for a message that works like this:               General::badlen = "The length of arguments `1`, `2`, `3` must be the same.";          myFun[x_List, y_List, z_List] /;        If[Length@x == Length@y == Length@z, True,          Message[myFun::badlen, x, y, z]; False] := x*y*z;


title: How to add an interpolating point to InterpolatingFunction?

text: Suppose I have an `InterpolatingFunction` object defined on $[\epsilon,1]$; let's call it `f`. What is the best way to construct a function `g` that coincides with `f` on $[\epsilon,1]$ and satisfies `g[0]==a` (`g` should be an `InterpolatingFunction`, too)? Small perturbations near $\epsilon$ are allowed. Any generalizations are welcome as well.


title: How to modify a matrix to satisfy a special condition?

text: I have matrix like this: ![enter image description here](http://i.stack.imgur.com/FZB6A.png) How do I modify this matrix to make it satisfy the following condition: > For each element _{i, j}_ in the matrix the sum of the elements of row _i_ > must be equal to the sum of the elements of column _j_. Zeros can't be > modified. ### Examples Good: ![good](http://i.stack.imgur.com/DuyAJ.png) Bad: ![bad](http://i.stack.imgur.com/mxTnm.png)


title: Launching remote slave kernels on HPC causes slow down

text: I had a post earlier regarding how I can run Mathematica at a HPC cluster (running `PBSpro`) by having a single `MathKernel` launch say 190 slaves remotely, in order so that I only utilise one `MathKernel` license. I now have this working and the schematic of my notebook looks like:               (*First launch all the remote kernels*)               Needs["SubKernels`RemoteKernels`"];     numCore=10; (*change to num cores per node*)     PBSKernel[host_String]:=     LaunchKernels[RemoteMachine[host,"ssh -x -f -l `3` `1` /path/to/math -mathlink -linkmode Connect `4` -linkname '`2`' -subkernel -noinit",numCore]];     hostfile = Environment["PBS_NODEFILE"];          If[hostfile =!= $Failed, (* if we're running in a batch job *)                 hosts = Import[hostfile, "List"];                 $ConfiguredKernels = Join[$ConfiguredKernels, PBSKernel /@ hosts];     ]          Export["kernels.txt", ParallelTable[{$MachineName, $KernelID}, {$KernelCount}], "Table"];          (*now they are launched evaluate something over them, first creating some output files,      one for each node, postfixed by the node name*)         $runningLogFile = "log.data";         ParallelEvaluate[If[FileExistsQ[$runningLogFile <> $MachineName],   Get[$runningLogFile <> $MachineName],Export[$runningLogFile <> $MachineName, "", "Text"];];]         DistributeDefinition[...];         ParallelMap[func, Tuples[{Range[1/10, 1, 1/10], Range[0, 59]}]];      where               func[{x_?NumericQ, y_?IntegerQ}] :=      Block[{},      g[r_,s_,t_]:= (*do calc, essentially an integral,altho it does depend on some prev defined helper functions external*);           (*func to write this to file in certain format*)      gwrite[r_, s_,t_] := gwrite[r,s,t] = g[r, s, t] /.  v_:>(PutAppend[Unevaluated[g[r, s, t] = v;], $runningLogFile <> $MachineName];  v);     (*write to file using CriticalSection to stop clashing*)     CriticalSection[{logfileLock}, Do[gwrite[tf, x, y], {tf, taufStart, taufEnd, intv}]];      ]      Now this all works fine, in the sense of giving me the expected data, but I have noticed a massive slowdown in computation speed this way, versus the previous method I was using whereby I launched 19 `MathKernels`, one for each compute node and then in turn had each launch 10 slaves locally. This is despite the fact that I have been assigned exactly the same resource, i.e. 19 nodes each with 10 cores and the same memory. The only difference is I am now only using one `MathKernel` instead of 19 to launch the slaves. The old method did my test portion of paramater space in 20mins vs over 4hrs for the current method, which is a drastic slow down that is too much for my longer jobs. What could be causing this slow down? Ideas: 1)Remote slave kernels need to communicate and this is causing some bottleneck. I can't see why this would be a problem for me as my `func` above is just an integral that should not depend on anything the other kernels are doing. 2) Somehow all the remote cores are not getting used? I know they are being Launched as they appear in the `Kernels.txt` file at the start of the calculation and also `$KernelCount=190`, but is my method of employing `ParallelMap` successfully using them all? Is there a way I can check this? I also know all the nodes are being used as I append the `$MachineName` to each data file and get 19 of them as expected. 3) Anything else? Mathlink just can't handle this many cores successfully? If I keep exactly the same script and clone it 19 times, then split the list of `Tuples` equally between the 19 new scripts, and feed these into 19 nodes with a `MathKernel` running on each (which launches 10 local cores) in serial then things massively speed up and I am back to 20mins. EDIT: after reviewing the output logs there is one odd message, although I'm not sure how relevant it is:               StringForm["From `1`:",  Parallel`Kernels`kernel[Parallel`Kernels`Private`bk[remoteKernel[SubKernels`RemoteKernels`Private`lk[LinkObject["49195@10.141.0.3,46878@10.141.0.3", 97, 97], "node012", {"node012", "ssh -x -f -l `3` `1` math -mathlink -linkmode Connect `4` -linkname '`2`' -subkernel -noinit"}, SubKernels`RemoteKernels`Private`speed$1327]], Parallel`Kernels`Private`id$1376, Parallel`Kernels`Private`name$1376], Parallel`Kernels`Private`ek[Parallel`Kernels`Private`nev$1377, Parallel`Kernels`Private`pb$1377, Parallel`Kernels`Private`rd$1377], Parallel`Kernels`Private`sk[Parallel`Kernels`Private`q$1378, Parallel`Kernels`Private`n0$1378, Parallel`Kernels`Private`n1$1378]]]      "DeleteFile::nffil: File not found during DeleteFile[/panfs/panasas01.panfs.cluster/username/out.datanode012]."


title: Simplification of double symbolic sums containing a DiscreteDelta without explicit summation range

text: I am trying to get Mathematica to automatically do simplifications like the following: $$\sum\limits_{q}^{q\in qV}\sum\limits_{q'}^{q'\in q'V}{f(q)g(q')\delta(q-q')}=\sum_{q}^{q\in qV}{f(q)g(q)}.$$ Where the range of values that $q$ and $q'$ are denoted by $qV$ and $q'V$ respectively and are the same. With pen an paper it is really easy to do that kind of simplifications because I know that both variables $q$ and $q'$ have the same range and that the Discrete Delta function $\delta(q-q')$ kills one of the sums. I have solved it, but it seems to me that I haven't done it in the most elegant/efficient way. My solution is:               Clear[qV, qpV, f,g, MyDiscreteDelta]     a = Sum[f[q] Sum[g[qp] MyDiscreteDelta[q, qp], {qp, qpV}], {q, qV}]     b = % /. {qpV -> {q, q1, q2, q3, q4}}     c = % /. {MyDiscreteDelta[x_, y_] -> HoldForm[If[x =!= y, 0, 1]]}     d = ReleaseHold[%]      $a$ defines the sum and I deliberately leave the sum range undefined as even though I know the real range of number I will be summing, it's too large to be put directly and more importantly, at this stage of the calculation I don't care about it, I am just interested in simplifying the result as much as possible, i.e., get rid of one of the sums by using the delta function. I am not using Mathematica defined function DiscreteDelta because a) it takes too long and b) when used symbolically it is left unevaluated. First thing I need to do to get rid of one of the sums is to give it a range to sum through, in doing so, I am actually giving it a list of possible variables that it could take(assuming that in the general case I will be having more than two summations and I will have more than one delta function in many of the variables). That's done in $b$. $c$ changes the empty definition of MyDiscreteDelta but I need to hold it as if I don't do so, it will all be evaluated to $0$. The last step, $d$ evaluates the expression and gives the results I want. This works, but I was wondering if there is a simpler way of doing it. Thanks in advance.


title: Polar Plot of a two variable function

text: I am trying to make a polar plot of a two variable function, which has a radial and an angular part. For instance `Y[r,φ]=r^2 Cos[φ]`. I can make this plot using `ContourPlot`, but is there a way to do it using `PolarPLot`? My `ContourPlot` is               Y[r_, f_] := r^2 Cos[f]     ContourPlot[Y[Sqrt[x^2 + y^2], ArcTan[x, y]], {x, 0, 1}(*Radius till 1*), {y, 0, 60}(*Angle from 0 to 60*)]      ![](http://i.imgur.com/h9I9tJs.png)


title: How to get the matrix coefficient from a QuadraticForm expression in a very fast way?

text: As we know, `QuadraticFormDistribution` is the distribution of `z.a.z + b.z + c` for a real-valued p x p symmetric positive definite matrix a, a length-p vector b, a scalar c, and a p-dimensional multivariate normal vector z I have a `QuadraticFormDistribution` expression `fexp`, which has a very complex form, and it can be simplified to the form `z.a.z + b.z + c`. Now, I want to get the matrix a, the vector b, and the scalar c from `fexp`. I think there are many methods for this simple problem. I always `ExpandAll` the expression `fexp` and use the `Coefficient` to get the matrix a, the vector b and the scalar c. The problem is: I will deal with very large numbers of these kinds of quadratic expressions, so I want to find a faster method to deal with this problem. Take an example,               fexp=12 - 9 z[3] - 16 z[4] - 25 z[5] - 36 z[6] - 49 z[7] - 64 z[8] -      81 z[9] - 100 z[10] - 121 z[11] - 144 z[12] +      k^2 (-1 + 2 z[1] + 2 z[2] + 2 z[3] + 2 z[4] + 2 z[5] + 2 z[6] +      2 z[7] + 2 z[8] + 2 z[9] + 2 z[10] + 2 z[11] + 2 z[12]) +      z[1] (-1 + 54 a[1, 2] z[2] + 96 Sqrt[2] a[1, 3] z[3] +      150 Sqrt[3] a[1, 4] z[4] + 432 a[1, 5] z[5] +      294 Sqrt[5] a[1, 6] z[6] + 384 Sqrt[6] a[1, 7] z[7] +      486 Sqrt[7] a[1, 8] z[8] + 1200 Sqrt[2] a[1, 9] z[9] +      2178 a[1, 10] z[10] + 864 Sqrt[10] a[1, 11] z[11] +      1014 Sqrt[11] a[1, 12] z[12]) +      k (1 + 4 a[1, 1] z[1]^2 + 3 z[3] + 4 z[4] + 5 z[5] + 6 z[6] +      7 z[7] + 8 z[8] + 9 z[9] + 10 z[10] + 11 z[11] + 12 z[12] +      z[1] (1 + 18 a[1, 2] z[2] + 32 a[1, 3] z[3] + 50 a[1, 4] z[4] +         72 a[1, 5] z[5] + 98 a[1, 6] z[6] + 128 a[1, 7] z[7] +         162 a[1, 8] z[8] + 200 a[1, 9] z[9] + 242 a[1, 10] z[10] +         288 a[1, 11] z[11] + 338 a[1, 12] z[12]) +        2 (8 a[2, 2] z[2]^2 + 18 a[3, 3] z[3]^2 + 49 a[3, 4] z[3] z[4] +         32 a[4, 4] z[4]^2 + 64 a[3, 5] z[3] z[5] +         81 a[4, 5] z[4] z[5] + 50 a[5, 5] z[5]^2 +         81 a[3, 6] z[3] z[6] + 100 a[4, 6] z[4] z[6] +         121 a[5, 6] z[5] z[6] + 72 a[6, 6] z[6]^2 +         100 a[3, 7] z[3] z[7] + 121 a[4, 7] z[4] z[7] +         144 a[5, 7] z[5] z[7] + 169 a[6, 7] z[6] z[7] +         98 a[7, 7] z[7]^2 + 121 a[3, 8] z[3] z[8] +         144 a[4, 8] z[4] z[8] + 169 a[5, 8] z[5] z[8] +         196 a[6, 8] z[6] z[8] + 225 a[7, 8] z[7] z[8] +         128 a[8, 8] z[8]^2 + 144 a[3, 9] z[3] z[9] +         169 a[4, 9] z[4] z[9] + 196 a[5, 9] z[5] z[9] +         225 a[6, 9] z[6] z[9] + 256 a[7, 9] z[7] z[9] +         289 a[8, 9] z[8] z[9] + 162 a[9, 9] z[9]^2 +         169 a[3, 10] z[3] z[10] + 196 a[4, 10] z[4] z[10] +         225 a[5, 10] z[5] z[10] + 256 a[6, 10] z[6] z[10] +         289 a[7, 10] z[7] z[10] + 324 a[8, 10] z[8] z[10] +         361 a[9, 10] z[9] z[10] + 200 a[10, 10] z[10]^2 +         196 a[3, 11] z[3] z[11] + 225 a[4, 11] z[4] z[11] +         256 a[5, 11] z[5] z[11] + 289 a[6, 11] z[6] z[11] +         324 a[7, 11] z[7] z[11] + 361 a[8, 11] z[8] z[11] +         400 a[9, 11] z[9] z[11] + 441 a[10, 11] z[10] z[11] +         242 a[11, 11] z[          11]^2 + (225 a[3, 12] z[3] + 256 a[4, 12] z[4] +            289 a[5, 12] z[5] + 324 a[6, 12] z[6] + 361 a[7, 12] z[7] +            400 a[8, 12] z[8] + 441 a[9, 12] z[9] +            484 a[10, 12] z[10] + 529 a[11, 12] z[11]) z[12] +         288 a[12, 12] z[12]^2 +         z[2] (1 + 25 a[2, 3] z[3] + 36 a[2, 4] z[4] +            49 a[2, 5] z[5] + 64 a[2, 6] z[6] + 81 a[2, 7] z[7] +            100 a[2, 8] z[8] + 121 a[2, 9] z[9] + 144 a[2, 10] z[10] +            169 a[2, 11] z[11] + 196 a[2, 12] z[12]))) +         2 (z[2] (-2 + 75 a[2, 3] z[3] + 108 Sqrt[2] a[2, 4] z[4] +         147 Sqrt[3] a[2, 5] z[5] + 384 a[2, 6] z[6] +         243 Sqrt[5] a[2, 7] z[7] + 300 Sqrt[6] a[2, 8] z[8] +         363 Sqrt[7] a[2, 9] z[9] + 864 Sqrt[2] a[2, 10] z[10] +         1521 a[2, 11] z[11] + 588 Sqrt[10] a[2, 12] z[12]) +         3 (49 a[3, 4] z[3] z[4] + 64 Sqrt[2] a[3, 5] z[3] z[5] +         100 Sqrt[2] a[4, 6] z[4] z[6] +         121 Sqrt[3] a[4, 7] z[4] z[7] + 169 a[6, 7] z[6] z[7] +         288 a[4, 8] z[4] z[8] + 196 Sqrt[2] a[6, 8] z[6] z[8] +         225 a[7, 8] z[7] z[8] + 169 Sqrt[5] a[4, 9] z[4] z[9] +         225 Sqrt[3] a[6, 9] z[6] z[9] +         256 Sqrt[2] a[7, 9] z[7] z[9] + 289 a[8, 9] z[8] z[9] +         196 Sqrt[6] a[4, 10] z[4] z[10] + 512 a[6, 10] z[6] z[10] +         289 Sqrt[3] a[7, 10] z[7] z[10] +         324 Sqrt[2] a[8, 10] z[8] z[10] + 361 a[9, 10] z[9] z[10] +         225 Sqrt[7] a[4, 11] z[4] z[11] +         289 Sqrt[5] a[6, 11] z[6] z[11] + 648 a[7, 11] z[7] z[11] +         361 Sqrt[3] a[8, 11] z[8] z[11] +         400 Sqrt[2] a[9, 11] z[9] z[11] + 441 a[10, 11] z[10] z[11] +         512 Sqrt[2] a[4, 12] z[4] z[12] +         324 Sqrt[6] a[6, 12] z[6] z[12] +         361 Sqrt[5] a[7, 12] z[7] z[12] + 800 a[8, 12] z[8] z[12] +         441 Sqrt[3] a[9, 12] z[9] z[12] +         484 Sqrt[2] a[10, 12] z[10] z[12] + 529 a[11, 12] z[11] z[12] +         z[3] (81 Sqrt[3] a[3, 6] z[6] + 200 a[3, 7] z[7] +            121 Sqrt[5] a[3, 8] z[8] + 144 Sqrt[6] a[3, 9] z[9] +            169 Sqrt[7] a[3, 10] z[10] + 392 Sqrt[2] a[3, 11] z[11] +            675 a[3, 12] z[12]) +         z[5] (81 a[4, 5] z[4] + 121 a[5, 6] z[6] +            144 Sqrt[2] a[5, 7] z[7] + 169 Sqrt[3] a[5, 8] z[8] +            392 a[5, 9] z[9] + 225 Sqrt[5] a[5, 10] z[10] +            256 Sqrt[6] a[5, 11] z[11] + 289 Sqrt[7] a[5, 12] z[12])));            vector z={z[1], z[2], z[3], z[4], z[5], z[6], z[7], z[8], z[9], z[10], z[11], z[12]};      The answer is very clear. I just want to know how to get the matrix a, the length-12 vector b, and scalar c in an efficient way.


title: How to add a common color legend to a grid of density plots?

text: I have a grid of density plots, say something like:               GraphicsGrid[Partition[Table[       DensityPlot[Sin[i x] Sin[y/i], {x, -4, 4}, {y, -3, 3}],       {i, 1, 4}], 2]]      I want to show a color legend bar next to this plot. I know that the option `PlotLegends -> Automatic` does this for a _single_ density plot, so if I do               GraphicsGrid[Partition[Table[       DensityPlot[Sin[i x] Sin[y/i], {x, -4, 4}, {y, -3, 3}, PlotLegends -> Automatic],       {i, 1, 4}], 2]]      I get the same color legend bar repeated four times, one for each density plot. But naturally what I want is to show a single bar common to all the density plots in the `GraphicsGrid`. How do I do this?


title: Nice formatting for the results of a dice game simulation

text: One of my _Mathematica_ interests is in problems related to the formatting of output. Getting nice formatting from _Mathematica_ is often frustrating, but occasionally _Mathematica_ will delight me by making a formatting task, which I anticipated would be difficult, easy. I want to share one such recent experience by posing it as a question. > Can you create a 10-row, 5-column table in _Mathematica_ presenting the > results from running a dice game simulation in which a pair of dice is > rolled 50 times? The table should look like this: ![dice.png](http://i.stack.imgur.com/YhzUZ.png) Your solution may show a different sequence of dice rolls. Note that rolls of seven are highlighted in red. That feature must be reproduced. I post this question, not because it is difficult (it isn't), but because I had fun with it myself. I will post my answer in 48 hours, providing nobody posts an equivalent answer beforehand.


title: Combining and separating sums

text: How can I combine or separate sums in _Mathematica_ in the way that `Together` or `Expand` work for rational expressions？ For example, how does one transform from $$\text{Sum}\left[\frac{a}{\sqrt{n!}},\\{n,0,\infty \\}\right]+\text{Sum}\left[\frac{b}{\sqrt{n!}},\\{n,0,\infty \\}\right]$$ to $$\text{Sum}\left[\frac{a+b}{\sqrt{n!}},\\{n,0,\infty \\}\right]$$ and vise versa?


title: Intermediate variables in NDSolve

text: I'm interested in using NDSolve as an integrator for a system of differential equations that looks like the following:               x'[t] == f[x[t], p[t]]     p'[t] == g[x[t], p[t]]      Both f and g are expensive functions to compute, but they both need to perform similar computations, so it is cheaper to write               vec'[t] == newf[vec[t]]      where vec is a list whose elements are `{x,p}`. Now, I have a quantity computed in newf that just depends on x and p at a given time (call it y), and I would like this quantity to be included in the solution. I could compute it separately as               vec'[t] == newf[vec[t]]     y[t] == calcy[vec[t]]      but this would be expensive, as I'd be computing it twice. I'd like to include it in the vector and just include it with newf's output, but as I'm specifying y and not y', this doesn't work. I tried making a list out of everything as               {vec'[t], y[t]} == newf2[vec[t]]      but _Mathematica_ spits back > NDSolve::underdet: There are more dependent variables, {vec[t], y[t]}, than > equations, so the system is underdetermined. Any suggestions? Here is a minimal working example, where I would like to have y return in the solution as well as the vector solution.               f[vec_List] := Module[{y},       y = Norm[vec];       {1, y}]          NDSolve[{vec'[t] == f[vec[t]], vec[0] == {0, 0}},       {vec},       {t, 0, 1}]          Plot[vec[t] /. %, {t, 0, 1}]


title: Optimizing the performance of an algorithm

text: First - a bit of an introduction. If you're only interested in the code, you can skip this section. The following question is drawn from Dennis E. Shasha's _Puzzling Adventures_ , and is listed under the question **Strategic Bullying**. We have a set of entities, each of them with a power represented by a number. A coalition's power will be the sum of its members' powers. Coalition A will destroy Coalition B if A's power is greater than that of B, but A will back down if B can obtain enough new allies such that B's power is greater than or equal to A. When a coalition is beaten, it is destroyed and its wealth (proportional to its power) is distributed among the winners. The winners do not lose any of their power. Each of the entities plays selfishly, striving above all to avoid certain destruction but also to acquire wealth. Based on the description of the question above, a set of entities is stable if there is a strong enough coalition to prevent any war and every member of that coalition has an interest in participating in the coalition. * * * Mathematically speaking, a set of entities _with distinct integer powers_ is:   * **Unstable** if the most powerful entity has more than half the total power of all entities, since it would destroy all other entities without being challenged;   * **Stable** if the most powerful entity has exactly half of the total power of all entities, since it would not be able to destroy the other entities (they have an incentive to form a coalition to stop the most powerful entity;   * and otherwise, **Stable** if and only if ALL subset of entities formed by removing any coalition of entities with less than half the power of all entities is **Unstable**. _(Thanks to Xerxes for pointing out the error in my original formulation)._ What is the most efficient code we can write to check the stability of a given set of entities? Now, I've written a function to check the stability of any given set of entities:               fstable = Which[        2*Max[#1] > Total[#1], 0,         2*Max[#1] == Total[#1], 1,        True, BitXor[Max[fstable /@ Subsets[#1, {Length[#1] - 1}]], 1]        ] &;      **Edit:** _The above code doesn't cover all the cases in the question as it assumed that one only needed to consider the stability of any subset formed by removing one element - look at Xerxes' answer below for a correct and faster implementation._ Examples of implementation:               fstable[{3, 4, 5, 6}]     0          fstable[{3, 4, 5}]     1      However, this function does not run very fast, and in fact even checking the stability of a set of 10 entities takes some time:               Timing[fstable[Range[10]]]     {30.966199, 0}      Is there any way we can improve this, or will we require some additional mathematical insight to allow us to eliminate specific cases? * * * _Note_ : The question the author actually asks is for us to find the largest set of entities with distinct integer powers and with the most powerful entity having a power of 21. To solve the problem by brute force, I would define the function using the idiom `f[x_]:=f[x]` to make use of _Mathematica_ 's ability to store function evaluations. However, the total number of function evaluations to be solved is quite large, so I'm not sure whether doing so will cause the kernel to run out of memory.


title: How to show the intersection points of three functions?

text: How can I show the intersection points of three functions in my `Plot`? I really can't figure it out.               LogPlot[{(4.058*10^31)/M, (1.647*10^60)/M^2 , M}, {M, 0, 2*10^20},              PlotRange -> {10^9, 10^21} , Filling -> 3 -> Axis,              AxesLabel -> {Subscript[M, U], Subscript["Λ", U]}]


title: Solving simple ODE for eigenvalues causes nightmares

text: I had some code which could solve an ODE for the function $u$ and the eigenvalue $E$. Unfortunately, $E$ is a function of another parameter $k$. This meant, that upon solving the equations, I had to carefully specify a valid range of $E$ and $k$ to search for solutions, or else _Mathematica_ would confuse different "bands" of $E(k)$ values. In an attempt to make the code more automatic, I sliced up the regions into thin rectangles in $E-k$ space (as specified by the tables below)               L = 8; Clear[k, e]; Clear[u];          Table[solution[n] =      NDSolve[{-D[u[x, e, k], x, x] + ((x - k)^2 - e) u[x, e, k] == 0.,       u[0, e, k] == 0, Derivative[1, 0, 0][u][0, e, k] == 1},      u, {x, 0, L}, {e, 0.1 (n - 1), .1 n}, {k, 0., 4.}], {n, 1, 40}];          Table[Eofk[n, k_] :=      e /. FindRoot[(u[L, e, k] /. solution[n]) == 0, {e, .1 n}], {n, 1,     40}];          Table[c[n] =      Plot[Eofk[n, k], {k, 0., 4.}, PlotRange -> {{-0.05, 5.}, {0, 12.}},      PlotStyle -> {RGBColor[1, 0, 0]}], {n, 1, 40}];      The first table solves a second order ODE for a function $u[x,e,k]$. The second table finds values of $ $such that $u[L,e,k]=0$ (i.e., solves the eigenvalue problem). The third should plot the various $e$ values found in the second table as a function of $k$. My original code is similar but it is messy because I have carefully chosen $e,k$ ranges for `NDSolve`. This new code does not work at all (try it). Can anyone tell me why it doesn't work and what I should do to fix it?


title: Problem with numerical evaluation of analytically solved integral, solution way off

text: The following command in Version 9.0.1:               N[Integrate[x^50*Sin[x], {x, 0, 1}]]      gives $1.4615\times 10^{48}$ which is way off from the correct solution which is between $0$ and $0.5$. Interestingly the indefinite integral comes out correct and if you increase the precision Mathematica gives an error and the following number: $0.0163$. The same result is given by `NIntegrate`. What exactly is the problem here, is this a known issue and is there a workaround?


title: How can I make labels with units and then position the label to the right of a bar?

text: I want to add units to my numbers to display them as the labels of my chart. I tried this:               sizes = {10718, 881, 395, 387, 375, 309, 297};     sizesStr = ToString[sizes] <> " kB";     BarChart[sizes, BarOrigin -> Left, ChartLabels -> Placed[sizesStr, Right]]      ![enter image description here](http://i.stack.imgur.com/cTir0.png) So instead of each bar displaying "xxx kB" I get the "kB" after the whole list. Should I blame `StringJoin` for not being listable? How do I get the numbers in each of the bars? Related: for the short bars the text will go beyond the vertical axis. Can I get the labels at the right of the bars, instead of inside them?


title: ReadList can't read from file containing high-precision Reals

text: I want to read a file which contains seven columns of real numbers. Mathematica seems to have a problem with the high precision I need. My code:               solFile = OpenRead["sol_Je9"];      data = ReadList[solFile, {Real, Real, Real, Real, Real, Real, Real}];     Close[solFile];      The error messages: > General::digit: Digit at position 1 in -239.916104094628030 is too large to > be used in base 10. ReadList::readn: Invalid real number found when reading > from sol_Je9. How can I fix this?


title: Speeding up code by avoiding repeated evaluation of a function

text: I want to speed up my code and I have two ideas, but I don't know how I can implement them. Here is a little part of my code, which I want to improve:                  R[FK_, f_] :=            FixedPoint[               Integrate[PDF[f, x]*x , {x, 0, id[Q[FK, #], #]}]               *If[CDF[f, id[Q[FK, #], #]] == 0., 0.1, CDF[f,id[Q[FK, #], #]]]               *id[Q[FK, #], #] &          , seedR] ;      1) You can see that the function is evaluated `id[Q[FK, #], #]` 3 times. That is unecassary. Normaly I would use a Module structure, define a new variable `idd = id[Q[FK, #], #]`, like                   R2[FK_, f_] :=            Module[{iddd},            FixedPoint[              iddd = id[Q[FK, #], #] &;              Integrate[PDF[f, x]*x , {x, 0, iddd}]             *If[CDF[f, iddd] == 0., 0.1, CDF[f, iddd]]             *iddd &           ,seedR]];      but this does not work. Do you have any idea? 2) Probably a easy one:                   If[CDF[f, id[Q[FK, #], #]] == 0., 0.1, CDF[f, id[Q[FK, #], #]]]      This involves the funtion `CDF[f, id[Q[FK, #], #]]` two times. One time should be enough. How can I do this? Thanks Peter


title: How do I show the transition probabilities in a graph of a Markov process?

text: How do I get `Graph` to display the transition probabilities for a Markov process as labels on the graph's edges? The information is clearly present in the graph, but only displays when I hover over the edges. Is there a way to get the information to display as edge labels (without going through complex machinations)? For example,               mp = DiscreteMarkovProcess[{1, 0, 0}, ({          {0.6, 0.1, 0.3},          {0.2, 0.7, 0.1},          {0.3, 0.3, 0.4}         })];     Graph[mp]      gives ![](http://i.stack.imgur.com/Opyht.png) with no labeling of edges, even though I can hover over any edge to see the associated transition probability (s a tooltip).


title: What are the standard colors for plots in Mathematica 10?

text: Mathematica 10 release appears to have changed the default styling of plots: the most visible changes are thicker lines and different default colors. Thus, answers to this stackoverflow question are only valid for Mathematica < 10\. For example, plots in this code will not give identical output in Mathematica 10, although they do in version 9:               fns = Table[x^n, {n, 0, 5}];     Plot[fns, {x, -1, 1}, PlotStyle -> ColorData[1, "ColorList"]]     Plot[fns, {x, -1, 1}]      ![enter image description here](http://i.stack.imgur.com/vBMw1.png) So, my question is: what is the new way of getting the default colors to reproduce for own uses?


title: No response from DSolve

text: Using optimal control theory, I am trying to find the optimal paths of d and k such that they will maximize profit as given by lc:               lc =        -35/2 d[t] (d[t] - 6/10) (d[t] - 9/10) - 20k[t] (k[t] - 5/10) (k[t] - 1) +          d[t]k[t] - u[t]^2 - w[t]^2 + md[t]u[t] +  mk[t]w[t]      with `d'[t] = u[t]` and `k'[t] = w[t]`. This gives me the following:               DSolve[{        d[0] == a,        k[0] == b,        md[20] E^(-r 20) == 0,        mk[20] E^(-r 20) == 0,        d'[t] == D[lc, md[t]],        k'[t] == D[lc, mk[t]],        md'[t] == -D[lc, d[t]] + r md[t],        mk'[t] == -D[lc, k[t]] + r mk[t],        D[lc, u[t]] == 0,        D[lc, w[t]] == 0} /. r -> 1/2,      {d[t], k[t], u[t], w[t], md[t], mk[t]}, t]      But this just returns the input. What am I doing wrong? My ultimate goal with this is to be able to obtain `d'[t] (=u[t])` given `d[0]` and `k[0]`. A symbolic solution in terms of `a` and `b` would be best, but giving numerical values for `a` and `b` and just getting a numerical value for `d'[0]` would also suffice.


title: Histograms with pre-counted data

text: I've got some large data sets which have been counted **but not binned** already - essentially, a list of pairs of values (not bins) and counts.* (Or, equivalently, it's been binned into too-small bins.) I want to plot histograms for them. I remember the deprecated version of Histogram from a separate package had a FrequencyData option, but that seems to have disappeared. Is there any built-in way to accomplish this now? (I'd like to still have all the fancy features of Histogram, i.e. I don't want to just rebin the data myself and plot it directly. Notably I'd like to still be able to use Histogram's automatic bin specification, or something like it.) *That is, my data is represented as `{{1, 6}, {2, 4}, {3, 2}, ...}` instead of `{1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, ...}`. (And before anyone suggests just expanding the data to the latter form to pass to Histogram: there are over 100K values, and the total count is over 100M.) Edit: okay, let me be really explicit. The perfect thing would be to be able to take the first representation of the data (`{{1,6}, ...}`), and get _exactly_ what Histogram would have produced had I given it the second version (`{1,1,1,1,...}`), without having to actually expand it to that form. (This includes being able to specify various options and extra arguments to Histogram.) I do not want a bar chart with 100K bars. I do not want to have to decide how many bins to make every time I do this, because I may do it many times with many varieties of data.


title: Binning events with weights

text: I have two vectors, one containing sampled energies and the other one for each of these energies a weight: `e(i)` and `w(i)`. I bin the energies `e(i)`. Now I want to compute the sum of weights `w` in each bin separately. I can bin the energies using `BinLists` and make a plot with `Histogram`. How do I perform the sum over weights in each bin? I would obviously need access to the indices `i` in each bin. Any elegant solution to this problem would be welcome.               e = {0.129454, 0.160294, 0.140456, 0.152359, 0.174006, 0.186969};     w = {2.12373*10^-6, 0.00029488, 5.68648*10^-6, 0.000100244, 2.28126*10^-6, 6.74131*10^-6};


title: Failure of CountryData due to connectivity?

text: Today I started researching the "Four colours theorem" and to do so I need data, regions for example. I tried using the following line:               CountryData["Ukraine", "Regions"]      But the result: > CountryData::dlfail: Internet download of data for CountryData failed. Use > Help > Internet Connectivity... to test or reconfigure internet > connectivity. >>               $Failed      I tried to do this next: ![CMD, Ping source](http://i.stack.imgur.com/4Gu8d.png) For those who don't understand (because the screen language is Russian) the default source does not respond. How I can change the default source? _Mathematica_ version – 9.0.1.0, OS – Windows 7.