# Tag Info

117

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

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

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

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

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

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

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

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

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

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

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

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

Could work with a low discrepancy rather than pseudorandom sequence. These tend to be better for avoiding approximation errors associated with "clumping" (and the corresponding "empty regions") that one gets with the latter. The code below will provide a sequence that is uniformly distributed modulo 1 in the unit square and seems to be low discrepancy, ...

4

EDITED Alex, I believe this is what you want: StartingWealth = 100; PercentChange = 0.5; WinProbability = .5; NumberOfProcesses = 2; Time = 5; processes = RandomFunction[BinomialProcess[WinProbability], {0, Time}, NumberOfProcesses]; paths = Table[ FoldList[Times, StartingWealth, If[Differences[Last[Transpose[processes["Path", x]]]][[#]] == ...

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

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

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

3

Your simulation's convergence to an incorrect result is because your comparison accepts any negative y, thus inflating the result. You are in essence computing the area of the union of the semi-unit circle that lies above the x-axis and the rectangle with top-left corner at {-1,0} and lower-right corner at {1,-1}, which is 2 + π/2 = 3.5708. This is the ...

3

Here is another way to solve this issue. I know, it is written everywhere, but this is a common mistake to think that the orifice equation is Cd * Ad * Sqrt[ 2/Rho * ( ps - p1[t] ) ] and could give complex results ! There is no physics going imaginary in the real world. The equation is just wrong. When the pressure reverse, the flow reverse, at least. ...

3

You might use rule replacement (instead of variable assignment) to inject specific values for the parameters. ee = C1 - 10 + PG/2 + GC1 + GL1 == 0 && C2 - 15 + PG/2 + GC2 + GL2 == 0 && GL1 - fd (a1) t 10 == 0 && GL2 - fd (a2) t 15 == 0 && GC1 - (t (1 - fd) 25 - PG) ((p1)/((p1) + (p2))) == 0 && GC2 - (t ...

3

To correctly compute the mean, try this: Manipulate[SeedRandom[seed]; meanvector := Mean[assets]; assets = Table[RandomFunction[GeometricBrownianMotionProcess[μ, σ, S0], {0, time, 0.1}]["Path"], {P}]; G1 := ListLogPlot[assets, GridLines -> {{}, {watermark}}, GridLinesStyle -> Directive[Green, Thick], Joined -> True, AxesLabel -> {"Time", "St"}, ...

3

Like this? NMaximize[{x[1] + x[1]^2 + #[[1]] x[2] + 2 Sin[Sin[x[3]]] + x[4] + Sin[Sin[x[5]]] + Sin[x[1]] + x[2] + Sin[x[3]] + x[4] + #[[2]] x[5], x[1] > 0 && x[2] > 3 && x[3] > 1 && x[4] > 0 && x[5] > 0 && x[1] <= 1 && x[2] <= 8 && x[3] <= 1 && ...

3

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

2

There is an open-source Modelica system called OpenModelica at www.openmodelica.org. Frankly, I am a little baffled at what you are calling a "downside". Modelica is an object oriented language, so only different classes need to be defined within the library. A single class can be implemented as many objects within your model. For example, a resistor ...

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