# How to run multiple random process and plot its average for a non-autonomous logistic model?

This problem is a continuation of the discussion from here,
with random impulse function as discussed in here to formed the following:

n = 200;
t1 = RandomVariate[NormalDistribution[5, 2], n];
t2 = 0.7*RandomReal[{0, 1}, n];
t1 = Accumulate[t1];
dat = Transpose[{t1, t2}];
train[t_] := Interpolation[dat][t];

Do[SeedRandom[1234];
at = train[t];
dath = {x'[t] == 5*x[t]*(1 - x[t]) - at*(1 + Sin[2*\[Pi]*t])};
psect[x0_] := Reap[NDSolve[{dath, x[0] == x0,
WhenEvent[Mod[t, 0.01] == 0, Sow[{t, x[t]}]]}, x, {t, 100}]][[2, 1]];
adata = Map[psect, {0.92,0.5}]; ,10];

ListPlot[adata, ImageSize -> Medium]


which produced

Do is introduced to allow the function repeats the process multiple times.
The average of this procedure is plotted to map the final result.

avg = Mean[Array[psect, 10]];
ListPlot[avg]


However, I'm not sure if this is how it should be.
Because I'm expecting Do to repeat the process (in this case 10 times)
and produce a final (perhaps single?) graph, as defined by function Mean.
This is not illustrated in the second picture, and it seems to be plot noise (maybe?).
But again, I'm not entirely sure because I'd planning to run the process potentially 100 or maybe 1000 times before taking the average. And for that, I'd like to ask for some guidance.

• You need to say what you mean by "messy", what a desirable outcome would be, and what problem you would like to fix. Asking for "guidance" is far too generic and broad. Commented Apr 8, 2022 at 14:13

## 1 Answer

There's a lot going on here. The main problem is that your replicate solutions are all forced by the same random function train[t] that you define before the loop. So right now, the only difference between your replicates is in the initial conditions, which quickly disappears as they all reach the same attractor (see the blue vs yellow line in your first figure near $$t=0$$). I think you want to re-initialize the random function before each replicate solution.

Let me make some other suggestions. First, since most of the time the system is on the attractor, which should be ergodic, you can focus on the statistics from a single longer run. Second, the random forcing is pretty ad hoc, so it isn't easy to connect with mathematical theory or others without Mathematica to replicate. How about a simpler random forcing, consisting of a random step function? Third, all the Sow and Reap seems extraneous - just solve it and make a table later.

As an example,

tmax = 100;
SeedRandom[1234];

(* random harvesting function *)
h = Interpolation[
Transpose[{Range[0, tmax], RandomReal[{0, 1}, tmax + 1]}],
InterpolationOrder -> 0];

(* solve DE *)
sol = NDSolve[{x'[t] == 5*x[t]*(1 - x[t]) - h[t], x[0] == 1}, x, {t, 0, tmax}][[1]];

Plot[{h[t], x[t] /. sol}, {t, 0, tmax}]

Histogram[
Table[x[t] /. sol, {t, 0, tmax, 0.01}], Automatic, "Probability",
PlotRange -> {{0, All}, All}]

Mean[Table[x[t] /. sol, {t, 0, tmax, 0.01}]]


(* 0.87401 *)

Addendum:

If you really want replicates, try something like this. I like using Table instead of Do since it the output is built-in.

tmax = 100;
nrep = 10;
SeedRandom[1234];
dat = Table[
h = Interpolation[Transpose[{Range[0, tmax], RandomReal[{0, 1}, tmax + 1]}], InterpolationOrder -> 0];
sol = NDSolve[{x'[t] == 5*x[t]*(1 - x[t]) - h[t], x[0] == 1}, x, {t, 0, tmax}][[1]];
Table[x[t] /. sol, {t, 0, tmax, 0.01}]
, nrep];

Show[
ListLinePlot[dat, DataRange -> {0, 100}],
ListLinePlot[Mean[dat], DataRange -> {0, 100}, PlotStyle -> {Black, Thick}]
]


If you increase the number of replicates, you'll see the average over replicates approach the average over time.

• I've tried replicating your code with train[t] as I'd like to repeatedly(and randomly) disturb the continuity of the system. But, not able to plot the train[t]. Yes, you're right @ChrisK. I want to re-initialize the random function before replicating each solution, then take the total average to plot the final. But, I'm not sure how to do the latter, which I tried with Do and potentially is not the best way to do it? Commented Apr 9, 2022 at 1:46
• I've replaced train[t] with the simpler h[t], which does randomly change over time. As for replicates, there's not much point because once the system is on its attractor, averaging over time should be the same as averaging over replicates. See en.wikipedia.org/wiki/Ergodicity Commented Apr 9, 2022 at 4:14
• I've updated my answer as a demonstration of how to get replicates. Commented Apr 9, 2022 at 4:32
• Agreed, h[t] is better to represent train[t]. After hours of trying to figure this out on a more complex system, thanks @ChrisK! Commented Apr 9, 2022 at 14:27
• NDSolve returns its answer as a list of solutions. Unless the system does not have unique solutions (a bad sign for time-dynamic systems in my opinion), there will only be one. [[1]] simply peels off that extra set of {} to extract the solution. Maybe not always necessary, but I don't want to end up with an extra level of {} later on, so I just do it. Commented Apr 11, 2022 at 15:21