Hot answers tagged simulation
17
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 &=& ...
15
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 = ...
14
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 @ ...
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 ...
10
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 = ...
7
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 -> ...
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
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
...
7
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 ...
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
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 ...
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
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
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 ...
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] = ...
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
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 = ...
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
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 ...
3
This is an alternative approach:
boundary = 5;
step = .1;
AbsoluteTiming[
rwtbl = NestList[
Mod[
# + RandomVariate[NormalDistribution[0, step], 3],
boundary] &,
{0, 0, 0},
10000
];
]
Graphics3D@Line[
rwtbl
]
(it's also 2400 times faster than your code on my machine).
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];
...
2
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:
...
2
I have code that works about 50x faster now. The trick was to use the GraphUtilities package function NeighborhoodVertices[] rather than VertexList[NeighborhoodGraph[]]. It sped up these two lines from about 0.5sec each to 0.007sec. Since these are called 1000 times or more, that's a big help.
I am still interested in a more functional way of doing this. ...
2
I didn't manage to make it work with RandomFunction, but if you would consider running simulations without it you could use TruncatedDistribution.
Define the parameters :
mu = 5/100; sigma = 15/100; x0 = 100;
xmin = 80; xmax = 200;
tmin = 1/100; tmax = 10; dt = 1/100;
Define the truncated distribution :
trDist[xmin_, xmax_, mu_, sigma_, x0_, t_] := ...
2
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 ...
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]
...
1
Here is something that works with the definitions you gave and does what I understand you try to do:
n = 10;
fullSystem = Join[
{
β*M*L^2*D[Θ[i][t],{t,2}] == M*g*μ*L*(-Sin[Θ[i][t]]) + τ[Θ[i+1][t],Θ[i][t]]
} /. i -> 0,
Table[
β*M*L^2*D[Θ[i][t],{t,2}]
==
M*g*μ*L*(-Sin[Θ[i][t]]) + τ[Θ[i-1][t],Θ[i][t]] + τ[Θ[i+1][t],Θ[i][t]],
...
1
Then, you must build the logic to handle the partial barriers.
There are lots of ways to proceed. I modified the visualization to illustrate how indexing the sample times can be used to restrict the barrier penetration test and select those paths that survive the first barrier.
Clear[PartialBarrierTest];
PartialBarrierTest[path_, iSt_, iFn_, U_, L_] :=
...
Only top voted, non community-wiki answers of a minimum length are eligible
