# System of SDE with complex noise

Hello!
Let’s say I want to solve numerically the following stochastic partial differential equation:
$$du(x,y,t)=\nabla^{2}u(x,y,t) dt+k\cdot dW^{r}(x,y,t)+i k\cdot dW^{c}(x,y,t))$$,
where u is a complex function in 2d, i is the imaginary unit, and k is a constant.
As one can see, we have a complex noise. In particular, the Winer process has real part W^{r} and imaginary part W^{c}, which are independent random variables.
Now, in order to solve the equation, I discretize in space, such that I obtain the following system of equations:

du[i,j]=F_{i,j}*u[i,j]dt+k*(dW^{r}[i,j]+idW^{c}[i,j])

I didn’t explicitely write the drift part (the laplacian) to make it shorter, and also because I want to concentrate on the stochastic part, as I think it is the cause of my problems.
So let me rewrite just the stochastic part of the system of equations:

du[i,j]=k\cdot dWr[i,j]+i k \cdot dWc[i,j]
(the “i” multiplying k outside the squared brackets is the imaginary unit and not an index!).
Now, reading the page of DifferentialEquations.jl ()DifferentialEquations.jl documentation, I understood that julia handles the stochastic part of an SDE in the following form:
$$du_{stochastic}=G*dW$$
where, if u is a vector, then G is a matrix such that the matrix-vector product G*dW gives the stochastic change of the solution du. (So I guess the solver multiplies G by an array of independent noises dW)

Now, in my case u is a matrix, so in order to put myself in the same scenario described in the documentation, I reshaped u from a matrix (N \times M) into a column vector (N M \times 1) ( usually I have N=M=300).
Having done that, the G matrix in my case will be N M \times 2N M, such that, if I multiply it by a vector of independent noises dW of dimension 2N M \times 1, I get a vector of dimension N M \times 1 corresponding to du. In particular, by looking at my system of equations, my G matrix will have the following form:
$$g=\begin{bmatrix} k & i \cdot k & 0 & 0 & 0 & 0 & \dots \\ 0 & 0 & k & i \cdot k& 0 & 0 & \dots \\ 0 & 0 & 0 & 0 & k & i \cdot k & \dots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \\ \end{bmatrix}$$.

Therefore, I implemented the noise prototype as follows:

function create_matrix_G_prototype(imax,jmax,p)
g_matrix_prototype=zeros(Complex{Float64},N*M,2*N*M)
for i in 1:imax*jmax
g_matrix_prototype[i,(2*i)-1]=k
g_matrix_prototype[i,(2*i)]=im*k
end
return  sparse(g_matrix_prototype)
end

G_matrix_prototype=create_matrix_G_prototype(imax,jmax,p)


In this way, I can implement the G function to pass to the solver as an in-place function:

function G!(du, u, p, t)
du.=G_matrix_prototype
ende


Then I define the problem and solve it with:

u_0=zeros(Complex{Float64},N*M)
tspan=(0.0,1.0)
problem_r=SDEProblem(F!,G!,u_0,tspan,p,noise_rate_prototype=G_matrix_prototype)
sol_r=solve(problem_r)e


The issue here is that, first of all, the code seems to be terribly slow. Moreover, I get a warning telling me that the dt chosen by the solver is smaller than the dt_{min}:

┌ Warning: dt(8.881784197001252e-16) <= dtmin(8.881784197001252e-16) at t=0.37906220534331775, and step error estimate = 1.1208004104491986. Aborting. There is either an error in your model specification or the true solution is unstable.


This is surprising to me since I didn’t expect the equation to be stiff (but at this point I’m starting to guess it is). Has anyone any idea of what is going wrong? what shall I do about the warning? maybe choose a specific solver for stiff equations?
Also, is there a better approach that I could take for this specific case in order to increase the performance and speed of the code? Is it a good idea to separate real and imaginary parts instead of dealing with complex values?

Thanks!

Did you try changing this to just a simple non-adaptive algorithm? I don’t think the implementation is doing anything bad persay here, but the adaptivity on complex SDEs would have to take into account the non-commutative nature of complex noise which may be messing up the error estimators.