Stochastic diff eq driven by 1/f noise

Specific Domains Numerics
,
1

Hi all, I’m new to DifferentialEquations.jl and looking for tips on how model a stochastic process driven by noise that has a 1/f spectrum. Admittedly, this question might be not only on how to use the Julia stochastic differential equations solvers, but also a more general question on how to model 1/f noise in a stochastic diff eq. More precisely, if we have a stochastic diff eq

d u(t) / dt = f(u, t) + g(u, t) xi(t),

with xi(t) the noise term, then xi(t) should have a spectrum S(f) ~ 1/f.

After reading the section on Noise Processes in the DifferentialEquations.jl docs, one possible way seems to be that I generate a discrete grid of sampled points of xi(t_j) and use the NoiseGrid type. I suppose then I should choose the step size of the solver to coincide with my grid? However I wonder if there is a better way. As 1/f noise is so common in physics and engineering I figured there might be someone out there that has already done this and has a good example on hand (using a pre-sampled grid of points or otherwise). Thanks in advance.

2

1/f noise obviously has issues near DC, f=0 what do you want the spectrum to look like exactly.

3

Yeah right, good point. I will assume a low-frequency cutoff at f_ir. The physical experiment should ultimately be insensitive to the precise shape below this cutoff, so it should be fine to take it e.g. to be a hard cutoff with S(f) = 0 of S(f) = const. for f < f_ir. If the hard cutoff is harder to model and there is some other analytical form that makes the modeling easier I’m fine with that too.

4

Just define an SDE for 1/f noise and then use that as the process in another part of the system.

5

Thanks Chris. That sounds like an elegant way to do it.

I have been doing some research into the paper you posted as well as other papers by the same authors. In particular I implemented the diff eq. given in Eq. (14) of this follow up: [1111.2995] Tsallis distributions and 1/f noise from nonlinear stochastic differential equations, i.e.,

dx = \sigma^2 ( \eta - \lambda/2) (x^2 + x_0^2)^{\eta - 1} x dt + \sigma (x^2 + x_0^2)^{\eta/2} dW

Unlike the other equations in this paper and the one you linked to, this is the only equation that allows negative x (according to the authors). I would like to have a noise process that describes symmetric fluctuations around a mean value, so this equation seems to be the best choice, for that reason.

However, I have some challenges in the numerical implementation. Occasionally, this equation will produce solutions with sudden extremely large (positive or negative) values, which causes very slow integration. My current implementation is so slow that it is not practical; I can’t produce a decent number of trajectories in a reasonable time. I’m wondering if the situation can be improved by choosing better solver options, or perhaps modifying the equation to prevent very large (absolute) values, without modifying the noise spectrum. Here is a minimal code example, with the solver options I am using:

λ = 3  # 1/f case

σ = 1.0
η = 5/2  # somewhat arbitrary choice, same as in the article
x0 = 1.0

f(x, p, t) = σ^2 * (η - λ/2) * (x0^2 + x^2)^(η - 1) * x
g(x, p, t) = σ * (x0^2 + x^2)^(η / 2)

prob = SDEProblem(f, g, 0.0, (0.0, 1.0))
sol = solve(prob, SOSRI(), maxiters=1e9, dtmin=1e-16);

In the article I linked to above, they don’t give a lot of details about the numerical solver, but in [1003.1155] Modeling scaled processes and 1/f^b noise by the nonlinear stochastic differential equations [Eq. (36)], the use an x-dependent time-step \Delta t = \kappa^2 / (x_0 + x)^l, with \kappa \ll 1 and l parameters. I also implemented such a function using a StepsizeLimiter (let me know if this looks right):

κ = 0.1
l = 2η - 2
dtFE(x, p, t) = abs(κ^2 / (x0 + x)^l)  # function to set step size with StepsizeLimiter
step_size_limiter = StepsizeLimiter(dtFE; safety_factor=1, max_step=true)
sol = solve(prob, EM(), maxiters=1e9, callback=step_size_limiter, adative=false, dt=1e-5);

seems to work, but I regularly get a warning Warning: Instability detected. Aborting. Not sure if I could simply discard the aborted trajectories, as one possible work around.

6

Rejecting trajectories could bias the distribution, so be careful about doing that.

adaptive.