# Tag Info

## Hot answers tagged simulation

19

This was a fun question to answer, even considering that I know nothing about general relativity. It's all a matter of translating the equations presented in this paper by Oliver James, Eugenie von Tunzelmann, Paul Franklin, and Kip Thorne into notebook expressions. Embedding Diagrams The paper gives some really cool figures to show the curvature of 4-...

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 ]; ListAnimate[%...

16

Yes you can. Below is a fairly general, Mathematica-compiled, fast and robust version. Examples 1. Michaelis-Menten kinetics Michaelis-Menten kinetics for enzyme-directed substrate conversion. The enzyme (e) converts the susbtrate (s) through an enzyme-substrate complex (c) to the product (p). For comparison, I've included the deterministic ODE system ...

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 = ConstantArray[0,...

12

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 - 1]/(...

12

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] =...

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 oneRound[...

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 -&...

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

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 ...

8

Here is a way to generate the mesh including region markers and different refinement in different regions: Needs["NDSolveFEM"] 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}}, ...

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[]}, ...

8

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 ...

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

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[ Graphics[...

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

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, 0....

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 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

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

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}, {1,...

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, ...

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 ...

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

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

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 ...

3

Here are solutions to both boundary protocols. They are built on the same basic framework -- mainly the function that generates the moves for the walker is what differs between the two. There is a little adjustment in the way the lines and walker point is drawn because of discontinuities in the path generated by the wrap-arround protocol, Path clips at the ...

3

To help get you started I attempted to translate equation (5) from your snapshot into Mathematica. Step 1 I(H) from the snapshot becomes: i[h_] := Piecewise[{ {Gamma[1 - 2 h]/h Sin[π/2 (1 - 2 h)], 0 < h < 1/2}, {Gamma[2 (1 - h)]/(h (2 h - 1)) Sin[π/2 (2 h - 1)], 1/2 < h < 1}, {π, h == 1/2} }] This is a function that can be ...

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