Random forcing

[Anna]

Hi Users,

I'm quite new to Dedalus.

While setting up a simple 1D system of nonlinear equations that are randomly forced, I've come across a weird issue with the forcing. To simplify, I use a General Function to define a forcing that could be time- or x- dependent and has a random amplitude that is close to 1, for example:

def my_forcing(*args):

    s = 0.1; mu = 1.0

    return np.random.normal(mu,s)

def Forcing(*args, domain=domain, F=my_forcing):

    return de.operators.GeneralFunction(domain, layout='g', func=F, args=args)

de.operators.parseables['F'] = Forcing

and later add it to my system with an additional equation:

problem.add_equation("u = F()")

I integrate the system in time and record time t and u that is applied inside the solver. After that, I compare u with a test forcing u_test = myforcing() computed independently after the simulation. Of course, they are not expected to be the same since the forcing is random, but surprisingly, the standard deviation of u recorded in the simulation is much larger than that of the a posteriori u_test (see the figure attached). Based on the parameters of my forcing (s = 0.1; mu = 1.0), u_test seems to be correct, and u wrong. If the randomness is removed from the forcing, everything works just fine and u and u_test are identical. 

What could have gone wrong here? I tested various variants of the forcing and came up with a minimal script that demonstrates the problem (also attached). 

I would be grateful for any help on this!

Thanks,

Anna

[Anna]

Hi all,

it seems that this issue is related to the timestepper. Timesteppers like, for example, SBDF2,

... =  (1 + \omega) F(Un+1) − \omega F(Un) + ...                     (eqn. 2.8 from Wang 2008), 

or RK443, will sum random forcing at more than one timestep / RK iteration. The resulting standard deviation will be higher. In the case of SBDF2, it is higher by about \sqrt(5), a factor that can be derived from the expression above if \omega=1 (constant step-size). 

A straightforward solution is to use schemes with 1st-order explicit terms, but besides that, is there a better way to incorporate random forcing into the solver? For more complex problems, a 1st-order explicit scheme may be not good enough. 

Best,

Anna

[Jeffrey S. Oishi]

Hi Anna,

Higher order stochastic timestepping is a sticky issue, one that is not (to my knowledge) well agreed upon in the stochastic differential equations community. As such, my advice would be to use first-order Euler-Maruyama timestepping if you have truly stochastic forcing. 

Perhaps others on the list will have other ideas, but I think any of the higher-order timesteppers in Dedalus will suffer from the same issue you have identified here.

Jeff

--
You received this message because you are subscribed to the Google Groups "Dedalus Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dedalus-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dedalus-users/45392e84-f200-4248-973e-a98422cbbb82n%40googlegroups.com.

[Anna]

Hi Jeff,

thank you for the comment, I see now that it is a tricky issue. I'll use a first-order scheme for the time being.

Anna