# How to deal with matrices involved in system of SDEs?

This question is in continuation of the the previous posts Solving Stochastic differential equation and Fast Simulations with Compile. What I want to do is numerically solving the epidemic model which can be defined as a system of stochastic differential equations (given below) by Euler-Maruyama method: $$\begin{array}{lll} dx_1&=& (-m_{12}x_1+m_{21}x_2)dt+\sqrt{\frac{m_{12}x_1+m_{21}x_2}{2}}(dw_1-dw_2),\\ dx_2&=&(m_{12}x_1-m_{21}x_2)dt+\sqrt{\frac{m_{12}x_1+m_{21}x_2}{2}}(-dw_1+dw_2). \end{array}$$ It can also be written in a more similar form as follows $$\left(\begin{array}{c}dx_1 \\dx_2\end{array}\right)=\left(\begin{array}{c}-m_{12}x_1+m_{21}x_2 \\m_{12}x_1-m_{21}x_2\end{array}\right)dt +\left( \begin{array}{cc} m_{12} x_1+m_{21} x_2 & -m_{12} x_1-m_{21} x_2 \\ -m_{12} x_1-m_{21} x_2 & m_{12} x_1+m_{21} x_2 \end{array} \right)^{1/2}\left(\begin{array}{c}dw_1 \\dw_2\end{array}\right),$$ where $0\leq t\leq T=100$, $x_1(0)=950$, $x_2(0)=50$, $m_{12}=\frac{0.04x_2}{x_1+x_2}$ and $m_{21}=0.01$. $w_1$ and $w_2$ are two independent standard Wiener processes.

I want to implement EM method or Milsetin method for finding some sample paths of this It$\hat{\text{o}}$ SDE, but there is a need for computing the (principal) matrix square root at each time discretization. For this, one may use the following code for computing matrix functions:

FunM[fun_, X_] :=
Module[{faux, dim, mataux, JordanD, sim, JordanF, eps, fdiag, diagQ,
fauxD}, (dim = Length@X;
faux[xx_, i_, j_] :=
Which[i <= j, 1/Abs[i - j]! (D[fun, {x, Abs[i - j]}]) /. x -> xx,
True, 0];
mataux[Y_] :=
Table[faux[Y[[i, j]], i, j], {i, 1, dim}, {j, 1, dim}];
JordanD = JordanDecomposition[X] // N; sim = JordanD[[1]];
JordanF = JordanD[[2]];
eps = 1*10^-10;
diagQ = Norm[JordanF - DiagonalMatrix[Diagonal[JordanF]]];
fauxD[xx_] := (fun) /. x -> xx;
fdiag := DiagonalMatrix[Map[fauxD, Diagonal[JordanF]]];
Which[diagQ < eps, sim.fdiag.Inverse[sim], True,
sim.mataux[JordanF].Inverse[sim]])]


Although I used this, I failed to implement EM and Milestin methods. I will be grateful if someone put some hints or answers for implementing EM for the above SDE.

-
–  Fazlollah Soleymani Aug 28 '14 at 14:27
It is now done. Can anyone give even some general hints or codes for dealing with system of SDEs? –  Fazlollah Soleymani Aug 29 '14 at 7:19