# System of SDE with complex noise

**URL:** <https://discourse.julialang.org/t/system-of-sde-with-complex-noise/111469>\
**Category:** Performance\
**Tags:** sde, differentialequation\
**Created:** [March 11, 2024, 3:23pm UTC](https://discourse.julialang.org/t/system-of-sde-with-complex-noise/111469 "2024-03-11T15:23:21Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![matt1](https://avatars.discourse-cdn.com/v4/letter/m/b782af/32.png) [@matt1](https://discourse.julialang.org/u/matt1)\
**Post date:** [March 11, 2024, 3:23pm UTC](https://discourse.julialang.org/t/system-of-sde-with-complex-noise/111469/1 "2024-03-11T15:23:21Z")

</div>

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](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/sde_example/#Example-4:-Systems-of-SDEs-with-Non-Diagonal-Noise), 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:

```julia
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:

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

```

Then I define the problem and solve it with:

```julia
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}:

```julia
┌ 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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 15, 2024, 11:00am UTC](https://discourse.julialang.org/t/system-of-sde-with-complex-noise/111469/2 "2024-06-15T11:00:46Z")

</div>

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.