# Matrix SDEs

**URL:** <https://discourse.julialang.org/t/matrix-sdes/89038>\
**Category:** Modelling & Simulations\
**Created:** [October 20, 2022, 9:55pm UTC](https://discourse.julialang.org/t/matrix-sdes/89038 "2022-10-20T21:55:44Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![simoneb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simoneb/32/43744_2.png) [@simoneb](https://discourse.julialang.org/u/simoneb)\
**Post date:** [October 20, 2022, 9:55pm UTC](https://discourse.julialang.org/t/matrix-sdes/89038/1 "2022-10-20T21:55:44Z")

</div>

Hi all,

I am trying to fit the following matrix SDE using the following code:

```nohighlight

G = [0 -1; 1 0];
p = (A=G, B=G)
function drift(du, u, p, t) ## drift function for the SDE
    du = t * p.A;
end # drift

function diffusion(du, u, p, t) ## diffusion function for the SDE
    du = t * p.B;
end

time_tot = 1.0;
tspan = (0.0, time_tot);
u0 = reshape(zeros(4), (2,2))

prob = SDEProblem(drift, diffusion, u0, tspan, p=p); ## setup SDE problem
sol = solve(prob, EM(), p=p, dt=0.001); ## Solve using E-M scheme NOT WORKING!

```

This doesn’t work at all. No errors are thrown but the solution never moves from its initial condition u0. Can someone please tell me how to correctly code this matrix SDE, or suggest some reading material for me. Is a mass matrix required? If so, what would it look like and how do I include it in the problem formulation?

Thanks,  
Simone.

Edit: What I really want is to constrain the solution to be skew-symmetric. How do I do that? The following code works but the solutions are not skew-symmetric. I guess you need the mass matrix to provide the constraints? If so, how?

```nohighlight

G = [0 -1; 1 0];
p = (A=G, B=G);
function drift(du, u, p, t) ## drift function for the SDE
    du[1] = t * p.A[1,1]; 
    du[2] = t * p.A[1,2];
    du[3] = t * p.A[2,1];
    du[4] = t * p.A[2,2];
end # drift

function diffusion(du, u, p, t) ## diffusion function for the SDE
    du[1] = t * p.B[1,1]; 
    du[2] = t * p.B[1,2];
    du[3] = t * p.B[2,1];
    du[4] = t * p.B[2,2];
end

time_tot = 1.0;
tspan = (0.0, time_tot);
u0 = reshape(zeros(4), (2,2));
prob = SDEProblem(drift, diffusion, u0, tspan, p=p); ## setup SDE problem
sol = solve(prob, p=p, dt=0.001);
sol.u
```

---

<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:** [October 21, 2022, 5:55am UTC](https://discourse.julialang.org/t/matrix-sdes/89038/2 "2022-10-21T05:55:40Z")

</div>

You defined a mutating form of the function but then you didn’t mutate, so the values are still zero. Actually mutate via `.=`, i.e.

```julia
function drift(du, u, p, t) ## drift function for the SDE
    du .= t .* p.A;
end # drift

function diffusion(du, u, p, t) ## diffusion function for the SDE
    du .= t .* p.B;
end

```

---

<div class="post-metadata">

**Author:** ![simoneb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simoneb/32/43744_2.png) [@simoneb](https://discourse.julialang.org/u/simoneb)\
**Post date:** [October 21, 2022, 6:20am UTC](https://discourse.julialang.org/t/matrix-sdes/89038/3 "2022-10-21T06:20:47Z")

</div>

Yep. I figured that out. Sorry!! The other problem was the skew-symmetric constraint. One of the symplectic solvers works well. It took a lot of trial and error. I am relatively new to julia and to DifferentialEquations.jl, so thanks for your patience.

Thanks for all your work and your help.

Simone.