1
$\begingroup$

I'm performing a stochastic evaluation, where i'm interested in the assymptotic behavior of the solutions, but my computer can't stand very large times. So I thought that I could evaluate a certain time $\Delta t$, save the last state of the system and clear the previous ones (in order to save memory) and keep this process $n$ times, but I don't have any ideia about how to implement this in my code. My code is the following,

a = .3;
μ = .1;
c = .2;
σ = 0.1;
n = 2500;
sol2 = RandomFunction[
   ItoProcess[{\[DifferentialD]s[t] == -a s[t] i[
        t] \[DifferentialD]t, \[DifferentialD]i[
        t] == (a s[t] i[t] - μ i[t] + 
          c (1 - s[t] - i[t]) i[t]) \[DifferentialD]t + σ i[
         t] \[DifferentialD]W[t]}, {s[t], i[t]}, {{s, i}, {0.3, 0.7}},
     t, W \[Distributed] WienerProcess[0, 1]], {0, 50, 0.01}, n];

I hope the problem is clear, any doubt about leave the comments. Thanks in advance!

$\endgroup$
1
$\begingroup$

It is indeed a bit unfortunate that the starting conditions have to be specified in ItoProcess and not in RandomFunction. So, to restart a simulation from a previously obtained dataset would require set up one new ItoProcess for each entry.

However, it is actually not necessary and (contrary to RandomVariate) also not benefitial to generate many samples at once: In the following you will find two way of generating 600 instance of random functions

a = .3;
μ = .1;
c = .2;
σ = 0.1;
proc = ItoProcess[{\[DifferentialD]s[t] == -a s[t] i[
       t] \[DifferentialD]t, \[DifferentialD]i[
       t] == (a s[t] i[t] - μ i[t] + 
         c (1 - s[t] - i[t]) i[t]) \[DifferentialD]t + σ i[
        t] \[DifferentialD]W[t]}, {s[t], i[t]}, {{s, i}, {0.3, 0.7}}, 
   t, W \[Distributed] WienerProcess[0, 1]];

aa = ParallelTable[RandomFunction[proc, {0, 50, 0.01}, 30], {20}]; // AbsoluteTiming // First
bb = RandomFunction[proc, {0, 50, 0.01}, 600]; // AbsoluteTiming // First

4.08794

14.5203

(I ran this on a Quad-Core CPU.) So, RandomFunction on its own seems to be not parallelized. Hence it is a two-fold good idea to compute smaller chunks in parallel over the whole time horizon and to save the results to file. This way, computations are i) faster and ii) less RAM is required.]]

Saving to file can be done like this:

i = 1;
Export["file" <> ToString[i] <> ".mx", aa];

Here i is a running index that should be incremented with each generated dataset.

Data can be retrieved later as follows:

aa = Import["file" <> ToString[i] <> ".mx"];

Word of warning

The mx-file format is very fast to save and load, but it is not appropriate for archiving: These files are OS- and also version-dependent.

$\endgroup$
  • $\begingroup$ Are they really OS-dependent? I work on Windows 10 64-bit, and I share my MX files freely with colleagues who run other operating systems. MX files are certainly version-dependent, but only in the backwards direction. A more recent version can always load MX files from an older version, but not the other way around. $\endgroup$ – Shredderroy Sep 3 '18 at 22:06
  • $\begingroup$ I actually haven't tried it. The MX documentation says: "MX files cannot be exchanged between operating systems that differ in $SystemWordLength." Well, there aren't that many computers that still run 32 bit OSes, but... $\endgroup$ – Henrik Schumacher Sep 3 '18 at 22:13
  • $\begingroup$ Your answer is Brilliant @HenrikSchumacher, much more faster! Just one more doubt, how can I put all the time series together to plot or build histograms? $\endgroup$ – Herr Schrödinger Sep 4 '18 at 5:01
  • 1
    $\begingroup$ Gern geschehen, Herr Schrödinger ;) For joining these datasets, try something like TemporalData[Join @@ Through[aa["ValueList"]], {aa[[1]]["Times"]}]. Of course, this is meant to work only if the datasets share the same vector of time steps. $\endgroup$ – Henrik Schumacher Sep 4 '18 at 5:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.