title: How to find the roots of an equation which is almost singular everywhere

text: I'm trying to find the roots of the equation as a variable of `k`:               a = 0.81; b = 1.94;     LegendreP[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, 1] == 0      The problem is that this function is almost divergent everywhere except only a few points, and FindRoot and Reduce seems to fail on this               FindRoot[LegendreP[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, 1] == 0, {k, 0}]      > FindRoot::nlnum: The function value {ComplexInfinity} is not a list of > numbers with dimensions {1} at {k} = {1.49012*10^-8}. >>               Reduce[LegendreP[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, 1] == 0, k]      > Reduce::inex: Reduce was unable to solve the system with inexact > coefficients or the system obtained by direct rationalization of inexact > numbers present in the system. Since many of the methods used by Reduce > require exact input, providing Reduce with an exact version of the system > may help. >>               With[{a = 81/100, b = 194/100},       Reduce[LegendreP[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, 1] == 0, k]]      > Reduce::nsmet: This system cannot be solved with the methods available to > Reduce. >> How would I solve this equation? The roots are `k = n*a`, where n is integer. * * * **Edit** The motivation of this problem is to determine the bound states of a one- dimensional potential well $$ -\frac{\hbar^2}{2m}\frac{d^2\psi}{dx^2}+V(x)\psi=-\frac{\hbar^2k^2}{2m}\psi $$ where $$ V(x)=-\frac{\hbar^2b}{m}sech^2(a x) $$ ,$k$ is a parameter. This equation has solution               ψ[x_]:=A*LegendreP[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, Tanh[a x]] +       B*LegendreQ[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, Tanh[a x]]          V[x_] := -((\[HBar]^2 b)/m) Sech[a x]^2     -(\[HBar]^2/(2 m)) ψ''[x] + V[x] ψ[x] == -((\[HBar]^2 k^2)/(2 m)) ψ[x] // FullSimplify     (*True*)      Quantum mechanics requires the wave function is finite, which gives a constrain on parameter `k` and thus the eigen-energies($\sim k^2$) are some discretized values. And I'm trying to find these eigen-energies. Mathematically, this is equivalent to finding the $k's$ that keep associate Legendre polynomials finite in region `0 <= x <= 1`               A*LegendreP[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, x]+B*LegendreQ[(-a + Sqrt[8 b + a^2])/(2 a),  k/a, x]      . **Edit2** For the wave function must be normalizable and square-integratable, here is the explanation in the Griffiths's Introduction to Quantum Mechanics: > For some solutions to the Schrodinger equation the integral is infinite; in > that case no multiplicative factor is going to make it 1. The same goes for > the trivial solution $\Psi = 0$. Such non-normalizable solutions cannot > represent particles, and must be rejected. Physically realizable states > correspond to the square-integrable solutions to Schrodinger's equation. For the case of free particle, he also explained why it's not normalizable: > In the case of the free particle, then, the separable solutions do not > represent physically realizable states. A free particle cannot exist in a > stationary state; or, to put it another way, there is no such thing as a > free particle with a definite energy.


title: Is everything in Mathematica ultimately stored as a rule?

text: My current understanding about Mathematica is that everything, at the lowest level, ends up as replacement rules.   First question: is this true?   Second question: does Mathematica "gloss over" these rules by representing them as functions? ![](http://i.stack.imgur.com/HViEe.png) doesn't seem to be expressed as a rule. Third question: does the by value/by reference question have any when you "pass" a value to a Mathematica function?


title: FindMinimum [] methods

text: Which algorithm is used when we do not use the option Method of `FindMinimum[]` function.               FindMinimum[f,{x,x0}]


title: How to plot the contour of f[x,y]==0 if always f[x,y]>=0

text: I've been stuck at this problem for weeks and I've asked several related questions here: ContourPlot shows only part of the contours [duplicate] How to plot the contour of the radius part of a complex function a lot of help has been kindly given but I think my problem is still not solved completely. In the first question, Jens pointed out that when a function only touch but not cross 0, `ContourPlot` would have difficulty finding the 0 contour. He then suggested to plot `f[x,y]==0.001` or some small number, instead of `f[x,y]==0.`. This seems works very well for some functions like `Abs[Cos[x]+Cos[y]]`. However it doesn't work on my function, and he also gave an explanation that the reason it doesn't work is because the maximals of my function change a lot in the plot region. Along this direction, I rescaled my function to make the maximals at the same order, but the `f[x,y]==0.001` trick still doesn't work. I've been experimenting with these for a while and found that the reason seem more than just the steep change in the maximals. And I also found some behavior of `ContourPlot` that are very confusing for me. Consider this function, the maximals of it changes quite of lot               f[x_, y_] := (Cos[x] + Cos[y]) Exp[I x y] Exp[x/2]               Plot3D[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi},       PlotPoints -> 40, PlotRange -> All]      ![enter image description here](http://i.stack.imgur.com/d8Hfm.png) Now if we plot the 0 contour, using different `PlotPoints` setting, we would get different results. In fact more plot points makes the contour worse. **Why?**               myPlot[sqs__] := Module[{},       Row@{ContourPlot[Sequence[sqs], ImageSize -> 400],          ListPlot[          Reap[ContourPlot[Sequence[sqs],              EvaluationMonitor :> Sow[{x, y}]]][[-1, 1]],           PlotStyle -> PointSize[Tiny], ImageSize -> 400, Axes -> False,           Frame -> True, AspectRatio -> 1]}       ]          myPlot[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi}]     myPlot[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi}, PlotPoints -> 50]     myPlot[Abs[f[x, y]] == 0.001, {x, 0, 4 Pi}, {y, 0, 4 Pi}, PlotPoints -> 100]      ![enter image description here](http://i.stack.imgur.com/uYbYW.png) If we increase the plot range, the `ContourPlot` just fails to give the result. **Why?** The MaxRecursion seems can help in this case, but the contour is still broken.               myPlot[Abs[f[x, y]] == 0.001, {x, 0, 16 Pi}, {y, 0, 16 Pi}]     myPlot[Abs[f[x, y]] == 0.001, {x, 0, 16 Pi}, {y, 0, 16 Pi}, MaxRecursion -> 3]      ![enter image description here](http://i.stack.imgur.com/XAVgT.png) Now compare with the second function               f[a_, a0_, k_, K0_] :=       a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k +           a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) +        2 k (I E^(           I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[a0 K0]) +           a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2])          With[{a0 = 10., a = 1.4},      Plot3D[Abs[f[a, a0, k, K0]] == 0.001, {K0, -2 π/a0,         2 π/a0}, {k, 0, 2}, ImageSize -> 400, PlotPoints -> 40]]      We can see that this function behaves well in terms of the changes in the maximals. All the maximals seems in the same order. And ContourPlot seems fails to produce the results, even plotting a very small region with `MaxRecursion->3`. Why? ![enter image description here](http://i.stack.imgur.com/fKTWe.png)               With[{a0 = 10., a = 1.4},      myPlot[Abs[f[a, a0, y, x]] == 0.001, {x, -π/a0, π/a0}, {y, 0,2}]]     With[{a0 = 10., a = 1.4},      myPlot[Abs[f[a, a0, y, x]] == 0.001, {x, -π/a0, π/a0}, {y, 0,2}, MaxRecursion -> 3]]      ![enter image description here](http://i.stack.imgur.com/rjt1J.png) **Summary of the Questions**   1. In the first function, why increase PlotPoints makes things worse, and why increases plot range resulted in almost no contour plotted(when plot range is 16pi)?   2. Why the second function, which seems behave better than the first, shows very poor contour results and how to improve it?


title: Problem with ContourPlot

text: I tried these               ContourPlot[ 1 + (1 + x)^3 - 3 (1 + x) y + y^3 == 0, {x, -3, 3}, {y, -3, 3}]     ContourPlot[(x + y + 2) == 0, {x, -3, 3}, {y, -3, 3}]     ContourPlot[(x + y + 2)^2 == 0, {x, -3, 3}, {y, -3, 3}]     ContourPlot[(x + y + 2)^3 == 0, {x, -3, 3}, {y, -3, 3}]      the third and fourth one produce graphs contrary to my mathematical common sense. Especially, the third one does not produce anything. Could someone explain why and how to fix it?


title: Plotting equipotentials gone wrong

text: I'm plotting the equipotential of the potential function: $$\phi(x,y,z):=\frac{1+x^{2}+y^{2}+z^{2}}{[\left(1+x^{2}+y^{2}+z^{2}\right)^{2}-4z^{2}]^{1/2}}$$ Mathematica code:               Φ[x_, y_, z_] := (1 + x^2 + y^2 + z^2)/((1 + x^2 + y^2 + z^2)^2 - 4 z^2)^(1/2)      So I write, for example               ContourPlot3D[Φ[x, y, z] == 1, {x, -100, 100}, {y, -100, 100}, {z, -100, 100}]      But, I get this image: ![enter image description here](http://i.stack.imgur.com/ElmCm.jpg) And, sometimes if I change the plot region, I get absolutely nothing. My **question** is, what did I wrong? I think I should obtain smooth 3D curves, but I don't.


title: Plotting a line segment from the degenerate case of an ellipse

text: An ellipse degenerates into a line segment when the defining constant distance from the two foci is the actual distance between the two foci. The ellipse closes into a line segment. The following code shows a family of concentric ellipses that should converge to a line segment when the constant is set to 5. However it doesn't plot the line segment even if I increase the MaxRecursion very high.               ContourPlot[Evaluate[Table[Sqrt[(x + 3)^2 + y^2] + Sqrt[x^2 + (y - 4)^2] == n, {n, 5,      6, .01}]], {x, -4, 1}, {y, -1, 5}]      Is there someway to get Mathematica to actually plot the line segment of the degenerate case of an ellipse?


title: ContourPlot shows only part of the contours

text: Consider these contour plots:               f[x_,y_]:=Cos[x] + Cos[y]          ContourPlot[f[x, y] == 0, {x, 0, 4 Pi}, {y, 0, 4 Pi}]      ![enter image description here](http://i.stack.imgur.com/bmuys.png)               ContourPlot[f[x, y]*f[x, y] == 0, {x, 0, 4 Pi}, {y, 0, 4 Pi}]      ![enter image description here](http://i.stack.imgur.com/fFQP9.png)               ContourPlot[Abs[f[x, y]] == 0, {x, 0, 4 Pi}, {y, 0, 4 Pi}]      ![enter image description here](http://i.stack.imgur.com/Nq2xe.png) **Question** : Why the second and third contour plot don't work? How to fix them to make them to reproduce the first contour plot? I tried to increase the `MaxRecursion` and it's supper slow and doesn't help much:               ContourPlot[Abs[f[x, y]] == 0, {x, 0, 4 Pi}, {y, 0, 4 Pi}, MaxRecursion -> 15]      ![enter image description here](http://i.stack.imgur.com/dCzko.png) * * * * * * This part is just my context of the problem, you can safely ignore it if you are not interested. Why do I want to plot the contour of `Abs[f[x,y]]==0` or `f[x,y]*Conjugate[f[x,y]]==0` if they are the same of `f[x,y]==0`? Because it would help a lot when `f[x,y]` is a complex function(x and y are real variables but `f[x,y]` maybe complex). For any function `f[x,y]`, it can be written in the form `r[x,y]*Exp[I θ[x,y]]` where `r[x,y]` and `θ[x,y]` are real function. So for `f[x,y]==0` what we really means is `r[x,y]==0`, and instead of contour plotting `f[x,y]==0` I should put `r[x,y]==0`. However, sometimes it would be difficult to get `r[x,y]` given the function `f[x,y]`, but fortunately we can just plot the contour of `f[x,y]*Conjugate[f[x,y]]==0` or `Abs[f[x,y]]==0` since which is mathematically equivalent to `r[x,y]==0`. For example, I have a function               a0 = 10.;a = 1.4;     f[k_,K0_]:=a^2 Sech[(a a0)/2]^2 (-2 I (1 + E^(2 I a0 k)) k + a (-1 + E^(2 I a0 k)) Tanh[(a a0)/2]) + 2 k (I E^(I a0 k) ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[a0 K0]) + a (-1 + E^(2 I a0 k)) k Tanh[(a a0)/2])      I know it can be separate into the form `r[k,K0]*Exp[I θ[k,K0]]`, but it's difficult to find what `r` and `θ` is. Since I want to plot the contour of `r[k,K0]==0`, instead I can plot the contour of `Abs[f[k,K0]]==0`, but I encountered the same problem as the simple example above.               ContourPlot[Abs[f[k, K0]] == 0, {K0, -2 π/a0, 2 π/a0}, {k, 0, 4}]      ![enter image description here](http://i.stack.imgur.com/u0KvT.png) However if we plot the real and imaginary part, we can get the correct contour(the diamond shape) but also introduced other artificial contours(horizontal lines) come from `Cos[θ[x,y]]==0` or `Sin[θ[x,y]]==0`.               Row@{ContourPlot[        Re[f[k, K0]] == 0, {K0, -2 π/a0, 2 π/a0}, {k, 0, 4}, PlotLabel -> "Re", PlotPoints -> 60],       ContourPlot[        Im[f[k, K0]] == 0, {K0, -2 π/a0, 2 π/a0}, {k, 0, 4}, PlotLabel -> "Im", PlotPoints -> 60]}      ![enter image description here](http://i.stack.imgur.com/lV3hz.png) Just for comparison, the above function can be separated into               2 I E^(I a0 k) (a^2 Sech[(a a0)/2]^2 (-2 k Cos[a0 k] + a Sin[a0 k] Tanh[(a a0)/2]) + k ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[a0 K0] + 2 a k Sin[a0 k] Tanh[(a a0)/2]))      so               r[k_,K0_]:=a^2 Sech[(a a0)/2]^2 (-2 k Cos[a0 k] + a Sin[a0 k] Tanh[(a a0)/2]) +       k ((a - k) (a + k) Cos[a0 k] + (a^2 + k^2) Cos[a0 K0] +          2 a k Sin[a0 k] Tanh[(a a0)/2])      and contour plot of r is               ContourPlot[r[k, K0] == 0, {K0, -2 \[Pi]/a0, 2 \[Pi]/a0}, {k, 0, 4}, PlotPoints -> 60]      ![enter image description here](http://i.stack.imgur.com/IcZfR.png)


title: Bar Chart Labels Mystery

text: I have some data which I want to plot in a simple bar chart:               example = {{"Knowledge", 2}, {"Comprehension", 4}, {"Application", 0}, {"Analysis", 1}, {"Synthesis", 0}, {"Evaluation", 2}}      The first element of each pair should be the bar label, and the second should be the height of the bar. I can get the desired result using:               BarChart[Last[#] & /@ example, ChartLabels -> {"Knowledge","Comprehension","Application","Analysis","Synthesis","Evaluation"}]      ![enter image description here](http://i.stack.imgur.com/y6FKi.png) But my real question is this: Though               First[#] &/@ example      Returns               {"Knowledge", "Comprehension", "Application", "Analysis", "Synthesis", "Evaluation"}      When I try:               BarChart[Last[#]&/@example, ChartLabels->First[#]&/@example]      I get a graph without labels: ![enter image description here](http://i.stack.imgur.com/y6ydg.png) This seems like strange behaviour, can anyone shed some light, or suggest a fix?


title: Numerical contour integrations in the complex plane - contour deformation gives different answer for analytic kernel

text: I am trying to do a contour integration in _Mathematica_ numerically. In particular, I'm checking the identity: $$ H_m^{(1)}(z) =\frac{i^{-m}}{\pi}\int_{-\pi/2 + i \infty}^{\pi/2 - i \infty} \exp[i m \beta + i z \cos\beta]\mathrm d\beta $$ In Stratton's _Electromagnetic Theory_ (I am a physicist), he claims (p. 368) that you can freely deform the contour since the kernel is an analytic function, assuming the contour stays within a certain region. Here's a diagram (from here): ![admissible region for Hankel function contour](http://i.stack.imgur.com/l0ESD.png) This is all fine, and I have no problem with this; however, when I try to evaluate the contour integration in _Mathematica_ with code that looks like:               Hmz[m_, z_] := Exp[-I m Pi/2]/Pi NIntegrate[Exp[I m B + I z Cos[B]],                      {B, -Pi/2 + I Infinity, -Pi/2 + I Pi/2, Pi/2 - I Pi/2, Pi/2 - I Infinity}]      or alternatively (along a different, but apparently valid, contour):               Hmzother[m_, z_] :=  Exp[-I m Pi/2]/Pi NIntegrate[Exp[I m B + I z Cos[B]],                                       {B,  -Pi/2 + I Infinity, -Pi/2, Pi/2, Pi/2 - I Infinity}]      I get different answers - at least for the imaginary part (particularly for small $z$). For example, evaluating for $m = 3$ and $z = 1$ gives (not actual lines of code):               Hmz[3, 1] == 0.0195634 - 2.72402 I     Hmzother[3, 1] == 0.0195634 - 0.248805 I      For reference, the true value should be `HankelH1[3,1] == 0.0195634 - 5.82152 I`. The only difference between the definition of `Hmz` and `Hmzother` here are the contours, yet they should give the same answer. Ideally, I'd like to use the contour `{B, I Infinity, 0, Pi, -I Infinity}`, but this gives an even worse disagreement. Does anybody know how to get _Mathematica_ to do the contour integral correctly in the complex plane? Thanks all.


title: Preventing notebook cells from being broken across two pages when printing

text: Sometimes I create a number of plots together, e.g. by `GraphicsGrid[]`. If it happens that I have to print the notebook, very often the cell is broken, half on one page and the remaining part on the next -- while I'd like _Mathematica_ to print it all in the next page if there is not enough space available. This also happens if I save the notebook as a PDF.


title: Confidence Regions for NonLinearFit with “complicated” integral model

text: I have a set of data $y_i (z_i)$ with errors $\Delta y_i$, that can be found here Here Thus first, I call the data               dataSN = Import["Cat3_0mod.txt", "Table"];     sigdata = dataSN[[All, 4]];     lendata = Length[zdata];          fitSN = Table[{dataSN[[i, 2]], dataSN[[i, 3]]}, {i, 1, lendata}];      I want to fit these data to a complicated model $$ Y(z) = 5\log \Bigl((1+z)\int_0^z \frac{dz'}{a+bz'+cz'^2}\Bigr)+25; $$ With help of a previous post NonLinearFit with "complicated" integral model I used 'NonLinerFit'               ModelSN1[z_] =        Simplify[5 Log[10, (1 + z) Integrate[1/(a + b x + c x^2), {x, 0, z}]] + 25,          Assumptions -> {z > 0, -b^2 + 4 a c > 0}];          fit =        NonlinearModelFit[fitSN, {ModelSN1[z], -b^2 + 4 a c > 0},         {{a, 0.000254}, {b, .000016}, {c, .000016}}, z, Weights -> 1/(sigdata)^2];      Even the fit looks fine ![enter image description here](http://i.stack.imgur.com/u0GAe.jpg) The reality is that the fit is complex ![enter image description here](http://i.stack.imgur.com/xxzF9.jpg) I Have some problems making the fit   1. Even I put conditions over the parameters $a,b,c$ the fitted model is complex. I tried to put as condition that the whole term inside $\log$ should be positive, but it does not work.   2. I want to make confidence regions on this fit using               band1SN[z_] = fit["MeanPredictionBands", ConfidenceLevel -> 0.68];       band2SN[z_] = fit["MeanPredictionBands", ConfidenceLevel -> 0.95];      but _Mathematica_ does not make it, an error related with the fact of assuming a constrained model appears. Of course, if the fit is complex, the confidence regions could not be calculated. I appreciate your help.


title: Displaying Values of Fuction at end of Period by Grid

text: My intention is to place the values of two function into a small box within my chart. To be more precise I have two vectors. A meanvector which shows the average value of all simulated stochastic processes and an other vector, namely thresholdedmeanvector which represents the average value only of those processes which remain beyond the in the function defined threshold. So I just want a box somewhere in the chart with the titles “mean of all stochastic processes” and below the displayed value of this vector at the end of the period (which is defined as “time” in the code) and the title “mean of all stochastic processes beyond the threshold” next to it, which also displays the value of the average of the stochastic processes which remained beyond the threshold at the end of the period (again defined as “time”) In the first step, I tried just to display the mean value of all stochastic processes somewhere in the chart. Therefore I uses the Grid Function, but I am struggling with it. I think I have placed it somewhere in a wrong position in the function, but I am not sure. Does anyone of you have a suggestion where I should place the Grid Function or an easier alternative? Thanks, and below you can find the underlying code.               Manipulate[SeedRandom[seed]; meanvector := Mean[assets];       assets = Table[RandomFunction[     GeometricBrownianMotionProcess[μ, σ, S0], {0, time, 0.1}]["Path"], {P}];                processes = Transpose[assets[[#]]][[2]] & /@ Range@Length[assets];       processesposition = Flatten[Position[     Min[processes[[#]]] & /@ Range@Length[assets], _?(# > threshold &)]];                      thresholdedassets = assets[[#]] & /@ processesposition;       thresholdedmeanvector = Mean[thresholdedassets];                (*vector of meanvector*)      vmv = TableView[{meanvector, NumberForm[meanvector, {4, 5}]}, {0, time}];            vmv = Prepend[        vmv, {"Mean Value of all Assets and the end of the period"}];           G1 := ListLogPlot[assets, GridLines -> {{}, {threshold}},         GridLinesStyle -> Directive[Green, Thick], Joined -> True,         AxesLabel -> {"Time", "St"},         PlotLabel -> Style["Forecasted Stock Price\n (Brownian Motion)", Bold],         PlotRange -> All, PlotStyle -> Directive[{Thin, Lighter@Gray}]];           G2 := ListLogPlot[{meanvector, thresholdedmeanvector},         Joined -> True,         PlotStyle -> {Directive[{Thick, Darker@Red}],           Directive[{Thick, Darker@Blue}]}];           (*introduction of Grid Box for meanvector*)           Grid[{Text@ Grid[vmv, Alignment -> Left, Frame -> All,ItemStyle -> {Automatic, {Bold}}],                       Show[G1, G2]}],           {{S0, 100, "Initial Stock Value"}, 1, 500, 0.5, Appearance -> "Labeled"},       {{μ, 0.08, "Drift μ"}, 0.01, 0.2, 0.01, Appearance -> "Labeled"},       {{σ, 0.2, "Standard Deviation σ"}, 0.01, 1, 0.05,        Appearance -> "Labeled"},      {{P, 6, "Paths"}, 1, 20, 1,Appearance -> "Labeled"},       {{time, 10, "Time t"}, 1, 20, 1, Appearance -> "Labeled"},        {{threshold, 75, "Debt Level"}, 1, 500, Appearance -> "Labeled"},        Delimiter, {{seed, 0, ""},Button["randomize", seed = RandomInteger[2^64 - 1]] &},             Button["Set Initial Values", {S0 = 100, μ = 0.08, σ =          0.20, P = 6, time = 10, threshold = 75}, ImageSize -> 150],       ControlPlacement -> Left]


title: Interpreting the interpolating function and saving data to plot with external program

text: So far I have been solving non-linear pdes with `NDSolve` and then plotting the result with the in-built `Plot3D` and `ContourPlots`. An example file can be found **here**. This is a two part question with part 1 leading in to part 2. ## Part 1 How do I interpret the interpolating function (how do I read it)? For instance, I can "expand" the `InterpolatingFunction` with > `InterpolatingFunctionGrid[hSol]` This obviously has all the points that `Plot3D` plots with. How do I read this, say, row-wise or column-wise to make more sense of whats going on? ## Part 2 How can I save this data to a `.csv` format or just some external file to use, say, gnuplot or some other plotting routine? I have trid `DumpSave` etc . but the resulting huge files are generally not cross platform compatible.


title: Multiple NIntegral with singularity

text: I need to compute the following integral `FF[0]`. When evaluating `FF[0]` _Mathematica_ simply takes too much time. The problem may arise due to singularity in the limit 1. I tried specifing different options of `NIntegrate`, but no luck. Any suggestions why the computation of `FF[0]` stalls ? I am working on the code below:               FF[z_] := NIntegrate[f1[v] NIntegrate[f4[u], {u, 1, v}], {v, 1, z}];      It involves the functions below:               f1[z_] := z^2/((1 + z + z^2) Sqrt[1 - z^3]);          f2[z_] := (-2 - 3 z^2 + 5 z^3 + 9/2 z^2 Sqrt[1 - z^3]                 NIntegrate[f1[u], {u, 1, z}])/(-1 + z^3);          f3[z_] := (z^2 (8 - 12 z - 2 z^3 + 6 z^4 + 3 Sqrt[1 - z^3] (-4 + z^3)                 NIntegrate[f1[u], {u, 1, z}]))/(-1 + z^3)^2;          f4[z_] := 1/((1 + z + z^2) Sqrt[1 - z^3]) (16 (-1 + z) (1 + z + z^2) Sqrt[1 - z^3] +                 6 (-4 - 4 z - 4 z^2 + z^3 + z^4 + z^5) NIntegrate[f1[u], {u, 1, z}] +                 Sqrt[1 - z^3] (-4 (-4 + 2 z + 2 z^2 + 3 z^3) +                4 (1 + z + z^2) NIntegrate[f2[v], {v, 1, z}] -                3 (1 + z + z^2) NIntegrate[f3[s], {s, 1, z}]));


title: Manipulate not showing anything

text: This works:               Manipulate[Plot[2*x + 3*y, {x, 3, 4}], {y, 3, 4}]      But this doesn't               f=2*x + 3*y     Manipulate[Plot[f, {x, 3, 4}], {y, 3, 4}]      My "f" is a painful long expression so I'd rather be able to get the latter method to work than work with the former.


title: Simulated Annealing Convergence

text: I am using _Simulated Annealing_ method for a simulation based optimization of a process that has 3 variables, using `NMinimize`. I print the input/output during every iteration using the `"Evaluation Monitor"`. I have noticed that after about `3000` iterations, _Mathematica_ gives me a convergence result. But checking the results from `"Evaluation Monitor"` shows that the ultimate convergence result is not the global minimal, but the process has visited a better minimal during one of the iterations. Why is this happening? Here is a sample code:               demand[n_,k_]:=Min[k Vf,n capacity];     supply[n_,k_]:=Min[(n Kj-k) w, n capacity];     flo[n_,Ku_,Kd_]:=Min[demand[n,Ku],supply[n,Kd]];     dx=Vf*dt; n=Round[Flen/dx]; m=Round[SimTime/dt];       p=Round[Rlen/dx];   θ=Vf/w; capacity=w*Vf*Kj/(Vf+w); α[a1_]:=1800.; β[a2_]:=0.1; L=3.;      Flen=8.; Rlen=3.; SimTime=30./60.; Kj=150.; w=20.; Vf=100.; dt=12./3600.; d=1.;      RMLocation=Round[(2/3) p]; j=0;          f[a1_,a2_,a3_]:=Module[{b1=a1,b2=a2,b3=a3,TT=0,NtwrkTT=0,j=0},        RM[x_,t_]:=Piecewise[{{100 b1,x<=3},{100 b2,3<x<=6},{100 b3,True}}];         NtwrkTT=0; Clear[k0,kr,k,γ];        k0=ConstantArray[0,n];         kr=Table[Table[0,{k,1,p}],{i,1,n}];         γ=ConstantArray[1,n];        For[i=2, i<n, i++, kr[[i,1]]=α[i dx] d/Vf];         TT=Plus@@(Plus@@kr); NtwrkTT=TT; k=k0;        While[TT>0, TT=0;          For[i=2, i<n, i++, FQin=If[i==2,Min[demand[L,k0[[i-1]]],supply[L,k0[[i]]]],FQout];       dem=demand[L,k0[[i]]]; dem=If[dem==0,0.001,dem];       γ[[i]]=Min[1,supply[L,k0[[i+1]]]/dem];       ϕ=γ[[i]] demand[1,kr[[i,p]]]/d;       Qr=(ϕ-β[i dx] FQin) dx;       FQout=Min[demand[L,k0[[i]]],supply[L,k0[[i+1]]]];       k[[i]]=k0[[i]]+(FQin-FQout+Qr)/Vf;       kr0=kr[[i]];       For[ir=2,ir<=p,ir++,        MR=If[ir==RMLocation+1,RM[i dx,j dt],capacity];        RQin=Min[MR,If[ir==2,flo[1,kr0[[ir-1]],kr0[[ir]]],RQout]];        MR=If[ir==RMLocation,RM[i dx,j dt],capacity];        RQout=Min[MR,If[ir<p,flo[1,kr0[[ir]],kr0[[ir+1]]],ϕ d]];        kr[[i,ir]]=kr0[[ir]]+(RQin-RQout)/Vf];       kr[[i,1]]=If[j<=m,α[i dx] d/Vf,0]];      TT=Plus@@(Plus@@kr);      TT+=Plus@@k;      k0=k;NtwrkTT+=TT;j++];     NtwrkTT dt]     NMinimize[{f[a,b,c],3<=a<=12&&3<=b<=12&&3<=c<=12&&Element[a|b|c,Integers]},{a,b,c},Method->{"SimulatedAnnealing","SearchPoints"->5^5},EvaluationMonitor:>Print["a = ",a," , b = ",b," , c = ",c," , f[a,b,c] = ",f[a,b,c]]]      Hope this helps. Ps. Any suggestions to improve the performance of this code will be greatly appreciated. Thank You. **edit:** As per @@Kuba's suggestion, I split the original question in to two. The second part is at Simulated Annealing Parameters and Results


title: GraphPlot with black edges no rendering function

text: Is there any way to make the edges of a GraphPlot black without tediously writing out a long "EdgeRenderingFunction"? I guess you could ask the same kind of question about vertex labeling. For example, if I write `EdgeRenderingFunction -> ({Black, Arrow[#1, 0.1]} &)` I don't even get the standard directed edges, I get "arrows". How do I just turn the standard directed edges black (instead of the default red)?


title: Is it possible to lower the bitwidth of (floating point) numbers?

text: The answer to this question will most probably be a definite no. Nevertheless, I would like to know whether it is possible to represent reals in the memory on less than the standard 32 or 64 bits (depending on machine). This limits my exploration to $\sim 2^{14} \times2^{14}$ matrices, though with half or quarter precision I could venture a bit further. The matrix itself cannot be broken down to smaller pieces for further compuation, and it cannot be reworked to a simpler one, but high precision is of no importance. The question is not about refactoring matrices but about number representation and memory handling by _Mathematica_. **Warning, the following computation may eat through your memory limits.**               n = 14;     all = (2^n)*(2^n);     m = RandomReal[{0, 1}, {2^n, 2^n}];     {all, (ByteCount@m/1024./1024) "MB", ((ByteCount@m/all*8.) // N) "bit"}     ClearAll[m];      > >     {268435456, 1024. "MB", 64. "bit"} >   I am using 64-bit Win7, with 2 GB RAM, and (expectedly) my system almost chokes on this code. And then of course I hadn't calculated _anything_ with the matrix yet...


title: Why doesn't ReplaceAll work on Image?

text: `ReplaceAll` seems not to reach inside `Image`, and I don't see why it shouldn't. Here is an example. Create a black and white image with four pixels:               image = Image[{{1, 0}, {0, 1}}]      The full form of `image` is:               Image[List[List[1.`, 0.`], List[0.`, 1.`]], "Real", Rule[ColorSpace, Automatic], Rule[Interleaving, None]]      The following line switches the pixels as desired:               Image[ImageData[image] /. n_?NumberQ -> 1 - n]      My question is, in the following usage, why doesn't `ReplaceAll` reach inside `Image` and accomplish the same?               image /. n_?NumberQ -> 1 - n


title: BarChart without axes neglects the last zero of a list

text: Creating a BarChart with {1, 2, 3, 0} as               BarChart[{1, 2, 3, 0}, Axes -> False, Frame -> True]      gives me only the first three bars, a chart identical to just               BarChart[{1, 2, 3}, Axes -> False, Frame -> True]        1. Is this an intended behaviour or a bug?   2. (Originally I asked for a workaround. I think I can just use `Axes` instead of `Frame`. Let me ask this instead for the sake of knowledge.) What to do if I need `Frame` but not `Axes`? I'm using Mathematica 9.


title: NDSolve for a large system of simple ODEs

text: I am solving a system of many (more than 100) ODEs. It is the kind of standard rate equation encountered in semiconductor physics. Here is the system:               eqns := Table[       P[k]'[t] ==          Sum[-P[k][t] P[q][t] a[k,q,t] - P[k][t] (1 - P[q][t]) b[k,q,t] +              (1 - P[k][t]) P[q][t] c[k,q,t] + (1 - P[k][t]) (1 - P[q][t]) d[k,q,t],              {q, Drop[Table[i, {i, 1, M}], {k}]}]          - e[k,t] (P[k][t])^2 + f[k,t] (1 - P[k][t])^2,       {k, 1, M}];      where `P[k][t]` represent probabilities (bounded between 0 and 1). The form of the solution for every `P[k][t]` is very close to an analytical one, and it is very smooth (basically has just one maximum, quite like a gaussian). I am interested in gettign `Sum[P[k][t], {k, 1, M}]`. I solve the system with               functionlist = Table[P[k], {k, 1, M}];     initialconds = Table[P[k][t0] == initialvalue[[k]], {k, 1, M}];     eqnstosolve = Join[eqns, initialconds];          NDSolve[eqnstosolve, functionlist, {t, t0, tfin},        Method -> {"EquationSimplification" -> "Solve"}];      Question: is there a better way? I tried several method and this seems the faster. However, for my large system of equation evaluation takes hours. I would like, for example, to use a very fast and sloppily method at the beginning and then to refine the peak. Or I would just like to solve faster than this method by adjusting the steps. So far everything I've tried just make things worse. I not even sure what `EquationSimplification->Solve` really does. Any help is appreciated. EDIT: since I have suggested to paste a running code, here is my actual code. The problem is what I stated previously in a more polished way:               Ns = 40;     \[Alpha] = 0.01;     T=0.1;     v=0.0001;          kEven[n_, N_] := 2 \[Pi] n/N  - (N - 1) \[Pi]/N;     kPos[Ns_] := Table[kEven[n, Ns], {n, Ns/2, Ns - 1}];     \[Epsilon][g_, k_] := Sqrt[1 + g^2 - 2 g Cos[k]];          Jf = 1.0;     \[Theta][k_, g_] := ArcTan[g - Cos[k], Sin[k]];     J[\[Omega]_, \[Alpha]_, s_, \[Omega]c_] :=        2 \[Alpha] \[Omega]^         s Exp[-\[Omega]/\[Omega]c] HeavisideTheta[\[Omega]];     gOhm[E_, T_] :=       If[E == 0,        2 N[\[Pi]] \[Alpha] T, \[Pi] ((J[(E + 10^(-20)), \[Alpha], 1,               10^6] -              J[-(E + 10^(-20)), 0.01, 1, 10^6]) Exp[(E + 10^(-20))/             T] )/(Exp[(E + 10^(-20))/T] - 1)]               h[t_] := HeavisideTheta[1 - v t] (1 - v t);     tin = -2/v;     relaxfact = 1;     tfin = relaxfact/v;              eqnsPOS := Table[P[k]'[t] ==((1/(2 Ns)) Sum[                 -P[k][t] P[q][t] (W[+1, +1, k, q] + W[+1, +1, k, -q])                  - P[k][t] (1 - P[q][t]) (W[+1, -1, k, q] + W[+1, -1, k, -q])                  + (1 - P[k][t]) P[q][t] (W[-1, +1, k, q] + W[-1, +1, k, -q])                  + (1 - P[k][t]) (1 - P[q][t]) (W[-1, -1, k, q] +                      W[-1, -1, k, -q])                 , {q, Drop[Table[i, {i, 1, Ns}], {k}]}]               - (1/(2 Ns)) W[+1, +1, k, k] (P[k][t])^2               + (1/(2 Ns)) W[-1, -1, k, k] (1 - P[k][t])^2)           , {k, 1, Ns}];                             kPos[2 Ns];         Clear[W];         W[s1_, s2_, k_, q_] :=           W[s1, s2, k,             q] = (1 -               s1 s2 Cos[                Sign[k] \[Theta][kPos[2 Ns][[Abs[k]]], h[t]] +                  Sign[q] \[Theta][kPos[2 Ns][[Abs[q]]],                    h[t]]]) (gOhm[(s1 \[Epsilon][h[t], kPos[2 Ns][[Abs[k]]]] +                 s2 \[Epsilon][h[t], kPos[2 Ns][[Abs[q]]]]), T]);              functionlist = Table[P[k], {k, 1, Ns}];              thermalPsFD[j_, g_] := 1/(Exp[\[Epsilon][g, kPos[2 Ns][[j]]]/T] + 1);                   thermalPsInit =           Table[1/(Exp[\[Epsilon][h[t], kPos[2 Ns][[k]]]/T] + 1), {k, 1,          Ns}] /. t -> tin;          initialconds = Table[P[k][tin] == thermalPsInit[[k]], {k, 1, Ns}];          eqnstosolve = Join[eqnsPOS, initialconds];          time = AbsoluteTiming[       spos = NDSolve[eqnstosolve, functionlist, {t, tin, tfin},           Method -> {"EquationSimplification" -> "Solve"}];       ];     Print["time NDSolve" <> ToString[time[[1]]/60.] <> ToString[" min"]];


title: How to format output when using Manipulate

text: 1. As in the title, I want to output the image of a 3D body with its volume and surface area shown nicely beside it. (Maybe in a table?) But how can I do it?   2. Is it possible to calculate the surface area for each specific facets?  For example, consider the simple case:               Manipulate[       {RegionPlot3D[x + y + z < A, {x, 0, 1 - A}, {y, 0, 1 - A}, {z, 0, 1 - A},           PlotRange -> {0, .5}],        NIntegrate[Boole[x + y + z < A], {x, 0, 1 - A}, {y, 0, 1 - A}, {z, 0, 1 - A}]},        {A, 0, .99}]      What I can get so far is showing both the plot and data but not formatting it well. ![Current progress](http://i.stack.imgur.com/Di2BY.png)


title: Define new number set

text: May seem strange at first, but can you make mathematica change the value of 0 to 10 and 11 to 1 for all occurrences. Let me explain: I'm looping through a matrix and for each element, I want something to happen depending on the values in each element above, below, to the left and to the right. When you're at the first element however there is nothing to the left, or above, so instead I would like to use the values on either side of the matrix as if you'd curled the matrix back round onto itself, as shown below. ![](http://i.stack.imgur.com/9i2Uo.png) X denotes what element you're looking at and then T,B,L & R, correspond to the cells at the Top, bottom, to the left & right respectively. Here's the code I wrote without thinking too much and it throws up an error because it asks for elements that do no exist. (e.g. "i-1" when i=1)               Do[       Do[     {n, m} = RandomChoice[{image[[i + 1, j]], image[[i, j + 1]],      image[[i - 1, j]], image[[i,j - 1]]} -> {{i + 1, j}, {i, j + 1}, {i - 1, j}, {i,j - 1}}];     newimage[[n, m]] = image[[i, j]];     , {j, columns}];       , {i, rows}];      So, without whacking loads of If statements in there, is there an elegant way to accomplish this? Hope I explained this half decently.


title: Stopping- threshold and stopping- time of stochastic processes.

text: Rod Lm helped me to generate the following code, which I want to extend a little bit. So my intention is to use the threshold (the green line) as boundary in the following way: as soon as the process hits the threshold (watermark) it stops and is also excluded from the calculation of the mean function. So all processes which equal or fall beneath the threshold should stop and all processes that remain beyond the threshold should continue and the mean should only be calculated from the processes that remain beyond the threshold. Is this or something similar even possible? Thank you. the code:               Manipulate[SeedRandom[seed];     meanvector := Mean[assets];     assets = Table[     RandomFunction[GeometricBrownianMotionProcess[μ, σ, S0], {0, time, 0.1}]["Path"], {P}];     G1 := ListLogPlot[     assets, GridLines -> {{}, {threshold}}, GridLinesStyle -> Directive[Green, Thick], Joined -> True, AxesLabel -> {"Time", "St"}, PlotLabel -> Style["Forecasted Stock Price\n (Brownian Motion)", Bold], PlotRange -> All, PlotStyle -> Directive[{Thin, Lighter@Gray}]];          G2 := ListLogPlot[     Mean[assets], Joined -> True, PlotStyle -> Directive[{Thick, Darker@Red}]];         Show[G1, G2],         {{S0, 25, "Initial Stock Value"}, 1, 500, 0.5, Appearance -> "Labeled"},         {{μ, 0.08, "Drift μ"}, 0.01, 0.2, 0.01, Appearance -> "Labeled"},           {{σ, 0.2, "Standard Deviation σ"}, 0.01, 1, 0.05, Appearance -> "Labeled"},             {{P, 6, "Paths"}, 1, 20, 1, Appearance -> "Labeled"},             {{time, 10, "Time t"}, 1, 20, 1, Appearance -> "Labeled"},             {{threshold, 25, "Watermark"}, 1, 500, Appearance -> "Labeled"},             {{seed, 1, "New Random Case"}, 1, 100, 1},     Button["Set Initial Values", {S0 = 25, μ = 0.08, σ = 0.20, P = 6, time = 10, threshold = 25}, ImageSize -> 150],                 ControlPlacement -> Left]


title: Running out of memory using FindMinimum with ParametricFunctions

text: I am running a pretty large optimization problem which actually consists of an outer an inner optmization problem. I use ParametricNDSolve for solving a ODE system. When the optimization is done, Mathematica uses almost all of my RAM memory. If I try Remove["Global`*"] the memory used is still the same. Is there an easy way to check where the RAM in used?


title: Plotting a Complex Valued Function

text: Suppose that a complex valued function `f[z]=f[x+i y]=u[x,y]+i v[x,y]`. I want to plot3d its real part by coloring it with respect to its imaginary part. For instance, if `m` and `M` are the smallest and greatest values of the imaginary part, I want to color the graphic from white to black depending on the value of the imaginary part (`m` is white, `M` is black). I would be very glad if you can help me.


title: Compile and "True should be a machine-size real number..." Error

text: I am trying to compile a little summation function. As far as I understand, because the built-in `Sum` function can also return `NaN` objects I have to make sure only numbers get out before returning.               X = Compile[{}, (        result = Sum[i, {i, 10}];        If[NumberQ[result], Return[result], Abort[]];        )];      The problem is that I keep getting the following two warnings... > CompiledFunction::cfse: Compiled expression True should be a machine-size > real number. CompiledFunction::cfex: Could not complete external evaluation > at instruction 2; proceeding with uncompiled evaluation. I guess the glitch is in the `NumberQ` function and its return type. Any ideas on how to handle this issue correctly and effectively? **EDIT:** When I omit the check for numeric values while inserting a Return[] statement Mathematica falls back to the uncompiled version (warning: Compile::cret).


title: Is there a way to recreate the typical Red/Blue-Postereffect using Mathematica?

text: I think nearly all of us have at least once seen an image like this: ![enter image description here](http://i.stack.imgur.com/Zxomk.jpg) So I wondered whether it is possible to recreate this effect using Mathematica because it seems like it wouldn't need that many processing steps. Any ideas? http://hplussummit.com/images/wolfram.jpg this can be used as a source image.


title: Can a Table iterator "leak" into a Module?

text: I want to repeatedly call a module by using a Table statement. When I use an `i` iterator in the Table command calling the module, I receive the error message "Table::itraw: Raw object 1 cannot be used as an iterator. >>" The module also contains a separate Table statement with the `i` iterator. I noticed that when I change the iterator in the Table statement calling the module to something else, say `v` I no longer get the error message above. This makes me think the `i` iterator is leaking into the module, so .... _Question 1_ : Am I correct that the value of the `i` iterator from the Table statement calling the Module is leaking into the module? _Question 2_ : Is this expected behavior? _Question 3_ : Why does this happen? _Question 4_ : Is there another way to solve this problem, other than to avoid the use of `i` as an iterator in both the Table calling the module and the module itself? Thank you for any advice! Todd * * * Here is a snippet of the offending code:               errorcorrection[signaltype_,chipid_,probesetpostocheck_]:=Module[{clusters,positions,loopindex,groupminuscurrent,similarprobesets},          loopindex=1; (* initialize while loop index to 1; this loop continues for the length of the number of probesets that need checked for unusual signal - which is the length of the probesetpostocheck *)          While[loopindex<=2(*Length[probesetpostocheck[chipid]]*),          clusters=Flatten[Drop[Sort[FindClusters[signaltype[chipid][[probesetpostocheck[chipid][[loopindex]]]],3,DistanceFunction->EuclideanDistance,Method->"Agglomerate"]],-1]]; (* locate probe signals within the current probeset being evaluated that are unreliable and need replaced *)          Print[clusters];          With[{i=i},positions=Flatten[Table[Position[signaltype[chipid][[probesetpostocheck[chipid][[loopindex]]]],clusters[[i]]],{i,1,Length[clusters]}]]]          Print[positions];          If[MemberQ[controlchips,chipid],          (groupminuscurrent=DeleteCases[controlchips,x_/;x==chipid];     similarprobesets=Table[signaltype[groupminuscurrent[[i]]][[probesetpostocheck[chipid][[loopindex]]]],{i,1,Length[groupminuscurrent]}]),          (groupminuscurrent=DeleteCases[experimentchips,x_/;x==chipid];     similarprobesets=Table[signaltype[groupminuscurrent[[i]]][[probesetpostocheck[chipid][[loopindex]]]],{i,1,Length[groupminuscurrent]}])]          Print[similarprobesets];          loopindex++]]      * * *               run2 = Table[       errorcorrection[pmsignal, experimentchips[[i]], pmpostocheck], {i,         1, 1}]               {11719.,1567.,2149.,2001.,1781.,1693.}          Table::itraw: Raw object 1 cannot be used as an iterator. >>          Table[Position[pmsignal[GSM356796][[pmpostocheck[GSM356796][[loopindex$112290]]]],clusters$112290[[1]]],{1,1,Length[clusters$112290]}]      * * *               run3 = Table[       errorcorrection[pmsignal, experimentchips[[v]], pmpostocheck], {v,         1, 1}]          {11719.,1567.,2149.,2001.,1781.,1693.}          {7,6,8,9,17,20}      * * * After considering everyone's responses below, I have been able to get my code to work better by replacing the `With[{i=i},positions=Flatten ....` statement inside the errorcorrection module with `Module[{i=i},positions=Flatten ....` As noted `Block` also works fine. I am still puzzled, however, why `i` takes on a value of 1 inside the `errorcorrection` module, when the module is called by this statement:               run2 = Table[errorcorrection[pmsignal, experimentchips[[i]], pmpostocheck], {i, 1, 1}]      An example of output showing the iterator `i` is here:               run2 = Table[errorcorrection14[pmsignal, experimentchips[[i]], pmpostocheck], {i, 1, 1}]               {11719.,1567.,2149.,2001.,1781.,1693.}      `i` equals 1               {7,6,8,9,17,20}      _Question 5_ : (continuing from numbering above): Doesn't the `Table` command in `run2=Table[errorcorrection14[pmsignal, experimentchips[[i]], pmpostocheck], {i, 1, 1}]` localize the iterator? _Question 6_ : If so, why does that value show up in the `errorcorrection` module? _Question 7_ : If `i` did NOT take on the value of 1 (or any value) in the `errorcorrection` module, would my original `With[{i=i},positions=Flatten` statement have worked ok? It seems to me (in my naivety) that the problem with using With in this case, is that `i` is taking on a value in the module from the `Table` statement calling the module .... but I thought Table should localize `i` and prevent its value from going into the module?? Thank you again, as I struggle to understand at a "deeper" level.


title: NDSolve Convergence test failure and significant effect of DifferenceOrder on eventual results

text: I am solving a non linear partial differential equation with what I call free boundary conditions (solid mechanicists would know this as _simply supported_ ). I realized that this boundary condition does induce a certain stiffness in my equation which is not so strong for, say, periodic boundary condition. How do I know this? The simulation runs longer for free boundaries with a stiff method such as `LSODA`. This tells me that the time step is being made smaller and smaller to accommodate this rapidly changing time variation. I also have various _structures_ forming as obvious from the plot of my result. So I use the `TensorProductGrid` spatial discretization method as I gather it would be useful to solve this equation with a preferential spatial grid treatment (denser grid in regions of necessity). As I increase the `DifferenceOrder` for my method/spatial discretization, I notice that I run into convergence issues `NDSolve::ndcf` at a later stage. Viz., for a smaller `DifferenceOrder` I have a convergence issue at a later time. This initially struck me as being a little strange as wouldn't a superior difference order help alleviate convergence? But then obviously, a greater difference order would catch the issue with convergence earlier! ## Here is my code:               m=0.05;     \[Epsilon]=10^-6;     Bo=1/300;     \[Delta]=10^-7;     Bi=1;     K1=1;     wn=2;          (*R=Ra/(Pr S1);*)     R=0*1/10;     {xMin,xMax}=wn{0,79.5788};     k=0.0677*92.389/(xMax/wn);     TMax=10000*100;     n=1;               uSolpbc[t_,x_]=u[t,x]/.NDSolve[{\!\(     \*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x]\)\)+\!\(     \*SubscriptBox[\(\[PartialD]\), \(x\)]\((     \*SuperscriptBox[\(u[t, x]\), \(3\)]\      \*SubscriptBox[\(\[PartialD]\), \(x, x, x\)]u[t, x])\)\)-Bo \!\(     \*SubscriptBox[\(\[PartialD]\), \(x\)]\((     \*SuperscriptBox[\(u[t, x]\), \(3\)]\      \*SubscriptBox[\(\[PartialD]\), \(x\)]u[t, x])\)\)+m  \!\(     \*SubscriptBox[\(\[PartialD]\), \(x\)]\((     \*SuperscriptBox[\((     \*FractionBox[\(u[t, x]\), \(K1 + Bi\ u[t, x]\)])\), \(2\)]\      \*SubscriptBox[\(\[PartialD]\), \(x\)]u[t, x])\)\) +\[Epsilon]/(Bi u[t,x] + K1) +\[Delta] u[t,x]^3/(Bi u[t,x] + K1)^3 D[D[u[t,x],x],x] + R \!\(     \*SubscriptBox[\(\[PartialD]\), \(x\)]\((     \*SuperscriptBox[\(u[t, x]\), \(3\)]\      \*SubscriptBox[\(\[PartialD]\), \(x\)]u[t, x])\)\)== 0,          u[0,x]==1-0.05 (Cos[k*x]),          (*PERIODIC BOUNDARY CONDITIONS-START*)     (*u[t,xMin]== u[t,xMax],     Derivative[0,1]u[t,xMin]==Derivative[0,1]u[t,xMax],     Derivative[0,2]u[t,xMin]==Derivative[0,2]u[t,xMax],     Derivative[0,3]u[t,xMin]==Derivative[0,3]u[t,xMax]*)     (*PERIODIC BOUNDARY CONDITIONS-END. Uncomment to run periodic case*)          (*FREE BOUNDARY CONDITIONS-START*)          Derivative[0,1][u][t,xMin]==0,Derivative[0,1][u][t,xMax]==0,Derivative[0,3][u][t,xMin]==0,Derivative[0,3][u][t,xMax]==0     (*FREE BOUNDARY CONDITIONS-START*)               },     u,     {t,0,TMax},     {x,xMin,xMax},     (*PrecisionGoal->Automatic,*)     Method->{"MethodOfLines","Method"->"LSODA","TemporalVariable"->t,"SpatialDiscretization"->{"TensorProductGrid","MinPoints"->10000,"MaxPoints"->10000,"DifferenceOrder"->1}}][[1]]      ## And the plotting function:               evapfbc2 = Plot[       uSolpbc[6981*100, x], {x, xMin, xMax},        PlotStyle -> {Black, Thick, Dashed},       AxesLabel -> {"\[Lambda]", "h"},       BaseStyle -> {FontWeight -> "Plain", FontSize -> 18},       PlotRange -> {{0, xMax}, {0, 3.5}}       ]      ## My questions are:   1. Should I be increasing the difference order for better results? (I understand that with a difference order set to `n`, it allows the stiff solver "elbow room" up to and including order `n`)   2. Is this the wrong method to use and is there something I could do in mathematica that would make this better? I am assuming that I would need to use a stiff solver for this problem. ## Note: With a difference order of 1, this runs for about 10-15 seconds on an i7/8GB RAM system. With a difference order of 12, it runs for about 20 **minutes** and consumes about 3-4 gigs of RAM.


title: How to display edge labels above edges?

text: The following code shows that edge labels "cover" edges in `GraphPlot3D`               GraphPlot3D[{{1 -> 2, "1\[Rule]2"}, 4 -> 1, {2 -> 4, "2\[Rule]4"},          1 -> 5, 2 -> 5, 5 -> 4},              EdgeRenderingFunction -> ({If[#3 =!= None,                  Text[#3, Mean[#1], Background -> Yellow], {}], Line[#]} &)]      But if I change `GraphPlot3D` to `GraphPlot`, the edges will go across the labels. How can I keep the labels with the highest z-index? Thank you. **Snapshot** `GraphPlot3D` ![enter image description here](http://i.stack.imgur.com/YHrfB.png) `GraphPlot` ![enter image description here](http://i.stack.imgur.com/k4Q3e.png)


title: Aligning a DensityPlot with a Plot

text: I am trying to show a gradient field behind a plot, and then its density in a separate plot below. I understand (as in this question) that this must normally be done using the ImagePadding directive. However, that does not quite align the plots in this case:               padding = {{20, 30}, {20, 30}};      GraphicsColumn[       {        Show[         DensityPlot[PDF[NormalDistribution[0, 1], Log[x/40]] ,            {x, 15, 60}, {y, 0, 30},            ColorFunction -> "SunsetColors"],          Graphics[{Dashed, Gray, Line[{{40, 0}, {40, 30}}]}],          Plot[If[x < 20, 10, If[x < 25, 2 x - 30, 20]], {x, 15, 60},            PlotRange -> {0, 30} ,            PlotStyle -> {RGBColor[0.5, .7, 0.2], Thickness[0.01]}],          ImagePadding -> padding, PlotRangePadding -> None        ],         Plot[PDF[NormalDistribution[0, 1], 5 Log[x/40]] , {x, 15, 60},            PlotStyle -> {RGBColor[1, 0.4, 0], Thickness[0.01]},           ImagePadding -> padding]       }     ]      ![Align me](http://i.stack.imgur.com/Nzm0g.png) Is there a way to get the abscissa ranges to match?


title: Random reals according to probability distribution

text: How do I generate a list a random real numbers that follow a probability distribution? If my pdf is $f(x)$, how do I "force" the random numbers to follow $f(x)$?


title: How to use ProbabilityDistribution for a posterior expression

text: I am performing a Bayesian updating on a Lognormal prior. After combining the prior and likelihood and dividing by the normalization constant I get the following expression as the posterior:               (0.805745 E^(-7.05*10^6 t - 0.255194 (15.576 + Log[t])^2))/t      I now want to use this as a distribution:               pdist = ProbabilityDistribution[posterior, {t, 10^-12, 10^-5}];      So far ok. I can now plot with the distribution               LogLinearPlot[PDF[pdist, t], {t, 10^-11, 10^-6}, PlotRange -> All, FrameLabel -> {"Failure Rate (f/hr)", "PDF"}]      But I cannot take the mean, Mean[pdist], or variance etc, they just return the input. Trying RandomReal[pdist] returns the following:               RandomReal::udist: The specification is not a random distribution recognized by the system      What am I missing in the assignment of the expression to ProbabilityDistribution[]. Clay


title: Random number generation from a given distribution

text: Is there any way to obtain random numbers from a user defined continuous distribution in Mathematica? Thank you very much


title: plotting hypergeometric functions

text: Does anyone know why a plot of a hypergeometric function turns out differently in _Mathematica_ than in _Maple_? The function I'm plotting between x=-30 and +30 is:               Re[Exp[-2 I p x]*(Conjugate[A]*Cosh[q - x]^(v + 1)*     Hypergeometric2F1[0.5 (v + 1 - I), 0.5 (v + 1 + I),       1/2, -Sinh[q - x]^2] +       Conjugate[B]*Cosh[q - x]^(v + 1)*Sinh[q - x]*      Hypergeometric2F1[0.5 (v + 1 - I) + 0.5, 0.5 (v + 1 + I) + 0.5,         3/2, -Sinh[q - x]^2])*(A*Cosh[q + x]^(v + 1)*       Hypergeometric2F1[0.5 (v + 1 + I), 0.5 (v + 1 - I),         1/2, -Sinh[q + x]^2] +       B*Cosh[q + x]^(v + 1)*Sinh[q + x]*       Hypergeometric2F1[0.5 (v + 1 + I) + 0.5, 0.5 (v + 1 - I) + 0.5,         3/2, -Sinh[q + x]^2])]      where `q=p=3`, `v=1`, `A=0.5+0.5I` and `B=-1+I`. I've attached both plots from _Maple_ (first image) and _Mathematica_ (second image) ![Maple plot](http://i.stack.imgur.com/SiiX1.jpg) ![Mathematica \[;pt](http://i.stack.imgur.com/5rreP.jpg)


title: Computing distance matrix for a list

text: Using functional programming in Mathematica, how can I compute a distance matrix for every element in a list of matrices... The distance would be computed between the item in the list and a "target matrix". For simplicity, let's consider that the list is a list of simple 3x2 matrices of integer with a range [0,4]. The target matrix is a matrix of that type as well. The reason I want to do that is to filter the list to find the matrix closest to the target one.


title: random number visualization without generator (Spectral Test)

text: Does anyone have codes for using only numbers to make this kind of random- number visualization and not requiring a number generator? http://demonstrations.wolfram.com/PoorStatisticalQualitiesForTheRANDURandomNumberGenerator/ http://en.wikipedia.org/wiki/Spectral_test


title: Pruning elements in an array based on the existence of similar elements in adjacent arrays

text: Imagine I have a series of particles floating around in $\mathrm{R}^3$, where, at any given time, particles can pop into and out of existence. Here, for some number of time points, $(t_1, \dots, t_N)$, I generate an unordered list of the positions of each of the particles present in the chamber: $(p_{(k,1)}, p_{(k,2)}, \dots) \in l_k$ s.t. I have a set of lists of particle positions, $(l_1, \dots, l_N) \in L$, each list corresponding to each time point. I would like to generate a set of pruned lists $(q_1, \dots, q_N) \in Q$, where each $q_i \subseteq l_i$ consists of the set of particle positions in $l_i$ that likely arose from measuring the same particle that   * existed and was measured at $\alpha$ previous time points. If, for example, we set $\alpha = 2$ and we measure a particle at time $t_i$, we would only place that measurement in $q_i$ if a very similar measurement was made at times $t_{(i-1)}$ and $t_{(i-2)}$.   * also existed and was measured for the next $\beta$ time points. If, for example, we set $\beta = 3$ and we measure a particle at time $t_i$, we would only place that measurement in $q_i$ if a very similar measurement was made at times $t_{(i+1)}$, $t_{(i+2)}$, and $t_{(i+3)}$.  If the particles were frozen in place, generating some $q_i \subseteq l_i$ would be a simple matter calculating a few intersections between $l_i$ (using the `Intersection[]` function in _Mathematica_ ) under the assumption that the probability one particle pops out of existence and another arises in exactly the same place in-between measurements is arbitrarily small. However, here the particles are moving, so we need to "guess" that some measured particle position at time $t_i$ and either times $t_{(i-\delta)}$ or $t_{(i+\delta)}$ (for small enough $\delta$) corresponds to the same particle if the Euclidean distance between the measured positions is below some threshold $\tau$. So if we set $\beta = 2$, we would need to see that the particle only moves a threshold distance $\tau$ from the current time point to the next, and then only moves a distance of at most $\tau$ in the interval between that time point to the one after (i.e., for some time point $t_k$ we require that some $\lvert p_{(k, \dots)} - p_{(k+1, \dots)}\rvert \leq \tau$ and that $\lvert p_{(k+1, \dots)} - p_{(k+2, \dots)}\rvert \leq \tau$). Said another way: Let $(l_1, l_2, \dots, l_N) \in L$ be an array of lists. For $\alpha = 1$ and $\beta = 1$, I'd like to remove all elements in each $l_i$ that are not within some threshold Euclidean distance $\tau$ of any element in the previous list, and any element in the next list in the array. If we set $\beta = 2$, I would like to prune away elements where there is no element, $p_j$ in list $l_{(i+1)}$ that is at most a Euclidean distance $\tau$ away, and where there is no element in list $l_{(i+2)}$ that is at most a Euclidean distance $\tau$ away from $p_j$. Provided a set of lists $(l_1, \dots, l_N) \in L$, and values for $\alpha$, $\beta$, and $\tau$, is there an elegant and fast way to generate my pruned lists $(q_1, \dots, q_N) \in Q$ (beyond writing a bunch of `Do` or `For` loops)? I am happy to handle the boundaries of the list, i.e., $l_1$ and $l_N$, in whatever manner is most convenient. ::: Let's generate a simple example ::: Let                l1 = {{87.02,892},{1,2},{65,23},{99,99}};      l2 = {{666,666},{777,777},{123,432},{876,28},{39,192},{87.01,892}};      l3 = {{87.0004,892},{8683,12},{7673,777},{7662,2120}};      l4 = {{40,40},{50,50},{87.002,892.001}};      Let $k = 3$, $\alpha = 2$, $\beta = 1$, and $\tau = 0.01$. Looking at list $l_2$, we notice that the only element in $l_3$ that is within a Euclidean distance $\tau$ of an element in $l_2$ is $p_0 = (87.0004, 892)$, displaced to the position $p_{(-1)} = (87.01, 892)$ in $l_2$. We now need to look at $l_1$ and check whether there is an element within a Euclidean distance of $\tau$ of $p_{(-1)} = (87.01, 892)$, and indeed we find that there is some $p_{(-2)} = (87.02, 892)$. Next, we need to check that there exists some $p_1$ in $l_4$ that is within a Euclidean distance $\tau$ of $p_0$, and indeed we find that $p_1 = (87.002,892.001)$ in $l_4$ satisfies this requirement. Looking through the rest of $l_3$, we find that no other elements satisfy the above set of requirements. Therefore, we can write $q_3 = (p_0) = ((87.0004, 892))$. We repeat this procedure to generate a pruned list $q_i$ for each $l_i$, and handle the boundary conditions in whatever manner is easiest. Please note that if we set $\tau = 0$, we could easily just write               q3 = Intersection[l1,l2,l3,l4]      To see this, try               q3 = Intersection[Round[l1], Round[l2], Round[l3], Round[l4]]      As expected, the result is               {{87, 892}}      If we want to grab the non-rounded value(s) in $l_3$ corresponding to this result, we can use the script               q3RealData = Flatten[Table[Nearest[l3, q3[[k]]], {k, 1, Length[q3]}], 1];      Which returns               {{87.0004,892}}      Dealing with $\tau > 0$ in the manner described above in an efficient way (rather than perhaps extending the rounding trick above) is the tricky part of this question.


title: Error/warning when using NSolve for simple equation

text: I am using NSolve to solve an equation, as shown here:               T == (0.000250709 E^(0.0269867/(0.0000734429 + 1. T)))/(1.00025 -          1. E^(0.0539734/(0.0000734429 + 1. T)))          NSolve[T == (0.000250709 E^(0.0269867/(0.0000734429 +                1. T)))/(1.00025 -           1. E^(0.0539734/(0.0000734429 + 1. T))), {T}, Reals]      It gives me the warning `Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>`, but still displays the result. Is it wrong to use NSolve for for such a general task, or am I missing a parameter?


title: Formatting results of integral approximation methods comparison.

text: The integral is `Integrate[Sqrt[1+x^3], {x, -1, 1}]`. I want to create a `Table` with three columns and four rows that shows the midpoint, trapezoidal, and Simpson's rule approximations of the integral with $n=4,8,16$ and $32$ subintervals. I defined the three approximations as:               f[x] := Sqrt[1 + x^3]          midptApprox[func_, {x_, a_, b_}, n_] := Module[{f, h = (b - a)/n},       f[t_] := func /. x -> t; h*Sum[f[a + i*h], {i, 0.5, n - 0.5}]]          trapex[func_, {x_, a_, b_}, n_] := Module[{f, h = N[(b - a)/n]},       f[t_] := func /. x -> t; (1/2)*h*         (f[a] + 2*Sum[f[a + i*h], {i, 1, n - 1}] + f[b])]          simpson[func_, {x_, a_, b_}, n_] := Module[{f, h = (b - a)/n, wt},       wt[i_] := 3. + (-1)^(i - 1); f[t_] := func /. x -> t;        If[EvenQ[n], (h/3)*(f[a] + Sum[wt[i]*f[a + i*h],                 {i, 1, n - 1}] + f[b]), "n must be even"]]      I was thinking to use `Grid` and `TableForm@Table[{n,f[n]}, {n,4,8,16,32}]` somehow but I don't know how to include the interval `[-1,1]`.


title: Put webMathematica folder outside Tomcat webapps folder?

text: Is it possible to put the `webMathematica` folder outside the usual Tomcat `webapps` folder and still have webMathematica work? _webMathematica_ does work OK if its whole folder is placed inside `webapps`. (This is with Tomcat 7 under OS X 10.8.3.) Configuration details would be greatly appreciated! In asking this I have two aims:   1. Granularize things so that I can just back up the webMathematica folder, along with other documents of mine, separately from backing up Tomcat's own files such as those in its distributed `webapps` folder.   2. Be able to access webMathematica (and, eventually, some specific apps using it) on localhost by means of a simple URL "alias", such as `webmathematica`, rather than the usual complete URL `localhost:8080/webmathematica` \-- just as I can now access Tomcat's webapps by entering the URL `tomcat.local`.


title: How to calculate this integral? Integrate[BesselJ[0, x - BesselJZero[0, 1]]/x, {x, -Infinity, Infinity}]

text: I tried to calculate the following integral, but it returned unevaluated.               In[1]:= Integrate[BesselJ[0, x - BesselJZero[0, 1]]/x, {x, -Infinity, Infinity}]     Out[1]= Integrate[BesselJ[0, x - BesselJZero[0, 1]]/x, {x, -Infinity, Infinity}]      Any ideas how to make Mathematica calculate it?


title: Does LinearModelFit perform an ordinary linear regression (least squares)?

text: Does `LinearModelFit` give an ordinary linear regression? I see lots of options, but nothing like "least squares" or OLR.


title: Finding n-th term of a sequence of functions

text: Lets say, I have the following series of functions, $f_1=a+1$   $f_2=a^2+2a+1$   $f_3=a^3+3a^2+3a+1$   Assume that $f_4, f_5,$ and $f_6$ are also known. However $f_n$ is not known as a function of $n$. From $f_1 ...f_6$, is it possible to reverse engineer the formula for $f_n$ using _Mathematica_? For the sample sequence given above, I should get, $f_n=(a+1)^n$. Is such a feature available in _Mathematica_? If so, can you provide links to the relevant documentation.


title: Manipulate doesn't work for plotting a region where a matrix is positive semi-definite

text: I'm facing the following problem and would very much like you help. Thanks is advance! I have defined the following 4x4 matrix: \begin{equation} S=\left( \begin{array}{cccc} \frac{1}{r^2} & 1 & 1 & t\sqrt{x}\\\ 1 & \frac{1}{t^2} & 1 & \frac{\sqrt{x}}{t}\\\ 1 & 1 & 1 & \sqrt{q} \\\ t\sqrt{x} & \frac{\sqrt{x}}{t} & \sqrt{q} & 1 \end{array} \right). \end{equation} where $(t,r)$ are two parameters (which I can choose arbitrarily to be any value in the segment $(0,1]$) and $(x,q)$ are the variables (each over the domain $[0,1]$). Now, for a fixed pair of $(r,t)$ I want to plot a region of $(x,q)$ over which the matrix $S$ is positive semi-definite (I want it to be a valid covariance matrix). The thing is that when I try to use the Manipulate function (manipulation on $(r,t)\in(0,1]^2$), and plotting the region of $(x,t)$, nothing shows. Namely, Iv'e written the following:               S={{1/r^2, 1, 1, t Sqrt[x]}, {1, 1/t^2, 1, Sqrt[x]/t},         {1, 1, 1, Sqrt[q]}, {t Sqrt[x], Sqrt[x]/t, Sqrt[q], 1}};          Manipulate[       RegionPlot[PositiveDefiniteMatrixQ[S],                   {x, 0, 1}, {q, 0, 1}], {r, 0, 1}, {t, 0, 1}]      and I get an empty graph no matter what values I set for $(r,t)$. You might say that maybe this matrix is never positive semi-definite, but the strange this is that if I manually input numerical values to $(r,t)$ and plot the region, everything works fine. I really need to manipulate over all possible parameter value of $(r,t)$. What is wrong here and how can I fix it? Moreover, is it possible for _Mathematica_ to parametrically find a local minima of a certain function (of the variable $q$ where $(r,t,x)$ serve as parameters) subject to the constraint that the above matrix must be positive semi-definite? If so, it would very much help me to know how to formalize such a command.


title: Find root approximants for low-precision numbers

text: I want to be able to find low-complexity algebraic approximants to decimal numbers. For example, 1.41 can be approximated as a solution to the polynomial `x^2==2`. There is a Mathematica function `RootApproximant` that seems like it should be able to do this approximation. However, I can't get it to work. `RootApproximant[N[Sqrt[2]]` returns $\sqrt{2}$ correctly, but `RootApproximant[1.41]` doesn't - it prefers to approximate it as `141/100`. Presumably I need to tell it to lower the accuracy it's aiming for. There is an option called "DegreeCost" that can be passed, but I can't work out how to use that to achieve my goal.


title: Solving PDEs with complicated boundary conditions

text: The system I'm trying to solve is $$\nabla^2 C_{(r,\theta)} =0$$ $$C_{(\infty,\theta)}=C_0$$ $$ [ \frac{\partial C_{(r,\theta)}}{\partial r} \cos(\theta)+\frac{1}{r}\frac{\partial C_{(r,\theta)}}{\partial\theta} \sin(\theta) ]|_{r=a*\cos(\theta)}=\frac{k}{D}C_{(r,\theta)}|_{r=a*\cos(\theta)}$$ and I have no idea how to enter the bottom equation into Mathematica. The problem I have is the inclusion of the boundary condition. How can I indicate to Mathematica that I want the derivatives to be taken and then have everything evaluated on the sphere $r=a*\cos(\theta)$ ?


title: Is there a way to separate variables between multiple notebooks?

text: I often correct homework by checking the calculations in Mathematica. Sometimes you would like to have two solutions open at once. However often defined symbols such as `f` or `phi` regularly overlap between the two notebooks.   * Is there a way to separate them other than using different symbols in every notebook?   * Is it possible to separate some variables yet share others between notebooks?


title: How to restrict global variables to within a given notebook?

text: > **Possible Duplicate:**   >  Is there a way to separate variables between multiple notebooks? I assigned a global variable like `a=10;` It works in all my notebooks, but I want to restrict that global variable to that notebook only. How to do that?


title: Given a sequence of expressions, how can I convert them all into functions?

text: I have a sequence of expressions               1, x, x^2, ..., x^5      I want to define a sequence of functions out of it, what should I do? Namely I want to define a vector of functions with elements: $$ f_1 = 1, f_2 = x, ..., f_6 = x^5 $$ Any help is appreciated.


title: Threads vs Cores: How many kernels need to launch in parallel Mathematica?

text: I have a dual-core intel machine. I would like to test how many parallel local kernels should be launch to get the optimized performance. Knowing that the HyperThreading technology of Intel could run 2 "Threads" on 1 core. That means we could run 4 threads on a dual-core CPU. For informations about "Threads vs Cores", look at here.. AMD says "no" to HyperThreading, and after them, cores is much more important. view amd vs intel: " _Based on the results of the Cinebench® and PovRay® benchmark tests, we find that physical cores outperform enhanced threading software in certain applications."_ Function to evaluate:               mlength[z_] :=       Length[FixedPointList[#^2 + z &, z, 20, SameTest -> (Abs[#] > 2 &)]]      **1\. Sequential Evaluation (No parallel)**               CloseKernels[]; Kernels[]     ParallelTable[        mlength[x + I y], {y, -1, 1, 0.005}, {x, -2, 1,          0.005}]; // AbsoluteTiming      > During evaluation of ParallelTable::nopar: No parallel kernels available; > proceeding with sequential evaluation. > > {8.814056, Null} ![enter image description here](http://i.stack.imgur.com/aOIVc.png) Although no parallel kernel is used, we can see that the second core on the right shares about 50% the work with the first core. I think Windows automatically shares the work for multicore. **2\. Using 1 Slave kernel**               CloseKernels[]; LaunchKernels[1]; Kernels[]     ParallelTable[        mlength[x + I y], {y, -1, 1, 0.005}, {x, -2, 1,          0.005}]; // AbsoluteTiming      > {9.110458, Null} ![enter image description here](http://i.stack.imgur.com/zb4xB.png) By asking Mathematica add one slave kernel, this slave kernel runs on the first core mainly (on the left) and we can see that the work is mainly done in the first core. There are still some works on the second core. **3\. Using 2 slaves kernel**               CloseKernels[]; LaunchKernels[2]; Kernels[]     ParallelTable[        mlength[x + I y], {y, -1, 1, 0.005}, {x, -2, 1,          0.005}]; // AbsoluteTiming      > {4.820431, Null} ![enter image description here](http://i.stack.imgur.com/D6Ffm.png) Now, with 2 slaves kernels, the work is evenly distributed on 2 cores. That's the reason why, we get the results done in half of the time. We can see that Mathematica parallelizes better than Windows. **4.Using 4 slaves kernel**               CloseKernels[]; LaunchKernels[4]; Kernels[]     ParallelTable[        mlength[x + I y], {y, -1, 1, 0.005}, {x, -2, 1,          0.005}]; // AbsoluteTiming      > {5.023232, Null} ![enter image description here](http://i.stack.imgur.com/pAeic.png) ![enter image description here](http://i.stack.imgur.com/Tm0qP.png) I'm trying to add more kernels by supposing that, one core could maintain 2 "threads". But the results is slower. I think Multithreading does not help in Mathematica. So the rule of thumb is: the number of parallel local kernels = number of cores (not threads). Increasing number of kernel just increases the overhead of data transfer between kernels. Any suggestions? ================================================================= **Udate to another computer Intel i7, 4 CORES, 8 THREADS (TRUE HT).**               maxKernel = 8;     data = Table[       CloseKernels[];         LaunchKernels[n];         ParallelTable[          mlength[x + I y], {y, -2, 1, 0.002}, {x, -2, 1, 0.002}]; //          AbsoluteTiming // First, {n, 0, maxKernel}]      > {40.332307, 44.634553, 22.154267, 17.925025, 14.795846, 13.280760, \ > 12.369708, 13.218756, 12.701726} ![enter image description here](http://i.stack.imgur.com/fYOo7.png) We can see that, the time decreases up to 4 Cores, and then, up to 6 threads, time optimized only 1-2 second.


title: Suppressing Null output with sparse data? A lot of zeros and little nonzeros

text: I have sparse situation with a lot of zeros, related chat. I want to see only the nonzero things or nonnull things. **Example** > **Input** >   >   >     hello[ii_, jj_] := If[ii == 1, Return[True], Return[False]] >     Table[If[hello[ii, jj], {ii, jj} -> 1], {ii, 1, 5}, {jj, 1, 5}] >   > > **Output** >   >   >     {{{1, 1} -> 1, {1, 2} -> 1, {1, 3} -> 1, {1, 4} -> 1, {1, 5} -> >        1}, {Null, Null, Null, Null, Null}, {Null, Null, Null, Null, >       Null}, {Null, Null, Null, Null, Null}, {Null, Null, Null, Null, >       Null}} >   > > **Intended Output** >   >   >     {{1, 1} -> 1, {1, 2} -> 1, {1, 3} -> 1, {1, 4} -> 1, {1, 5} -> 1} >   where the full form of the intended output could be a 5x5 matrix with 1st row with ones and otherwise zeros.


title: Maintaining map dimensions after using Colorize

text: I am trying to beautify a map by using `Colorize` as demonstrated in the _Mathematica_ documentation. First, I import a shapefile from DIVA-GIS for France which contains polygons for the country, its regions and departments.               mapdata = Import["http://gadm.org/data/shp/FRA_adm.zip", "Data"];      I then extract the region data and (poorly) place a bounding region around the map:               Graphics[Rectangle[{-6, 41}, {9.7, 51.5}]];     Graphics@{EdgeForm[Thick], FaceForm[White], ("Geometry" /. mapdata[[2]])};     t1 = Show[%%, %]      I can then `Colorize`               t2 = Colorize[MorphologicalComponents[t1, CornerNeighbors -> False], ColorFunction -> "LightTerrain"]      which gives me the following `Image` (note, it is not a `Graphic`) ![France in terrestrial colors](http://i.stack.imgur.com/1CNXP.png) I'd now like to take the department information and lay it over the colored regional map               t3 = Show[Graphics@{Black, EdgeForm[Thin], FaceForm[], ("Geometry" /. mapdata[[1]])}]      ![enter image description here](http://i.stack.imgur.com/tbq72.png) However the `MorphologicalComponents` function rescaled the regional map, preventing me from overlapping t2 and t3 above. Is there a way to get the same or similar effect as obtained through `MorphologicalComponents` and `Colorize` without losing the x,y coordinate data?


title: Elegant operations on matrix rows and columns

text: ## Question The Mathematica tutorial has a section 'Basic Matrix Operations', describing operations like transpose, inverse and determinant. These operations all work on entire matrices. I am missing a section on basic operations on matrix rows / columns. For example:   1. Extracting a row from a matrix   2. Inserting a row into a matrix   3. Adding two rows within a matrix together   4. Swapping two rows   5. Multiplying a row with a number And similar for columns. What is the most elegant way to implementation of these operations? Speed is not important for me, but simplicity is. ## Summary Here I summarize my personal taste. I will update it whenever someone suggests a way I like more.               m = Range@12 ~Partition~ 3;     m // MatrixForm      $\begin{pmatrix} 1 & 2 & 3 \\\ 4 & 5 & 6 \\\ 7 & 8 & 9 \\\ 10 & 11 & 12 \end{pmatrix}$ Insert a column at position 2:               v = Range[21, 24];     Insert[m // Transpose, v, 2] // Transpose // MatrixForm      $\begin{pmatrix} 1 & 21 & 2 & 3 \\\ 4 & 22 & 5 & 6 \\\ 7 & 23 & 8 & 9 \\\ 10 & 24& 11 & 12 \end{pmatrix}$ ### Extract row / column Extract row 2:               m[[2]]      $(4,5,6)$ Extract column 2               m[[All, 2]] // MatrixForm      $\begin{pmatrix}2\\\5\\\8\\\11\end{pmatrix}$ ### Insert a row / column Insert a row at position 2:               v = Range[13, 15];     Insert[m, v, 2] // MatrixForm      $\begin{pmatrix} 1 & 2 & 3 \\\ 13 & 14 & 15 \\\ 4 & 5 & 6 \\\ 7 & 8 & 9 \\\ 10 & 11 & 12 \end{pmatrix}$ ### Adding two rows / columns column 3 = column 3 + column 1:               m2 = m;       m2[[All, 3]] += m2[[All, 1]];     m2 // MatrixForm      $\begin{pmatrix} 1 & 2 & 4 \\\ 4 & 5 & 10 \\\ 7 & 8 & 16 \\\ 10 & 11 & 22 \end{pmatrix}$ row 2 = row 2 + row 3:               m2 = m;     m2[[2]] += m2[[3]];     m2 // MatrixForm      $\begin{pmatrix} 1 & 2 & 3 \\\ 11 & 13 & 15 \\\ 7 & 8 & 9 \\\ 10 & 11 & 12 \end{pmatrix}$ ### Swapping rows / columns Swap row 1 and row 3:               m2 = m;     m2[[{1, 3}]] = m2[[{3, 1}]];     m2 // MatrixForm      $\begin{pmatrix} 7 & 8 & 9 \\\ 4 & 5 & 6 \\\ 1 & 2 & 3 \\\ 10 & 11 & 12 \end{pmatrix}$ Swap column 1 and 3:               m2[[All, {1, 3}]] = m2[[All, {3, 1}]];     m2 // MatrixForm      $\begin{pmatrix} 3 & 2 & 1 \\\ 6 & 5 & 4 \\\ 9 & 8 & 7 \\\ 12 & 11 & 10 \end{pmatrix}$ ### Multiplying rows / columns Multiply row 2 with 2:               m*{1, 2, 1, 1} // MatrixForm      $\begin{pmatrix} 1 & 2 & 3 \\\ 8 & 10 & 12 \\\ 7 & 8 & 9 \\\ 10 & 11 & 12 \end{pmatrix}$ Multiply column 1 with 5:                ((m // Transpose)*{5, 1, 1}) // Transpose // MatrixForm      $\begin{pmatrix} 5 & 2 & 3 \\\ 20 & 5 & 6 \\\ 35 & 8 & 9 \\\ 50 & 11 & 12 \end{pmatrix}$ ## References   * What is the most efficient way to add rows and columns to a matrix?   * Thanks to nikie for suggesting Matrix and Tensor Operations tutorial   * Chris Degnen pointed out http://stackoverflow.com/questions/7537401/how-to-insert-a-column-into-a-matrix-the-correct-mathematica-way


title: Construct matrix B from columns of matrix A?

text: Starting with Mathematica and working with matrices, I am looking for a way to construct a new matrix B from columns of matrix A. For example, matrix a has 8 columns and 4 rows, I want matrix B to be columns 1,2, 7 and 8. Got this so far:               A12 = A[[All, 1;;2]];     A78 = A[[All, 7;;8]];      ... but have not yet found a way to properly join these two... I did find a way to insert a column in a matrix in the Mathematica documentation, by using a nested Transpose-Insert-Transpose, but did not return the correct new matrix B. Any suggestion is much appreciated. J.


title: Manipulating lists with constant values

text: let us assume we have some data in the following form.               Table[data[[1,i,{1,2}]],{i,1,Length[data[[1]]]}]      This would give all $x-y$ pairs of the first file. What I want to do now is to multipliy this table with constant factors but different factors for $x$ and $y$. I found a solution for this with               MapThread[ Composition[Flatten, List], {xlist, #}] & /@ ylist      where `xlist` is                Table[data[[1,i,1]],{i,1,Length[data[[1]]]}]*c1      and `ylist` is                {Flatten[Table[data[[1,i,2]],{i,1,Length[data[[1]]]}]]}*c2      Is there a more comfortable way to do this? Thank you in adavance Sincerely * * * EDIT:                test= {{0.0015856, -1486.76}, {0.00157776, -1483.45}}      I need :                test2= {{1.5856, -148.676}, {1.57776, -148.345}}      So the $x$ values should be multiplied with $1000$ and the $y$ values multiplied with $0.1$ EDIT2: Transpose is exactly what I need, thank you very much


title: Convenient methods to sum select rows (or cols) of a matrix?

text: What are some convenient methods to sum selected rows or columns in a matrix? This question Summing along rows or columns of a matrix only addresses how to sun entire rows or columns ie, by level. For example, suppose one wants to add columns 2 and 4 of:               M = Column@Table[RandomInteger[10], {3}, {4}];      Neither MapAt, nor ReplacePart appear useful in this regard. Of course `#[[2]] + #[[4]] & /@ M` works for summing columns and `Transpose` can be used for rows, but that seems like a hack to specify using such syntax. Any better solutions?


title: How to extract data from a list?

text: I have a long list of coordinates defined in some 8D space. I need to extract the first three coordinates only. How ? For example :               coords = {{0.1, 0.4, -0.2, 1.2, 3.2, 0.5, 0.8, 2.2}, ...};      I need to define the list of all first three components :               reducedCoords := {{0.1, 0.4, -0.2}, ...}      Please, notice that I'm working with _Mathematica 7_.


title: Swap elements in list without copy

text: I have a lists with 100,000 elements and would like two swap to entries without dublicating the list Example:               a = {11,2,7,31};     swap[2,4]; (* how to implement this? *)      a (* should print {11,31,7,2} *)      I can do               swap[i_,j_]:=(a=ReplacePart[a,{i->a[[j]],j->a[[i]]}])      but this copies the list and if there are many swaps this takes way too much time... **Edit:** The important point of my question is how to do it without having Mathematica creating a copy of the list. This is in no way answered by the alleged duplicate.


title: How to switch elements i,j in a list?

text: Are there any direct function to switch elements i,j in any list? or is necessary to create one with Do / If ? (I don´t ask for making me this function if it doesn´t exists) My target is can switching positions i,j in {1,4,3,2,5} for example I tell switch(2,5) an obtain {1,5,3,2,4}


title: List operation on specific elements

text: Let's consider two tables given by               T1 = Table[{i,j,f1[i,j]},{i,0,2},{j,0,3}];     T2 = Table[{i,j,f2[i,j]},{i,0,2},{j,0,3}];      The functions `f1` and `f2` are expensive to compute, so that I want to compute the tables `T1` and `T2` only one time. I now want to obtain the table `T3` defined as               T3 = Table[{i,j,f1[i,j] f2[i,j]},{i,0,2},{j,0,3}];      However, as I already know `T1` and `T2`, I want to estimate `T3` without any reevaluation of `f1` and `f2`, but using only operations on the already known tables `T1` and `T2`. How would you perform such an operation on the lists ?


title: Arithmetic operations on columns

text: I have imported a file, data.txt, then plotted columns 2 and 3 as function of column 1 using the commands below.               data = Import["/home/data.txt", {Data}]          ListlogLogPlot[       Transpose /@ {{data[[All, 1]], data[[All, 2]]}, {data[[All, 1]], data[[All, 3]]}}]      Now I have to add/subtract column 1 and column 2, then plot them together. Please help me to proceed


title: Column Operation - Subtract part of nested list from part of another nested list

text: I have two lists with x- and y-coordinates with different z-coordinates of the form               {{x0, y0, z0_1}, {x1, y1, z1_1}, ..., {xn, yn, zn_1}}        and               {{x0, y0, z0_2}, {x1, y1, z1_2}, ..., {xn, yn, zn_2}}.        I want to construct the list > >     {{x0, y0, z0_1 - z0_2}, {x1, y1, z1_1 - z1_2}, ..., {xn, yn, zn_1 - > zn_2}}   >   i.e. construct a list containing the difference of only the z-coordinates, for each x- and y-coordinate. Is it possible to do this with a few simple commands, or should I construct a loop and create a new list from scratch?


title: Subtracting list from list of same dimensions

text: This might seem like a basic operation I have already searched for similar answers but found none. I want to subtract the second elements of my list ct from the corresponding element in lt, and letting the first element to remain as they are the same.               lt={{{1, 1}}, {{2, 1}}, {{3, 1}}, {{2, 3}}, {{5, 1}}, {{2, 4}, {3, 2}}, {{7, 1}}, {{2, 7}}};               ct = {{{1, 1}}, {{2, 1}}, {{3, 1}}, {{2, 2}}, {{5, 1}}, {{2, 1}, {3, 1}}, {{7, 1}}, {{2, 3}}};      as you can see the first elements are the same I want them to remain. desired result               nt={{{1, 0}}, {{2, 0}}, {{3, 0}}, {{2, 1}}, {{5, 0}}, {{2, 3}, {3, 1}}, {{7, 0}}, {{2, 4}}};      also note that my resulting list nt is my list lt with the last previous first element. like element 4 in my nt is just my element 2 in lt. and where the second element is 0 that means there is no previous first element of the same value, that is a next way to look at it.


title: Efficient priority queues?

text: I'm trying to figure out how to transfer the concept of a priority queue to the functional world. Searches have turned up some implementations that use Append and other expense list copying techniques. I'm guessing there is a better way. An example of what I am trying to solve is consider the products of all pairs of N digit numbers in descending value order. For small N I can do something like...               Reverse[         Cases[             SortBy[                 Flatten[                     Table[{i, j, i*j}, {i, 1, 9}, {j, 1, 9}],                 1],              Last[#] &],          {i_, j_, k_} /; i <= j]     ]      Alternative solutions to the problem in particular are welcomed, but I am really looking for a generic answer of how to apply the priority queue concept to the functional world.


title: How do I plot a system with gain, first-order lag and time delay?

text: I want to plot the output response of a system with gain, first-order lag and time delay, where the transfer function is               tf = 2 Exp[-0.1 s]/(.1 + s)      To plot a time delay, we need to use an approximation, as _Mathematica_ 's `TransferFunctionModel` is restricted to rational transfer functions. We can use a Padé approximation               pa = PadeApproximant[tf, {s, 0, 5}]      and then determine the output response and plot it using               model = TransferFunctionModel[pa, s]     output = OutputResponse[model, UnitStep[t], t];     Plot[output, {t, 0, 10}, PlotRange -> All]      My output response looks very weird: ![time delay plot](http://i.stack.imgur.com/h9f6F.png) Can somebody please explain to me what I'm doing wrong?


title: Draw generalized k-ellipse using "the locus of points so that sum of the distances to the foci is constant" deffinition

text: My goal is to draw a family (for different distances) of _3-ellipse_ for equilateral set of foci. _3-ellipse_ definition is: > "the locus of points so that sum of the distances to the three foci is > constant" Eventually I want to generalize this to **k-ellipse**. Thanks for your help I am bloked.


title: How can I highlight and individually color the connected components of a graph?

text: I have a graph `G`, and I wish to individually highlight and color the different connected components of a graph. It is simple enough to do this for one component:               HighlightGraph[G, ConnectedComponents[G][[2]]]      But how does one do this for multiple components? I'd ideally like to be able to pass a color string to specify each highlighted component, something like `componentColorString = { Blue, Green, Red,...}`. Is this possible?


title: Permutations with restrictions

text: I have a problem that can be modeled as: there are 8 boxes in total, and 60 different items. I want to put all items into boxes (all of them can be put into a single box). Now I want to find: 1) the number of total posssible situations 2) the number of posssible situations when certain number of items are put into certain box (such are 6 different items are put into box #2). **Edit:** I found a somehow similar problem, but I think my data would be too large to enumerate. http://stackoverflow.com/questions/20721161/combinations-between- two-lists


title: Append CellTag to list of Cells

text: As suggested in this earlier post, I restate my question here about appending `CellTags` option to existing list of `Cells`. What I have for the moment:               listOef = {};          SelectionMove[this, Before, Notebook];     SelectionMove[this, Next, Cell];     SelectionMove[this, Next, Cell];     Do[      (       SelectionMove[this, Next, CellGroup];       AppendTo[listOef, NotebookRead[this]]       )       , {i, aantalVragen}]      This copies cells I give in the `EvaluationNotebook[]` to the list `listOef`. I now want to write this list to the new notebook `nb` after permuting its elements randomly with               allPermutations = Table[RandomPermutation[aantalVragen], {aantalToetsen}];      like this:                Do[       (NotebookWrite[nb,         Permute[listOef, allPermutations[[toetsNummer]]][[vraagNummer]]];        ),      {vraagNummer, aantalVragen}]      In this `Do[]` command, I would like to add Labels to every cell during the `NotebookWrite` command: its label would be `"vraagNummer"`, the counter which goes from `1` to `aantalVragen`, which is given at the beginning of the document. I'd like to have something like this:               SetOptions[ Cell, CellLabel->"vraag"<>ToString[vraagNummer] ]      but this doesn't seem to work. How can I address the cell which has just been written by `NotebookWrite`? Or how would you suggest to do this using this structure? Thanks for all help!


title: Algorithm in RandomChoice[{w1,w2,...}->{e1,e2,...},n]?

text: I wonder if someone happens to know the actual algorithm implemented in the particular overload of `RandomChoice` that handles explicitly presented discrete distributions, namely               RandomChoice[{w1, w2, ...}->{e1, e1, ...}]         gives a pseudorandom choice weighted by the w1.      It's blazing fast, as illustrated in my test notebook here: https://www.dropbox.com/s/rwm80fut60v926b/FastNonUniformPseudoRandoms002.cdf (btw, such things are really useful in, say, particle filters). I know of at least four algorithms for solving this problem. A greedy-space array of outcomes has O(1) time but O(S) space where S is the sum of the frequency counts (and only works when the weights are integers or can be scaled up to integers). Two algorithms explicitly invert the CDF, either by linear search -- O(N) time, O(N) space, or binary search -- O(log N) time, O(N) space. The best I know is Walker's Alias Method http://en.wikipedia.org/wiki/Alias_method which is O(N) space, O(1) time. I implemented this in my notebook and it beats the pants off the others (as expected) but is still 30 times slower than the built-in (though still just O(1)). I wonder if the built-in is just an optimized Walker's or is it something better?


title: ParallelTable and Precision

text: I'm using `ParallelTable[]` to calculate a function over a range of my parameters , ($\omega,\ell$). This seems to be working well (in terms of speed increase) except for some strange Precision issues. I set the `$MinPrecision=40` at the start of the notebook as well as my working precision, precision goal and accuracy goal for some computations later in the notebook as                wp=$MinPrecision;      ac=$MinPrecision-8;      pg=wp/2;      The rest of my code basically defines some helper functions before going on to the main function (that I feed into `ParallelTable`). This main function uses the helper functions along with `NDSolve` and `NIntegrate` to do some computations, and it's in these that I feed the `Working Precision->wp` ,`AccuracyGoal->ac`, `PrecisionGoal->pg`. I have used `DistributeDefinitions` on all my variables, helper functions, and the main function, and even on `$MinPrecision`, also all my initial conditions for NDSolve have `N[...,wp]` wrapped around them, but still when I use `ParallelTable[...]` I get errors:                Precision::precsm. Requested precision 39.153977328439204` is smaller than $MinPrecision.       I don't get any such error when I simply use `Do[...]` or `Table[..]` on my main function, and indeed the results differ at the 33rd decimal place, which is the decimal the above number is precision too. I have no idea what this number is unfortunately, it doesn't look like any of my outputs or inputs. I just don't understand how this could not crop up with the non-parallelized forms, but crop up with parallelized? **Minimal Example** This is much simpler than my code, but I think it still captures it and the problem: Definitions:               M = 1;     $MinPrecision = 40;         wp = $MinPrecision;     ac = $MinPrecision - 8;     pg = wp/2;     rinf = 15000;      The main function to be parallelized:               dGenBessE[\[Omega]_?NumericQ, l_?IntegerQ] :=      Block[{\[CapitalPhi]out, init, dinit},          init = 0.00006630728036817007679447778124486601253323 +      6.913102762021489976135937610105907096265*10^-6 I;     dinit = -6.958226432502243329110910813935678705519*10^-7 +      6.631148430876236565520382557577187147081*10^-6 I;            \[CapitalPhi]out[\[Omega], l] = \[CapitalPhi] /.        Block[{$MaxExtraPrecision = 100},            NDSolve[{\[CapitalPhi]''[r] + (2 (r - M))/(                r (r - 2 M)) \[CapitalPhi]'[                 r] + ((\[Omega]^2 r^2)/(r - 2 M)^2 - (l (l + 1))/(                  r (r - 2 M))) \[CapitalPhi][r] ==               0, \[CapitalPhi][rinf] ==               N[init, wp], \[CapitalPhi]'[rinf] ==               N[dinit, wp]}, \[CapitalPhi], {r, 30, 40},             WorkingPrecision -> wp, AccuracyGoal -> ac,             PrecisionGoal -> pg, MaxSteps -> \[Infinity]]][[1]];          Print["For \[Omega]=", \[Omega], " and l=", l, ": "];          Print["Kernel ID: ", $KernelID];      Print["Precision of init: ", Precision[init]];      Print["Precision of dinit: ", Precision[dinit]];      Print["\[CapitalPhi]out at 35=" ,      N[\[CapitalPhi]out[\[Omega], l][35], wp] ];      Print["Precision of \[CapitalPhi]out: ",       Precision[\[CapitalPhi]out[\[Omega], l][35]]];     ]      Do Paralleization prereqs:               LaunchKernels[4]     DistributeDefinitions[M, rinf, wp, ac, pg, $MinPrecision,dGenBessE];      Attempt to run it with Paralleize for two different values:               Parallelize[{dGenBessE[1/10, 0], dGenBessE[1/10, 1]}] // AbsoluteTiming      **Result** For me this leads to a               Precision::precsm: Requested precision 38.95475956978393` is smaller than   $MinPrecision. Using $MinPrecision instead.


Is there are any command to expand quaternion multiplication as below STEP by STEP.   Below is equation from my text.   I'm self learner and for some reason I can't wrap my head around quaternion multiplication.   I really need to see step by step workout.   I remember seeing step by step workout in wolframalpha but that's for calculus integration.   Just wondering if the same can be done in mathematica.   $$ \begin{equation} \begin{split} w =& qvq^* \\\ =& (q_0+\vec{q})(0+\vec{v})(q_0-\vec{q}) \\\ =&(2q_0^2-1)\vec{v}+2(\vec{q}\cdot\vec{v})\vec{q} + 2q_0(\vec{q}\times\vec{v}) \end{split} \end{equation} $$


title: Why can't I pass in a product of matrices in my modified matrixform?

text: In attempt to build a function that would work like `MatrixForm` but would wrap a matrix with brackets, I tried the following function:               matrixform[mat_] := TraditionalForm[DisplayForm[RowBox[{     StyleBox["[", SpanMaxSize -> \[Infinity]],      GridBox[mat],       StyleBox["]", SpanMaxSize -> \[Infinity]]}]]]       This does what I want on an ordinary matrix, say `a={{1,2},{3,4}}`, as well as on `2a` and `a+a`. But if I pass in a product such as `a.a` it does not work properly. I have tried wrapping the product in `Evaluate`, and tried executing `a.a` separately and then passing in `%` and this has made no difference. Note that a better function `Composition[bracketF,Grid]` suggested here lacks these issues. So the question is why can we not pass in a product of matrices into `matrixform`? This is my first time working with functions such as `StyleBox, GridBox, RowBox`, etc.


title: Slider with varying scales

text: I am trying to implement a slider whose control range varies depending on the current value:               {Dynamic@asd,      Slider[Dynamic[asd],      Dynamic @Which[asd < 10, {1, 10}, asd >= 10 , {9, 50, 1}]]}      However I now have the problem that the slider remains in place once the control scales change and thus the value "leaps" as soon as the scales change. To Illustrate: once i reach ten the scale changes but the slider remains in place and thus the value goes immediately to 50 (the "new" ten). Is it possible to have the slider jump instead of the value? Perhaps by briefly locking the value at ten?


title: plot Contourplot simple

text: demElev[path_, samples_] :=       Module[{elevUrl, x, y, markers, dem},        elevURL =         "http://maps.google.com/maps/api/elevation/json?sensor=false"; {x,          y} = path[[1]] - path[[2]];        pnts = Table[         path[[2]] + {i, j}, {i, 0, x, x/samples}, {j, 0, y, y/samples}];        markers =         StringReplace[           ToString[#], {"}, {" -> "|", "{" -> "", "}" -> "",             " " -> ""}] & /@ pnts;        dem = ToExpression@           StringCases[Import[elevURL <> "&locations=" <> #],             NumberString] & /@ markers; Partition[Flatten[dem], 4]]     path = {{41.334387, 104.574955}, {51.02568, 120.563405}}     dem = demElev[path, 120]          ListContourPlot[dem[[All, {3, 2, 1}]],       ContourStyle ->        Directive[GrayLevel[0], Opacity[0.946], AbsoluteThickness[2.5`]],       Contours -> {688, 1101}, PlotLegends -> Automatic,       ContourShading -> None, ColorFunction -> "Thermometer",       Epilog -> {AbsolutePointSize[10]}]      How to draw contour plot not a sharp angle this picture?? 2D![enter image description here](http://i.stack.imgur.com/Owq6o.png) ![enter image description here](http://i.stack.imgur.com/vtCBj.png)


title: Table with multiple iterators running parallel

text: I want to use `Table` to generate a list, but it only allows one iterator to increment at a time. So in               Table[i*r, {i, 0, 5}, {r, 0, 5}]      `r` will run `0` through `5` for `i = 0`, then change `i` and run through `r` again etc., giving 36 elements in a 6 x 6 array. But I want a list of only 6 elements, so that `i` and `r` are changed at the same time, in parallel. So it would give `0*0`, then `1*1`, `2*2` etc. Giving the output of:               {0, 1, 4, 9, 16, 25}  (* <-This is what I want. *)      Is this possible with `Table` or anything else? Also I know I could do:               Table[i*i, {i,0,5}]      But this is just an example, and in my project this way won't work. Also if at all possible I would like to use a list for one iterator's values, but this isn't so important.


title: Solving/Reducing equations in $\mathbb{Z}/p\mathbb{Z}$

text: I was trying to find all the numbers $n$ for which $2^n=n\mod 10^k$ using _Mathematica_. My first try:               Reduce[2^n == n, n, Modulus -> 100]      However, I receive the following error:               Reduce::nsmet: This system cannot be solved with the methods available to Reduce. >>      using $n^2$ instead of $2^n$ works just fine, where is the problem with $2^n\;$? On top of that, how do I keep the modulus $10^k$ variable and therefore, receive a solution dependent on $k\;$?


title: Verifying and deriving basic (block) matrix identities

text: How can I use the new symbolic matrix/tensor capabilities to verify matrix identities, such as ![enter image description here](http://i.stack.imgur.com/rbU71.gif) (1) or ![enter image description here](http://i.stack.imgur.com/YwVpc.gif) (2) Even better, how can I ask Mathematica to derive expressions for X, Y, Z, and U like the ones shown on the RHS of (2) (many equivalent forms exist) (and as a bonus, display them similarly).