I am trying to simulate a simple stochastic pde using the DifferentialEquations package, the equation I want to simulate is a simple conservation equation with conservative noise in 2 dimensions, i.e.

\partial_t \rho = D \nabla^2 \rho + \xi ; where the correlation of the noise \xi is

<\xi(x,t) \xi(x^\prime,t^\prime)> = -D \nabla^2 \delta(x-x^\prime) \delta (t-t^\prime).

Here the \nabla^2 in the noise correlation function stems from the conservation of the total density.

As a starting point I decided to adopt a very naive strategy, namely define \rho on a fixed N*N lattice, hence (\rho(x,y) is a matrix \rho_{ij}), and compute the Laplacian numerically. Defining an appropriate matrix K, the discrete Laplacian can be written as \rho.K+ K.\rho.

Though, I do not now how include the Laplacian in the correlation of the noise, probably I have to define a new noise dWc = dW.K+K.dW, where dW is matrix of N*N Wiener Processes. However I am not sure how to do that. From the documentation of the DifferentialEquations package it is clear how to implement colored noise or non diagonal noise (namely implementing something like M.dW); but I didn’t find a way to implement something like dWc = dW.K+K.dW. In general what I would like to do is to define a new noise which is a function of the old one dWnew = F(dW), is that possible?