euler maruyama scheme for stochastic delay differential equation

[Ziaeemehr]

Hi, As I know there is only Runke kutta 23 method available for DDE. Is it possible for you to add also Euler maruyama or if you don't have time give some guide how to add it?

best, Abolfazl

[Wrzlprmft]

Adding a new integrator to JiTCDDE is far from easy, because the existing integrator is strongly interwoven with the algorithm to gain efficiency. Also, I would be quite hesitant to add any non-adaptive integrators. I see two general options for you:

• Implement Euler–Maruyama with a fixed step size for DDEs. It’s not that difficult to do, since many of the things that make JiTCDDE complex are gone (storing a variable-step past, interpolating it, error estimates, and step-size adaption). Of course, it will be slow and you have to ensure that your time step is sufficiently small yourself (instead of relying on step-size adaption).

• Use jitcdde_input to define your noise(s) as an input. This limits the frequency of the noise (i.e., requires some smoothness), but depending on your application this may be more realistic than white noise anyway. If you do this, ensure that max_step is not larger than the time scale of your noise. This example essentially does this (although the random data is thought as a placeholder for real data there). Technically, you can also do this with callback_functions, but this requires more care and is slower, so I advise against this.

Also see Issue #27.

If you can tell me more about the nature and purpose of your noise, I may also give more specific advice.

[Ziaeemehr]

Hi, For some optimization reason, I implemented the SDDE code in C++. this is not directly an issue about jitcdde, but I hope to get some hints if you could find some free time.

I show some simple implementation in python just for demonstration which is very slow, but I actually need to get some hints, on what can I use in C/C++.

In the Euler scheme, delayed points can be accessed by index.

system of equations are:
    dy0dt = -y0 + y1(t-1) + y1(t-10) + sigma * dW
    dy1dt = y0 + y1(t-1) - y1
    dy2dt = y1 - y1(t-10)

D = [int(delays[i]/dt) for i in range(len(delays))]

def f_sys_index(t, n):
        f0p = -y[0][n] * y[1][n-D[0]] + y[1][n-D[1]]
        f1p = y[0][n] * y[1][n-D[0]] - y[1][n]
        f2p = y[1][n] - y[1][n-D[1]]

def euler(f_sys):
        for _ in range(nstart, nstart+num_iterations):
            dy = dt * f_sys(t[-1], -1)
            append(y[:, -1] + dy, t[-1] + dt)

To reduce the memory usage I keep only int(maxdelay/dt) elements for each component and roll the coordinates at each time step.

def append(new_y, new_t):
        """append an element"""
        y = np.roll(y, -1)
        y[:, -1] = new_y
        t = np.roll(t, -1)
        t[-1] = new_t

does jitcdde use a circular buffer? thank you for any guide in advance.

[Wrzlprmft]

JiTCDDE uses a doubly-linked cursor list. This is necessary because of the adaptive step size, covering for time-dependent delays, etc. (see Appendix A of the paper for details). If you are working with a constant step size and fixed delays, a circular buffer is a reasonable choice.