Matrix SDEs

Specific Domains Modelling & Simulations
1

Hi all,

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


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?


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
2

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.

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
1 Like
3

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.

1 Like