Does DifferentialEquations.jl support matrix-shaped noise

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}

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?

At this time it’s not supported. You’d have to reshape in the g function, but that reshape is free

1 Like

But how can I do the reshape inside the G? If G returns a vector then I think every component will just get a single noise, hence I obtain a diagonal noise.