Hot answers tagged

16

I am delighted by this problem mostly because I was not aware of the underlying physical model of phase separation (the Cahn–Hilliard equation)! Anyway, here is an approximation of a somewhat similar behavior to start the discussion: NestList[ Sharpen@Threshold[GaussianFilter[#, 2], {"Hard", {"BlackFraction", 0.1}}] &, img, 50 ]; ...


12

Your code shows you are somewhat confused about the indices. Here is some code that is more Mathematica idiomatic and which makes keeping the indices straight much easier. I am running the simulation with a 24 hour step to cut down the data plotted. Should work just a well with your time step. maxk = 10; mink = 2; steps = 365; batterylevel = ...


12

Let's say that both players start with 10 dollars. We can represent the game state with a list: start = {10, 10}; Create a function which plays one round of the game. There are only two possible outcomes, we can use RandomChoice to pick randomly between them: oneRound[{x_, y_}] := RandomChoice[{{x + 1, y - 1}, {x - 1, y + 1}}] For example ...


11

Taken from the Matlab code here http://www.math.utah.edu/~eyre/computing/matlab-intro/ch.txt with slight modifications: m = n = 256; delx = 1/(m - 1); delx2 = delx^2; x = Range[0, 1, delx]; delt = 0.00000001; ntmax = 250; epsilon = 0.01; eps2 = epsilon^2; a = 2; lam1 = delt/delx2; lam2 = eps2*lam1/delx2; leig = Transpose@Table[2 Cos[π Range[0., n - ...


11

In Mathematica it is natural to approach such a task with list operations and pattern matching. dice1 = RandomInteger[{1, 6}, 100]; dice2 = RandomInteger[{1, 6}, 100]; Count[dice1 + dice2, 2 | 3 | 4 | 5 | 6] You seem to be a very new beginner, since you are using x[i] and y[i] as if these are vectors, when they are in fact not, in Mathematica. Mathematica ...


11

There are many ways you can do this, e.g. ri = RandomInteger[{1, 6}, {100, 2}]; SortBy[Normal@GroupBy[ri, Total, Length@#/100. &], First] yielding: {2 -> 0.01, 3 -> 0.05, 4 -> 0.11, 5 -> 0.08, 6 -> 0.13, 7 -> 0.12, 8 -> 0.17, 9 -> 0.14, 10 -> 0.11, 11 -> 0.07, 12 -> 0.01} rules linking sum to frequency. You can also exploit ...


11

Here's a version using CellularAutomaton to simulate the Ising model (with Glauber dynamics), which can be viewed as a lattice version of the Cahn–Hilliard equation. It's not quite as smooth as the GIF, but it can probably be improved with some tweaks: With[{L = 100, T = 0.1, t = 100, dt = 10, r = 0.1}, ArrayPlot[#, ColorRules -> {-1 -> Black, 1 ...


10

Here is an example with a force that pushes down between the stretch defined by 2<=x<=3 && y==1 and a fixation at the point y==1 && x==5. Applying a force over a stretch is done by a NeumannValue and setting a fixed displacement on a point is done by a DirichletCondition. {uif, vif} = NDSolveValue[{planeStrainOperator[10^3, 33/100] ...


8

There is an example in the documentation which may get you started: pairs = RandomReal[{-1, 1}, {10000, 2}]; 4 Count[Map[Norm, pairs], _?(# <= 1 &)]/10000. (*3.1248*) You can plot see this as: Graphics[{PointSize[Small], Blue, Point@Select[pairs, Norm[#] <= 1 &], Gray, Point@Select[pairs, Norm[#] > 1 &], Red, Thick, Circle[]}, ...


7

The Finite Element solver in Mathematica does run in parallel, both element computation and the linear solve process are spread over the CPU cores available. Additionally, the option "MeshElementBlocks" for ToElementMesh splits the mesh elements in blocks which could be used for a domain decomposition. To get a more detailed answer you'd need to clarify ...


7

Here is a modified version of your code. On my PC it completes your example run in about 5 seconds. I won't try to describe every change but will point out the major features. Some of the changes are stylistic rather than performance-based. This is not a criticism of your style but a reflection of the way I broke the original code down in order to ...


7

When building a simulation like yours you should test the performance of the individual components before incorporating them into the simulation. That is, you should know the cost of the components as well as the values they return. Here is a simple example based on your code. You use Clip in a couple of places to limit values from below. That is suspect, ...


7

While the answers so far have covered a lot of ground already I have not seen EmpiricalDistribution. I would like to build upon this observation by providing a couple of general considerations that I have found to be useful when doing statistical experiments using Mathematica. What users of Mathematica may take for granted may surprise newcomers: You can ...


7

Manipulate[ Module[{ eqs = {(t vx0 + x0) Cos[t/2] + t (vy0 + x0/2) Sin[t/2], t (vy0 + x0/2) Cos[t/2] - (t vx0 + x0) Sin[t/2]}}, max = (If[NumberQ[#1], #1, 40] &)[ Quiet@FindRoot[Evaluate[Plus @@ (eqs^2) == 1], {t, 155, 0, 5000}, MaxIterations -> 50][[1, 2]]]; If[max == 0, max = .01]; Show[ParametricPlot3D[Append[eqs, ...


7

Here is a way to generate the mesh including region markers and different refinement in different regions: Needs["NDSolve`FEM`"] L = 1; i1 = 0.625; i2 = 0.25; (bmesh = ToBoundaryMesh[ "Coordinates" -> {{0, -L}, {10, -L}, {10, -L + i1}, {11, -L + i2}, {11, L - i2}, {10, L - i1}, {10, L}, {0, L}, {10, -5}, {12, -5}, {12, 5}, {10, 5}}, ...


7

Here is an approach using RandomPoint and graphics primitives: pts = RandomPoint[Rectangle[], 10^6]; (* generate random points on the unit square *) rm = RegionMember[Disk[{0.5, 0.5}, 0.5]]; (* RegionMemberFunction for an embedded Disk *) Now we count the number of points that fall in the circle and divide that by the total number of points. That should ...


6

init = RandomReal[1, {200, 200}]; alpha = -1; beta = -1; gamma = 1; s = 50; replenish[n_] := RandomReal[1, {n}] step[a_] := MapThread[ Prepend, {MapThread[ Append, {Insert[ Table[a[[i, j]] + alpha (a[[i, j + 1]] - a[[i, j]] + a[[i, j - 1]]) + beta (a[[i + 1, j + 1]] + a[[i - 1, j + 1]] + a[[i + 1, j - 1]] + ...


6

You could use Region functionality (for simpler regions), e.g. rw[pt_, s_, n_, reg_] := Module[{ch = {{0, 0}, {1, 0}, {-1, 0}, {0, 1}, {0, -1}}, np, st}, st = RandomChoice[ch, n]; FoldList[If[RegionMember[reg, #1 + s #2], #1 + s #2, #1] &, pt, st] ] an[p_, step_, num_, regn_] := With[{pnts = rw[p, step, num, regn]}, ListAnimate[ ...


6

This allows you to watch two pendulums swing in time, g = 9.81; pendulumTrajectory[L_, theta0_, t_] := With[ {theta = theta0 Cos[Sqrt[g/L] t]}, {L Sin[theta], -L Cos[theta]} ]; Manipulate[ point1 = pendulumTrajectory[2.01, 15 Degree, t]; point2 = pendulumTrajectory[3.01, 25 Degree, t]; line2 = Line[{{0, 0}, point2}]; Graphics[{{Red, ...


6

Here is a generalization to many pendula. The phase differences in this case are revealed by the special points in time at which "revivals" occur, i.e., where the pendula that have all been launched with the same angle all group back together: pendulum[α_, l_] := Module[{radius = .03, p = -l {Sin[α], Cos[α]}}, {{Brown, Line[{{0, 0}, p}]}, ...


6

This is not answer, but an extremely long comment. I find this problem very interesting, but haven't been able to solve it. In my attempts, I developed a tool to visualize the the two-walker random walk. I am posting this tool because I think it might be useful to the OP or anyone else looking this problem for exploring what's going on. steps = {{0, 1}, ...


5

First, you're not using a fixed step method. (An Euler scheme may be applied to any step size and to one that varies.) To get a true fixed step method you have to turn off "DiscontinuityProcessing" when you have a discontinuous ODE; otherwise, NDSolve will try to adapt the steps to account for the discontinuity. The "DiscontinuityProcessing" stage resets ...


5

You go to the casino. You bet a dollar the ball will land on red. If it does you get two dollars if it does not you get zero. Likewise if you bet the ball will land on black. You cannot bet it will land on green. Use this previous answer to help do the pattern matching. Note to the original poster, I found that using Google and searching for mathematica ...


5

You can make a function for the trajectory of a projectile, which takes the initial position, velocity, angle, and time as arguments. projTrajectory[{x0_, y0_}, v0_, theta_, t_] := Module[{T = -(v0 Sin[theta] + Sqrt[-2 g y0 + v0^2 Sin[theta]^2])/g}, {x0, y0} + If[t < T, {v0*Cos[theta]*t, v0*Sin[theta]*t + 0.5 g*t^2}, {v0*Cos[theta]*T, ...


4

Here's a q-d implementation: With[{hits = Cases[Partition[RandomChoice[{17/36, 17/36, 1/18} -> {1, 2, 3}, 1000000], 5, 1], {x_, x_, x_, x_, y_}]}, Count[hits, {x_, x_, x_, x_, y_} /; y != x] - Length@hits] It will run through 1 million spins, checking for any 4-runs and counting those where the subsequent spin differs, less the count in total (so ...


4

I misunderstood and commented only about computing the Cesaro means as per the question In s I have all partial sums, but I do not know how to divide them by corresponding n. The desired scatter plot of 500 Monte Carlo attempts with samples of increasing length could be obtained with something like ticks = Range[-0.06, 0.06, 0.02]; s = Table[u = ...


4

In your example plot, the Monte Carlo integral is computed afresh for each new amount of sampling points: s[n_] := Total[Sqrt[1 - #^2] & /@ RandomReal[1, n]]/n Then take one random point, find the mean, take two new random points, find their mean, and so on until 500. The result is: ListPlot[Table[s[n], {n, 500}], PlotRange -> {0.7, 0.9}] The ...


4

dist = TransformedDistribution[x + y, { Distributed[x, DiscreteUniformDistribution[{1, 6}]], Distributed[y, DiscreteUniformDistribution[{1, 6}]]}]; SeedRandom[1] For small sample sizes, the match to the theoretical values is poor. data = Total /@ RandomInteger[{1, 6}, {100, 2}]; Show[ Histogram[data, {1.5, 12.5, 1}, "PDF", ...


4

You can use RandomVariate to sample from a DiscreteUniformDistribution and then add up the pairs, calculate the probability of the sums observed, and then extract the probabilities of interest. (#/100. & /@ Counts[Plus @@@ RandomVariate[ DiscreteUniformDistribution[{1, 6}], {100, 2}]] )[#] & /@ Range[2, 6] Hope this helps.


4

This is not a direct answer to your problem, but rather a generalization of @Jens' code from double to n-tuple pendulums. Meaning you can also use it for double pendulums if you like. I'm providing it due to popular demand. Needs["VariationalMethods`"] n = 10; (* number of pendulum segments *) rate = 10; (* animation frame rate *) Clear[s, ϕ, t, g, m]; ...



Only top voted, non community-wiki answers of a minimum length are eligible