I have the following SDE

du = f(u,p,t)dt + sum_i g_i(u,p,t) dW_i

where du, f, and g_i are all matrices, while W_i are all independent Wiener processes. I have no problem in setting up the ODE when g_i is zero, but I can’t figure out how to write the above problem in DifferentialEquations.jl in the general case, where all functions g_i have to be packed in a single function g.

Is g going to be a vector of matrices (each element being g_i), or a 3-index tensor such that g[:,:,i] = g_i[:,:] ? In any case it is not clear how can I set the noise_rate_prototype, or whether I should use a completely different approach.

Any ideas?

Thanks

# SDE with matrix values, question on DifferentialEquations.jl

**benkj**#1

Did you see this part of the tutorial which “translates” between the matrix and summation forms?

**benkj**#3

Yes, indeed my “poor man” solution is based on vectorising everything: matrix u is transformed into a vector, so is f and g_i. In that case g becomes a matrix (whose columns are g_i), so I can use the example that you mentioned. However, in f and g I have to “devectorise” u to transform it back into a matrix, are my f and g are doing some non-trivial matrix operations that are not straightforward in the vectorised notation. So I think that my poor-man solution is not optimal, and I don’t see any reason why in julia it shouldn’t be possible to work directly with matrix-valued functions by playing with data types, rather than using my ugly trick of vectorising and devectorising the function each time.

Multiplication between a 3-tensor and 2-tensor isn’t well-defined though. Generalizing this definition in a way that doesn’t require a bunch of `g`

functions passed in is hard, but having a bunch of functions rather than one in-place function would be very very inefficient. That’s the main issue that keeps it a noise matrix for now. But you can define it as any type you want and then define `*`

against `dW`

to have the correct interpretation. That would be a useful `AbstractArray`

to work out and add to the library if you’re up for it.

Using `reshape`

is cheap/free though since it’s just a view. There shouldn’t be a cost to this.

**benkj**#5

Thanks, I like your idea of defining * against a type and dW. I will see if I manage do to something, but otherwise I’ll keep my poor man solution.

Cheers