# Tag Info

## Hot answers tagged simulation

125

I did a very simple (in fact over-simple) snowflake simulator with CellularAutomaton years before. It's based on the hexagonal grid: and range-1 rules: Initial code First we'll need some functions to display our snowflakes: Clear[vertexFunc] vertexFunc = Compile[{{para, _Real, 1}}, Module[{center, ratio}, center = para[[1 ;; 2]]; ratio = ...

45

========== update =========== Remember guys how we can cut out a snowflake from a sheet of paper carving 12th folded part? Like the image below. So I decided to write an app to imitate the process. It also can be used to make random snowflakes (similar to to @bill s' but with reflection to imitate real cutting paper process and reflective symmetry of ...

29

The only reason I am attempting to answer this is to perhaps get a Reversal badge. There you go... We will go slowly and this answer is the basis for what comes next. Let's start with two dimensions. You'll see why. We create a rectangular region: Needs["NDSolveFEM"] mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, ...

20

The programming style you are using is not very fitting for Mathematica. Here's a better way (shorter, much faster): n = 1000000; (* number of points to use *) octantVolume = N[ Total@UnitStep[1 - Norm /@ RandomReal[1, {n, 3}]]/n ] The reason why you get the error you mention is that for some x, y, the expression 1 - x^2 - y^2 is negative, thus its ...

20

Here is a simple method that begins with an $n$-sided polygon (defined by the $n$ points in tab), then rotates the polygon and superimposes it six times to achieve the six-fold symmetry. The makeFlake function is: makeFlake[n_] := Module[{tab, rot}, tab = RandomReal[{-1/2, 1/2}, {n, 2}]; rot = RotationMatrix[Pi/3]; Graphics[{Hue[RandomReal[]], ...

16

Not so much snowflakes as random artworks with the same symmetry as snowflakes, but I wanted to join in the festive fun! These are generated with a "randomart" package I wrote a while ago (code at the bottom of the answer). It uses a kind of non-linear iterated function system to generate random images. Here's a grid of random images with snowflake ...

15

Method of random number generation is also significant: Default: n = 10^6; AbsoluteTiming[N@Mean@UnitStep[1. - Total[RandomReal[1, {3, n}]^2]] - π/6] {0.197896, 0.000649224} Niederreiter low-discrepancy sequence (see "methods" here): SeedRandom[Method -> {"MKL", Method -> {"Niederreiter", "Dimension" -> 3}}]; ...

15

When n is large it's much faster to operate on a 3 x n array than to process each of the n 3-vectors separately. This is one of the standard "tricks" to speed things up. n = 10^6; (* Isn't that easier to read than 1000000 ? *) AbsoluteTiming @ N[ Total@UnitStep[ 1. - Norm/@RandomReal[1,{n,3}] ]/n ] (* {4.555842, 0.524302} *) AbsoluteTiming @ N[ ...

14

Here is an un-golfed and simplified version of an L-System production based on a previous answer of mine: f1[initState_, rotAngle_, prodRules_, iters_] := Module[{currAngle = 0, currPos = {0, 0}, res = {}}, (res = {res, Line@{currPos, currPos += {Cos@currAngle, Sin@currAngle}}}; If[NumericQ@#, currAngle += I^# rotAngle]) & /@ ...

13

The separation-of-variables solution you quoted has two indices appearing in it: n and j (the subscripts of the coefficient $A_{nj}$). Here, n is azimuthal mode order, i.e. it counts the number of nodes along the direction in which the polar-angle $\theta$ varies (divided by 2). The index j is needed because the wave is supposed to satisfy the boundary ...

13

The printed version of the 2002 edition was printed 3 times and sold out 3 times; Springer and Google recently started selling it (book only) as a PDF eBook (no software) on the Springer and Google sites for \$79. I know other authors (e.g. here) have gone to some trouble to make their books available here on stack exchange ... We are delighted to be able ...

12

One way is to set up a DAE: See tutorial/DSolveExamplesOfDAEs and example/ModelConstrainedSystemsAsDAEs. The constraint that the driver (bottom rotating link) has a fixed length is taken care of by initial conditions and the DE. There are two possible starting positions for the driven link. One might have to inspect the result of Solve to determine which ...

12

Well I guess one more couldn't hurt. Using an iterated matrix-replacement scheme and some fancy opacity: powzerz = 2; width = 550; primitive = Scale[Cuboid[], 0.99999]; matrix0 = {{{1}}}; matrixT = CrossMatrix[{1, 1, 1}]; rules = {0 -> (0 #1 &), 1 -> (#1 &)}; iterate[matrix0_, matrixT_, rules_, power_] := Nest[Function[prev, ...

11

This is simplest implementation. If a new crater gets closer than 30 to some old craters, only closest old crater is getting replaced with new one. You can built on this example something more sophisticated. craters = {{0, 0}}; number = {1}; Dynamic[new = RandomReal[{-250, 250}, 2]; near = Nearest[craters, new][[1]]; Row[{ Graphics[{PointSize[.05], ...

11

A smooth changing fractal snowflake: {s, d, t} = {0, 1, 3}; Dynamic@Graphics@ Polygon@Reap[ If[# != 0, t += 8.^-5; Do[#0[# - 1]; Sow[d = Sign@d #; {Re[s += d], Im@s}] & /@ (# E^(I t #) &@ Range@6/(5^(4 - #))); d *= E^((\[Pi] - 63 t)/3 I), {6}]] &@ 3][[2, 1]]

11

"But, I just need how the crank slides along the rod. The rest, I can try it as an exercise" I'll help you out with this detail to get you started. It's really just a case of adding vectors together: The disk moves along a circle inside a circle, so if the outer circle has radius 1 the equation for its movement might be 0.8 {Cos[theta], Sin[theta]}. ...

9

I did a solution with contour tracing on the distance function. It gets pretty unstable sometimes, but it's a fun question to experiment with interactivity. DynamicModule[{p1 = {0, 2}, p2 = {1, 3}, angles = {0, 0}, distance, grad, tangent}, distance[a1_, a2_] := Norm[{Cos@a1, Sin@a1} - (Norm[p2 - p1] {Cos@a2, Sin@a2} + p1)]; grad = ...

8

A couple things: You don't update velocity (obx and oby appear to be unused?). You do your timestep wrong (you want x+=velocity*timestep, not x=velocity*time!). You also have your direction of your force wrong. You should take the negative gradient of potential. The gradient of 1/r is -1/r^2. The negative of that is 1/r^2. So you get the common sense result ...

7

Very similar to Vitaliy's answer, but deleting all craters within the critical distance, and somewhat more compact: craters = {{0, 0}}; number = {1}; Dynamic[ (craters = #; Row[{Graphics[{PointSize@.05, Point@#}, ImageSize -> 230, PlotRange -> 300, Frame -> True], ListLinePlot[AppendTo[number, Length@#], PlotRange -> All, ImageSize ...

7

With eqn[{k_, r_, H0_, P0_}] := {H'[t] == r (1 - H[t]/k) - d H[t] P[t], P'[t] == -s P[t] + e H[t] P[t], H[0] == H0, P[0] == P0} d = 0.01; s = 0.3; e = 0.02; I would define one simulation as sim := Module[ {k = RandomVariate[NormalDistribution[150, 20]], r = RandomVariate[NormalDistribution[0.4, 0.003]], H0 = RandomVariate[UniformDistribution[{50, ...

6

You can do this without a NDSolve by calculating the distance from the follower cranks joint to the end of the driving cranks end. Then use this distance with law of cosines to calculate the deviation angle. This is also pretty easy to implement on ANY hardware capable of doing a ArcCos and Atan2 operation (note that in c atan2 parameters are swapped). ...

5

With f[a_, {x_, y_}] := Piecewise[{{a x, x < y}, {x, x == y}, {1 - a + a x, x > y}}, 0] one simulation can be defined by sim[length_] := Module[{rv = RandomVariate[BetaDistribution[3, 1], length], y, yBar}, y[1] = First@rv; yBar[t_Integer] := yBar[t] = 1/t * Sum[y[i], {i, 1, t}]; y[t_Integer] := y[t] = f[0.5, {rv[[t]], yBar[t - 1]}]; ...

4

You could try specifying the event differently (this will trigger when the absolute difference between p1 and p2 is smaller than some value): Ad1=10/1000^2; Ad2=1.5*1000^(-2); Ad3=1.5*1000^(-2); Cd1=0.67; Cd2=0.67; Cd3=0.67; V1=10/1000; V2=10/1000; Rho=875; beta=1000*10^6; ps=100*10^5; Q1=Ad1*Cd1*Sqrt[(2/Rho)*(ps-p1[t])]; ...

4

The problem with WhenEvent has to do with the OP's DE. For an event to be detected, there has to be a point at which the condition is crossed, that is, changes from False for t < t0 to True for t > t0. NDSolve then applies a root-finding algorithm to approximate the value of t0 at which the event occurs. In your DE, the solution p1[t] theoretically ...

4

You can display zero crossing using MeshFunctions. Here is a clumsy exploitation from created graphic. The half-periods (difference between consecutive points) are displayed below with mean in red. x1plot = ListPlot[x1data, AxesLabel -> {"t", "x1"}, Joined -> True, MeshFunctions -> (#2 &), Mesh -> {{0.}}, MeshStyle -> {Red, ...

4

With 1. notes in the answer of your latter question, 2. initial values for FQout and RQout(the values are carelessly chosen real numbers because it's bound to be overlaid), the modified code is about 40 times faster now: f = With[{dx = 1/6, n = 48, m = 300, p = 36, capacity = 2500, A = 18., B = 0.1, L = 3., RML = 30, Vf = 100, Kj = 150, w = 20, ...

4

I think your core idea about post-processing the splittings is a good one. It can be expressed somewhat more elegantly (although elegance is a subjective matter) with the help of linked lists rather than Reap/Sow: ClearAll[ll]; SetAttributes[ll,HoldAllComplete]; toLinkedList[l_List]:= Fold[ll[#2,#1]&,ll[],Reverse[l]] and ClearAll[process]; ...

4

It looks like I'm late to the party but here is my approach. Common functions and data: trinarize[a_List, γ_?NumericQ] := UnitStep[a - γ] + UnitStep[γ + a] labels = {0 -> "bottom", 1 -> "middle", 2 -> "top"}; SeedRandom[1] sampleSimulation = RandomChoice[{-1, 0, 1}, 20] {0, -1, 0, 0, -1, -1, -1, 0, -1, -1, -1, -1, 1, -1, 0, 1, -1, -1, 0, 0} ...

4

I do not know why, but the following are much faster than Walk: ClearAll[walkF,walkF2]; walkF[x_List, a_List, b_List, sigma_] := With[{f = TruncatedDistribution[Transpose[{a - x, b - x}], ProductDistribution @@ ConstantArray[NormalDistribution[0, sigma], Length[x]]]}, x + RandomVariate[f]] walkF2[x_List, a_List, b_List, sigma_] := ...

3

2 problems: You were comparing, in the WhenEvent, solution, which had complex value at that t, to real numbers. I used Abs. If this does not work for you, you can use Re, but can't compare complex number to real number using >. Second, your system is stiff, need to use StiffnessSwitching to help NDSolve. d1 = 10/1000^2; Ad1 = 10/1000^2; Ad2 = ...

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