# Tag Info

116

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

42

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

20

You are trying to implement Euler-Maruyama simulation method for a 2-stage short-term interest rate model which is given by the following system of SDEs: \begin{eqnarray} \mathrm{d} \theta_t &=& -\lambda_\theta \left( \theta_t - \bar\theta\right) \mathrm{d}t + \sigma_\theta \mathrm{d}W_{\theta,t} \\ \mathrm{d} \pi_t &=& ...

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

19

I'm not sure what your goal is, exactly, but here is a simulation I cooked up. It should give you some ideas: metersToAU[m_] := m/(1.496*10^11) ; orbit = First@AstronomicalData["Earth", "OrbitPath"]; earthCurrentPosition = AstronomicalData["Earth", "Position"] // metersToAU; radiusEarth = AstronomicalData["Earth", "Radius"] // metersToAU; radiusSun = ...

17

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

16

This answer is going to be a bit of a sprawl. Please read on. I am going to present several methods of simulation, hopefully in increasing order of performance. Method 1 We can carry out the filling of seats, at least as I understand the puzzle, quite literally like this: fillseats[seats_List] := ReplacePart[seats, {{1}, {2}} + RandomChoice @ ...

14

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]) & /@ ...

14

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

13

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}}]; ...

12

For the example you give there is no reason you can't use NestList, you just need to make two simple changes: Don't use the side effect in deltaπt to get the value for θnow, give it as an explicit second argument Then you just do: NestList[{deltaπt[#[[1]], #[[2]]] + #[[1]], deltaθt[#[[2]]] + #[[2]]} &, {2, 2}, noYear] Similar changes would allow ...

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

I would proceed like the following. It will be natural to propose that the win-event occurs following BinomialDistribution with probability $p=0.8$ so that we can use the built-in BinomialProcess to simulate the win and losses in $20$ time steps for $50$ sample paths. timstep = 20; win = BinomialProcess[.8]; samplepaths=50; process = ...

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

10

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

8

Just the code:) Animate[Show[ Graphics[Translate[Rotate[{Circle[], Thick, Blue, Line[{{0, 0}, {0, -1}}], Red, PointSize[.02], Point[{0, -1}]}, -t], {t , 0}], PlotRange -> {{0, 4 Pi}, {-2, 2}}, ImageSize -> {Large, Tiny}, Axes -> {True, False}, AxesOrigin -> {0, -1}], ParametricPlot[{(a - Sin[a]), (-Cos[a])}, {a, 0, t}, PlotStyle -> ...

8

I managed to double the speed by simply moving some multiplications and divisions... AbsoluteTiming[ or2 = -2./Omegar^2.; oz2 = 2./Omegaz^2.; p = p^2.; Io*Table[ Exp[(p[[i, t, 1]] + p[[i, t, 2]])*or2 - p[[i, t, 3]]*oz2], {t, 1, Deltat + 1}, {i, n}]; ] A dot product is 30 to 40% faster still. AbsoluteTiming[ or2 = ...

7

There are two areas for optimization that I see here. The first, if possible, is to generate all your random data in advance and then access it with an incrementing index, e.g. list[[i++]]. The second is to partially evaluate the definitions of thetaNext and piNext for a given set of parameters. A note: Random has been deprecated for some time now and may ...

7

Rapid calculations are afforded by Accumulate to generate the walk and Mod to implement the periodicity. Scaling the entire thing to the unit cube simplifies the code a little. Furthermore, don't generate normally distributed displacements: uniform displacements will do when the increments are small. With these efficiencies we may generate the coordinates ...

7

I'll switch it up a bit: I'll give you somewhat simplified code, and your task is to figure out what I'm trying to do: With[{frames = 15}, Animate[ ParametricPlot[{u - Sin[u], 1 - Cos[u]}, {u, -\$MachineEpsilon, t}, Axes -> None, Epilog -> {Line[{{t, 1}, {t - Sin[t], 1 - Cos[t]}}], ...

7

As I mentioned in the comments, using Mod[] is one good way to enforce your periodic boundary conditions: just generate the random walk as usual, and then apply Mod[] to bring back inside the sections that are outside your box. Here's an example of what I'm describing: n = 5*10^3; (* number of steps *) s = 20; (* cube edge length *) h = 1/10; (* step bound ...

7

The problem seems to be that outArea = Last /@ ComponentMeasurements[points, "Area"] // Total; estimates the area of the whole square image, not of the disk or the points. For instance, with SeedRandom[1]; points = Show[ Graphics[{Pink, Point /@ Select[Partition[RandomReal[{-1, 1}, 100000], 2], ({a, b} = #; a^2 + b^2 <= 0.98) ...

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

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

6

The limitation you quote is not a general limitation of Modelica. It is possible to define a Modelica component that has a variable number of inputs/outputs. Typically the number of inputs/outputs is then given by a parameter to that component. For example, the following component has one input but 2 outputs, varied with the parameter nout: model SIMO ...

6

You can get the state values for every with data["States"], which you can then easily feed into a indicator function. data = RandomFunction[ GeometricBrownianMotionProcess[0.01, .15, 100], {0, 1, .01}, 100]; corridorIndicator[data_, upperBound_, lowerBound_] := Boole[Max@# < upperBound && Min@# > lowerBound] & /@ data ...

6

There's no need to iterate over i and t to build particleintensity, just do the whole dot product in one go. The Exp moves outside the dot product (so we use it once as a Listable function rather than 5 million times individually). particleintensity = Io Exp[(p^2).{-2/Omegar^2, -2/Omegar^2, -2/Omegaz^2}]; int = Total[particleintensity]; This runs ...

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