let’s say I want to solve the following system of SDEs:

\begin{align*}
du_{1} &= f_{1}(u,du,p,t) dt+c_{11}dW_{1}+c_{12}dW_{2} \\
du_{2} &= f_{2}(u,du,p,t) dt+c_{23}dW_{3}+c_{24}dW_{4} \\
du_{3} &= f_{3}(u,du,p,t) dt+c_{35}dW_{5}+c_{36}dW_{6} \\
du_{4} &=f_{4}(u,du,p,t) dt+ c_{47}dW_{7}+c_{48}dW_{8}
\end{align*}

In Julia, the implementation of the g function (which gives the change of u due to the noise) looks like this:

function g!(u,du,t,p)

du[1,1] = c_{11}

du[1,2] = c_{12}

du[1,3] = 0

du[1,4] = 0

du[1,5] = 0

du[1,6] = 0

du[1,7] = 0

du[1,8] = 0

du[2,1] = 0

du[2,2] = 0

du[2,3] = c_{23}

du[2,4] = c_{24}

du[2,5] = 0

du[2,6] = 0

du[2,7] = 0

du[2,8] = 0

du[3,1] = 0

du[3,2] = 0

du[3,3] = 0

du[3,4] = 0

du[3,5] = c_{35}

du[3,6] = c_{36}

du[3,7] = 0

du[3,8] = 0

du[4,1] = 0

du[4,2] = 0

du[4,3] = 0

du[4,4] = 0

du[4,5] = 0

du[4,6] = 0

du[4,7] = c_{47}

du[4,8] = c_{48}

end

and then I define the problem and solve it as:

prob=SDEProblem(f!,g!,ones(4),(0.0,1.0),noise_rate_prototype = zeros(4,8))

where I assumed the initial condition u_0=[1 1 1 1 ], and I let the system evolve form 0.0 to time 1.0. In this example, inside g, I define a matrix 4x8 (containing the noise rates), such that g*dw gives the change du due to the stochastic part of the equation. since du is a 1d vector of length 4, dW is a vector of independent noises of length 8.

This is fine.

But I am wondering if there exist the possibility in julia to have instead of a 1d vector dW, a matrix. This would be useful in the case where the object u that I’m solving for is a matrix and the different components get independent noises.

For example let’s say u and du are 2x2 matrices. Then, is there a way to have a g function given by a 2x4 matrix and the dW by a 4x2 matrix so that g*dW is 2x2 and it correctly gives independent noises for the different comonents of du?