# Tag Info

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

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

13

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

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

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

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

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

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