# Tag Info

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

12

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

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

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

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

You could have the velocity at which it spins decay randomly in time. For example, here I have the angle at time $t$ satisfy $\theta''(t)=-f(t)\theta'(t)$ with $f(t)$ a random(ish) function (see later). If $\theta(0)=1$ and $\theta'(0)=0$ then $\theta'(t)\rightarrow 0$ for large $t$. So now we just need the $f(t)$, which I construct by interpolating a ...

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

5

I suggest a different user interface, for fun. The user clicks down the mouse and drags. An arrow will be drawn that indicates the magnitude and direction of a "flick" of the spinner. When the mouse is released ("MouseUp"), a random destination is computed and the index of the spinner spins to it, decelerating as if under constant deceleration. The ...

5

Here's a way to animate the pointer: Manipulate[ang = 0.9*ang + dest; Column[{Show[PieChart[{1, 1, 1, 1, 1}, ChartStyle -> "DarkRainbow", ChartLabels -> Placed[{2, 4, 6, 8, 10}, "RadialCenter", Style[#, 24] &], LabelingFunction -> None], Graphics[{{PointSize[0.05], Point[{0, 0}]}, {Thick, Arrow[{{0, 0}, 0.7 ...

4

You can do something like the following: First define the DE in every node node[1] = β*M*L^2*D[Θ[i][t],{t,2}]==M*g*μ*L*(-Sin[Θ[i][t]])+τ[Θ[i+1][t],Θ[i][t]]; node[n] = β*M*L^2*D[Θ[i][t],{t,2}]==M*g*μ*L*(-Sin[Θ[i][t]])+τ[Θ[i-1][t],Θ[i][t]]; node[i_Integer/;1$<$i$<$n] = ...

4

Here is a faster way. I've only coded a portion of your problem to provide the flavor. Roll your own GBM GBMPathCompiled = Compile[{{S0, _Real}, {drift, _Real}, {diff, _Real}, {nSteps, _Integer}}, FoldList[(#1 drift Exp[ diff #2]) &, S0, RandomVariate[NormalDistribution[0, 1], nSteps]]]; Let's first visualize the situation Module[{S0 = ...

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

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

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

Too long for comment: Look at the simulation demonstrations and pick one which makes you think "That's the stuff I'd like to do!" For instance predator-pray ecosystem there is a link to Download Author Code on each demonstration page, do that and open it in Mathematica and if it asks you to enable dynamics, do it. Now you have something like this: ...

3

Here are some tweaks, plus a way to generate non-adjacent blocked sites. (I wasn't sure whether you wanted two free sites between blocked sites or just one. There's just one below, and comment where it could be changed to two.) The run time (per hop) is shorter by about 2/3, only a modest improvement. kBack1[k0_] := k0 Exp[-2]; kFor1[k0_] := k0 Exp[-2]; ...

3

Regarding your first point: The planetary orbital planes are indeed inclined to that of the Earth and you have the inclinations. However, you also need to know the azimuthal locations of the ascending and descending nodes which together define the line about which to pivot the orbital ellipse. Regarding your second point: If you assume, as you say, that ...

3

This might not be the canonical way, but the way I would do it is to represent the system as a pair of state-updating equations, translate that into matrix form and use NestList to show the time path. This is of course assuming that you are happy to work in discrete time. Consider, for example, where the reproduction rate at time $t$ is a positive function ...

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

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

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

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

2

Of course the true exponential function is the solution of a differential equation in which the time steps are taken to be infinitesimal. However, it seems like you're looking for an implementation that's as literal as possible and uses discrete time steps: coli[rate_, nStart_, n_] := NestList[Function[{coli}, Round[coli + rate coli]], nStart, n] ...

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