Hello!

Let’s say I want to solve numerically the following stochastic partial differential equation:

\begin{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))
\end{equation},

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:

\begin{equation}
du_{stochastic}=G*dW
\end{equation}

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:

\begin{equation}
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}
\end{equation}.

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!