is it efficient to define random number generator inside function integrator?

user_784433 0
Viewed 59

Consider an Euler integrator that solves a stochastic differential equation:

void euler(vector<double> &x0,
           vector<double> &dxdt,
           const double dt)
{
    std::random_device rd;
    std::mt19937 rng(rd());
    std::normal_distribution<> dist(0, 1);
    f(dxdt, t, dt)        
    for (int i=0; i<x0.size(); ++i)
       x0[i] += dxdt[i] * dt + sqrt(dt) * 0.01 * dist(rng);
}

is this efficient to define the random generator for each time step integration? and probably there is a better option? another problem with this method is that when I try to fix the random seed

const unsigned int seed = 2;
std::mt19937 rng(seed);

for each time step, I get the same random numbers and this affects the answer.

c++
1 Answers

Seeding a PRNG if often costly and should usually only be done once during the whole program run so, no, this is not efficient.

I suggest that you break the creation of the PRNG out into a separate function that has a static PRNG (only initialized once).

std::mt19937& rng() {
    static std::mt19937 instance{std::random_device{}()};
    return instance;
}

or, if the PRNG is going to be used from multiple threads simultaneously:

std::mt19937& rng() {
    thread_local std::mt19937 instance{std::random_device{}()};
    return instance;
}

You can then use rng() in every function that needs a PRNG:

void euler(std::vector<double>& x0,
           std::vector<double>& dxdt,
           const double dt) 
{
    // Now use rng() instead of rng in here
}
user_75822470
Related