Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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, 

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.

share|improve this question
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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.