Extension of: Numerically solving a system of SDE's with Levy noise?

Question

A great answer by Alex is to be found here for my original question: Numerically solving a system of SDE's with Levy noise?

Now Let us perturb this system with time delays so the system is:

\begin{align} dS&=\big[A-(\mu_1+\kappa) S(t)-\beta S(t) D(t)+ \kappa S(t-\tau_1)e^{-\mu\tau_1}+\gamma I(t-\tau_2)e^{-\mu\tau_2}\big]dt+\sigma_1 S(t) dW_1(t)+\int_U \lambda_1(u)S(t-)\tilde{N}(dt,du)\\[1ex] dI&=\big[\beta S(t) D(t)-(\mu_2+\gamma)I(t)\big]dt+\sigma_2 I(t) dW_2(t)+\int_U \lambda_2(u)I(t-)\tilde{N}(dt,du)\\[1ex] dD&=\big[\eta(I(t)-D(t)]dt+\sigma_4 D(t) dW_4(t)+\int_U \lambda_4(u)D(t-)\tilde{N}(dt,du) \end{align}

where $\kappa$ is a postive and $\tau_1$ and $\tau_2$ are nonnegative with $t\geq -\tau$, where, $\tau=\max\lbrace \tau_1, \tau_2\rbrace$.

Code:

SeedRandom[1234];
A (*Recruitment rate*)= {0.9, 0.3, 0.6, 0.6};
µ1 (*Natural mortality rate of S*)= {0.3, 0.3, 0.4, 0.4};
\[Beta] (*Transmission rate*)= {0.07, 1.3, 0.35, 0.8};
\[Gamma] (*Recovered rate*)= {0.05, 0.05, 0.2, 0.3};
µ2 (*General mortality of I*)= {0.5, 0.5, 0.3, 0.3};
\[Eta] (*Exponentially fading memory rate*)= {0.09, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= {0.15, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= {0.25, 0.25, 0.15, 0.15};
\[Sigma]4 (*Intensity of W4(t)*)= {0.27, 0.27, 0.13, 0.13};
\[Lambda]1 (*Jump intensity of S*)= {0.2, 0.2, 0.5, 0.5};
\[Lambda]2 (*Jump intensity of I*)= {0.23, 0.23, 0.3, 0.3};
\[Lambda]4 (*Jump intensity of D*)= {0.1, 0.1, 0.7, 0.7};
\[Kappa](*----*)= {0.2, 0.1, 0.7, 0.7};
\[Tau]1(*----*)= {1, 0.1, 0.7, 0.7};
\[Tau]2(*----*)= {1, 0.1, 0.7, 0.7};

tmax = 301; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1]; pWe2 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW2 = 
 Interpolation[Table[{(j - 1), pWe2[[j]]}, {j, Length[pWe2]}], 
  InterpolationOrder -> 1]; pWe4 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW4 = 
 Interpolation[Table[{(j - 1), pWe4[[j]]}, {j, Length[pWe4]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1.], {0, tmax}]; pL2 = 
 RandomFunction[PoissonProcess[1.1], {0, tmax}]; pL4 = 
 RandomFunction[PoissonProcess[.9], {0, tmax}]; ListStepPlot[{pL1, 
  pL2, pL4}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences; dpL2 = 
 pL2["SliceData", Range[tmax]] // First // Differences; dpL4 = 
 pL4["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;
L2[t_] := If[t < 1, 0, dpL2[[Round[t]]]]/tmax;
L4[t_] := If[t < 1, 0, dpL4[[Round[t]]]]/tmax;

eq1 = -s'[t] + (a - (mu1 + kappa) s[t] - beta s[t] d[t] + 
     kappa s[t - tau1] Exp[mu tau1] + 
     gamma i[t - tau2] Exp[mu tau2]) + sigma1 s[t] dW1[t] + 
   lambda1 s[t] L1[t];
eq2 = -i'[t] + (beta s[t] d[t] - (mu2 + gamma) i[t]) + 
   sigma2 i[t] dW2[t] + lambda2 i[t] L2[t];
eq3 = -r'[t] + (gamma i[t] - mu3 r[t]);
eq4 = -d'[t] + eta (i[t] - d[t]) + sigma4 d[t] dW4[t] + 
   lambda4 d[t] L4[t];

ic = {s[0] == 0.6, i[0] == 0.3, d[0] == 0.05};
rul[j_] := {a -> A[[j]], beta -> \[Beta][[j]], gamma -> \[Gamma][[j]],
    eta -> \[Eta][[j]], mu1 -> µ1[[j]], mu2 -> µ2[[j]], 
   sigma1 -> \[Sigma]1[[j]], sigma2 -> \[Sigma]2[[j]], 
   sigma4 -> \[Sigma]4[[j]], lambda1 -> \[Lambda]1[[j]], 
   lambda2 -> \[Lambda]2[[j]], lambda4 -> \[Lambda]4[[j]], 
   kappa -> \[Kappa][[j]], tau1 -> \[Tau]1[[j]], tau2 -> \[Tau]2[[j]]};

eqn[j_] := {eq1, eq2, eq4} /. rul[j];

sol[j_] := 
  NDSolve[{eqn[j] == {0, 0, 0}, ic}, {s, i, d}, {t, 0, tmax - 1}];

With[{sol = sol[1]}, {Plot[Evaluate[s[t] /. sol], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> Green, 
   FrameLabel -> {"Time t (Days)", "S"}], 
  Plot[Evaluate[i[t] /. sol], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> Blue, 
   FrameLabel -> {"Time t (Days)", "I"}], 
  Plot[Evaluate[d[t] /. sol], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> Red, 
   FrameLabel -> {"Time t (Days)", "D"}]}]

How can we numerically solve this?

Edit: This is the code for the second papers equations(see comments to Alex's answer):

A (*Recruitment rate*)= {1, 0.3, 0.6, 0.6};
mu0 (*Natural mortality rate of S*)= {0.09, 0.3, 0.4, 0.4};
\[Beta] (*Transmission rate*)= {0.18, 1.3, 0.35, 0.8};
\[Gamma] (*Recovered rate*)= {0.55, 0.05, 0.2, 0.3};
\[Delta] (*Recovered rate*)= {0.44, 0.05, 0.2, 0.3};
\[Alpha]1 (*General mortality of I*)= {0.4, 0.5, 0.3, 0.3};
\[Alpha]2 (*General mortality of I*)= {0.02, 0.5, 0.3, 0.3};
\[Epsilon] (*Exponentially fading memory rate*)= {0.6, 0.09, 0.7, 0.2};
\[Rho] (*Exponentially fading memory rate*)= {0.2, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= {0.2, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= {0.85, 0.25, 0.15, 0.15};
\[Sigma]3 (*Intensity of W4(t)*)= {0.5, 0.27, 0.13, 0.13};
\[Lambda]1 (*Jump intensity of S*)= {0., 0.2, 0.5, 0.5};
\[Lambda]2 (*Jump intensity of I*)= {0., 0.23, 0.3, 0.3};
\[Lambda]3 (*Jump intensity of D*)= {0., 0.1, 0.7, 0.7};
\[Tau]1 (*Jump intensity of S*)= {0.5, 0.2, 0.5, 0.5};
\[Tau]2 (*Jump intensity of I*)= {0.5, 0.23, 0.3, 0.3};
\[Tau]3 (*Jump intensity of D*)= {0.5, 0.1, 0.7, 0.7};

tmax = 5001; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1]; pWe2 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW2 = 
 Interpolation[Table[{(j - 1), pWe2[[j]]}, {j, Length[pWe2]}], 
  InterpolationOrder -> 1];
pWe3 = RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 
   1]]; dW3 = 
 Interpolation[Table[{(j - 1), pWe3[[j]]}, {j, Length[pWe3]}], 
  InterpolationOrder -> 1]; pWe4 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW4 = 
 Interpolation[Table[{(j - 1), pWe4[[j]]}, {j, Length[pWe4]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1.], {0, tmax}];
pL2 = RandomFunction[PoissonProcess[1.1], {0, tmax}];
pL3 = RandomFunction[PoissonProcess[1.1], {0, tmax}];
pL4 = RandomFunction[PoissonProcess[.9], {0, tmax}];
ListStepPlot[{pL1, pL2, pL3, pL4}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences;
dpL2 = pL2["SliceData", Range[tmax]] // First // Differences;
dpL3 = pL3["SliceData", Range[tmax]] // First // Differences;
dpL4 = pL4["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;
L2[t_] := If[t < 1, 0, dpL2[[Round[t]]]]/tmax;
L3[t_] := If[t < 1, 0, dpL3[[Round[t]]]]/tmax;
L4[t_] := If[t < 1, 0, dpL4[[Round[t]]]]/tmax;

eq1 = -s'[t] + (a - (mu + rho) s[t] - beta s[t] i[t] + 
     rho s[t - tau1] Exp[-mu tau1] + 
     gamma i[t - tau2] Exp[-mu tau2] + 
     epsilon q[t - tau3] Exp[-mu tau3]) + sigma1 s[t] dW1[t] + 
   lambda1 s[t] L1[t];
eq2 = -i'[t] + (beta s[t] i[t] - (mu + alpha1 + delta + gamma) i[t]) +
    sigma2 i[t] dW2[t] + lambda2 i[t] L2[t];
eq3 = -q'[t] + (delta i[t] - (mu + alpha2 + epsilon) q[t]) + 
    sigma3 q[t] dW3[t] + lambda3 q[t] L3[t];;
eq4 = -r'[t] + (gamma i[t] + rho s[t] + epsilon q[t] - mu r[t] - 
     rho s[t - tau1] Exp[-mu tau1] - gamma i[t - tau2] Exp[-mu tau2] -
      epsilon q[t - tau3] Exp[-mu tau3]) + sigma4 d[t] dW4[t] + 
   lambda4 d[t] L4[t];

ic1 = {s[t /; t <= 0] == 5, i[t /; t <= 0] == 15, q[t /; t <= 0] == 5};
rul[j_] := {a -> A[[j]], beta -> \[Beta][[j]], rho -> \[Rho][[j]], 
   gamma -> \[Gamma][[j]], epsilon -> \[Epsilon][[j]], 
   delta -> \[Delta][[j]], alpha1 -> \[Alpha]1[[j]], 
   alpha2 -> \[Alpha]2[[j]], sigma1 -> \[Sigma]1[[j]], 
   sigma2 -> \[Sigma]2[[j]], sigma3 -> \[Sigma]3[[j]], 
   sigma4 -> \[Sigma]4[[j]], lambda1 -> \[Lambda]1[[j]], 
   lambda2 -> \[Lambda]2[[j]], lambda3 -> \[Lambda]3[[j]], 
   lambda4 -> \[Lambda]4[[j]], tau1 -> \[Tau]1[[j]], 
   tau2 -> \[Tau]2[[j]], tau3 -> \[Tau]3[[j]], mu -> mu0[[j]]};

eqn[j_] := {eq1, eq2, eq3} /. rul[j]; 

sol1[j_] := 
  NDSolve[{eqn[j] == {0, 0, 0}, ic1}, {s, i, q}, {t, 0, tmax - 1}];

With[{sol = sol1[1]}, {Plot[Evaluate[s[t] /. {sol}], {t, 0, tmax - 1},
    PlotRange -> All, Frame -> True, PlotStyle -> {Green}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["S(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18], 
  Plot[Evaluate[i[t] /. {sol}], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> {Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["I(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18], 
  Plot[Evaluate[q[t] /. {sol}], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> {Red},
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["Q(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18]}]

EDIT 2:

   A (*Recruitment rate*)= {1, 0.3, 0.6, 0.6};
mu0 (*Natural mortality rate of S*)= {0.09, 0.3, 0.4, 0.4};
\[Beta] (*Transmission rate*)= {0.39, 1.3, 0.35, 0.8};
\[Gamma] (*Recovered rate*)= {0.55, 0.05, 0.2, 0.3};
\[Delta] (*Recovered rate*)= {0.44, 0.05, 0.2, 0.3};
\[Alpha]1 (*General mortality of I*)= {0.4, 0.5, 0.3, 0.3};
\[Alpha]2 (*General mortality of I*)= {0.02, 0.5, 0.3, 0.3};
\[Epsilon] (*Exponentially fading memory rate*)= {0.6, 0.09, 0.7, 0.2};
\[Rho] (*Exponentially fading memory rate*)= {0.2, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= {0.2, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= {0.4, 0.25, 0.15, 0.15};
\[Sigma]3 (*Intensity of W3(t)*)= {0.4, 0.27, 0.13, 0.13};
\[Sigma]4 (*Intensity of W4(t)*)= {0.5, 0.27, 0.13, 0.13};
\[Tau]1 (*Jump intensity of S*)= {0.5, 0.2, 0.5, 0.5};
\[Tau]2 (*Jump intensity of I*)= {0.5, 0.23, 0.3, 0.3};
\[Tau]3 (*Jump intensity of D*)= {0.5, 0.1, 0.7, 0.7};

tmax = 201; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1]; pWe2 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW2 = 
 Interpolation[Table[{(j - 1), pWe2[[j]]}, {j, Length[pWe2]}], 
  InterpolationOrder -> 1];
pWe3 = RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 
   1]]; dW3 = 
 Interpolation[Table[{(j - 1), pWe3[[j]]}, {j, Length[pWe3]}], 
  InterpolationOrder -> 1]; pWe4 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW4 = 
 Interpolation[Table[{(j - 1), pWe4[[j]]}, {j, Length[pWe4]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1.], {0, tmax}];
pL2 = RandomFunction[PoissonProcess[1.1], {0, tmax}];
pL3 = RandomFunction[PoissonProcess[1.1], {0, tmax}];
pL4 = RandomFunction[PoissonProcess[.9], {0, tmax}];
ListStepPlot[{pL1, pL2, pL3, pL4}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences;
dpL2 = pL2["SliceData", Range[tmax]] // First // Differences;
dpL3 = pL3["SliceData", Range[tmax]] // First // Differences;
dpL4 = pL4["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;
L2[t_] := If[t < 1, 0, dpL2[[Round[t]]]]/tmax;
L3[t_] := If[t < 1, 0, dpL3[[Round[t]]]]/tmax;
L4[t_] := If[t < 1, 0, dpL4[[Round[t]]]]/tmax;

eq1 = -s'[t] + (a - (mu + rho) s[t] - beta s[t] i[t] + 
     rho s[t - tau1] Exp[-mu tau1] + gamma i[t - tau2] Exp[-mu tau2] +
      epsilon q[t - tau3] Exp[-mu tau3]) + sigma1 s[t] dW1[t];
eq2 = -i'[t] + (beta s[t] i[t] - (mu + alpha1 + delta + gamma) i[t]) +
    sigma2 i[t] dW2[t];
eq3 = -q'[t] + (delta i[t] - (mu + alpha2 + epsilon) q[t]) + 
  sigma3 q[t] dW3[t]; eq4 = -r'[t] + (gamma i[t] + rho s[t] + 
    epsilon q[t] - mu r[t] - rho s[t - tau1] Exp[-mu tau1] - 
    gamma i[t - tau2] Exp[-mu tau2] - 
    epsilon q[t - tau3] Exp[-mu tau3]) + sigma4 d[t] dW4[t];

eq1o = -s'[t] + (a - (mu + rho) s[t] - beta s[t] i[t]) + 
   sigma1 s[t] dW1[t];
eq2o = -i'[
     t] + (beta s[t] i[t] - (mu + alpha1 + delta + gamma) i[t]) + 
   sigma2 i[t] dW2[t];
eq3o = -q'[t] + (delta i[t] - (mu + alpha2 + epsilon) q[t]) + 
  sigma3 q[t] dW3[t]; eq4o = -r'[t] + (gamma i[t] + rho s[t] + 
    epsilon q[t] - mu r[t]) + sigma4 d[t] dW4[t];

eq1d = -s'[t] + (a - (mu + rho) s[t] - beta s[t] i[t] + 
     rho s[t - tau1] Exp[-mu tau1] + gamma i[t - tau2] Exp[-mu tau2] +
      epsilon q[t - tau3] Exp[-mu tau3]);
eq2d = -i'[
     t] + (beta s[t] i[t] - (mu + alpha1 + delta + gamma) i[t]);
eq3d = -q'[t] + (delta i[t] - (mu + alpha2 + epsilon) q[t]);
eq4d = -r'[t] + (gamma i[t] + rho s[t] + epsilon q[t] - mu r[t] - 
     rho s[t - tau1] Exp[-mu tau1] - gamma i[t - tau2] Exp[-mu tau2] -
      epsilon q[t - tau3] Exp[-mu tau3]);

ic1 = {s[t /; t <= 0] == 10, i[t /; t <= 0] == 0.1, 
   q[t /; t <= 0] == 0.1};
rul[j_] := {a -> A[[j]], beta -> \[Beta][[j]], rho -> \[Rho][[j]], 
   gamma -> \[Gamma][[j]], epsilon -> \[Epsilon][[j]], 
   delta -> \[Delta][[j]], alpha1 -> \[Alpha]1[[j]], 
   alpha2 -> \[Alpha]2[[j]], sigma1 -> \[Sigma]1[[j]], 
   sigma2 -> \[Sigma]2[[j]], sigma3 -> \[Sigma]3[[j]], 
   sigma4 -> \[Sigma]4[[j]], tau1 -> \[Tau]1[[j]], 
   tau2 -> \[Tau]2[[j]], tau3 -> \[Tau]3[[j]], mu -> mu0[[j]]};

eqn1[j_] := {eq1, eq2, eq3} /. rul[j];
eqn2[j_] := {eq1d, eq2d, eq3d} /. rul[j];

sol1[j_] := 
  NDSolve[{eqn1[j] == {0, 0, 0}, ic1}, {s, i, q}, {t, 0, tmax - 1}, 
   MaxSteps -> 10^6];
sol2[j_] := 
  NDSolve[{eqn2[j] == {0, 0, 0}, ic1}, {s, i, q}, {t, 0, tmax - 1}, 
   MaxSteps -> 10^6];



With[{sol = sol1[1], 
  sol2 = sol2[1]}, {Plot[
   Evaluate[s[t] /. {sol, sol2}], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> {Green, Black}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["S(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, 
   PlotLegends -> 
    Placed[LineLegend[{Green, Black}, {"S(t) with delays", 
       "S(t) without noise"}, LegendFunction -> Framed], {.8, .75}]], 
  Plot[Evaluate[i[t] /. {sol, sol2}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Blue, Black}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["I(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, 
   PlotLegends -> 
    Placed[LineLegend[{Blue, Black}, {"I(t) with delays", 
       "I(t) without noise"}, LegendFunction -> Framed], {.8, .75}]], 
  Plot[Evaluate[q[t] /. {sol, sol2}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["Q(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black}, {"Q(t) with delay", 
       "Q(t) without noise"}, LegendFunction -> Framed], {.8, .75}]]}]

EDIT 3:

As mentioned in the comments of Alex's answer, when we simulate this paper: https://www.researchgate.net/profile/Dianli-Zhao/publication/332241627_Threshold_dynamics_of_the_stochastic_epidemic_model_with_jump-diffusion_infection_force/links/5ca85d24a6fdcca26d013e72/Threshold-dynamics-of-the-stochastic-epidemic-model-with-jump-diffusion-infection-force.pdf, our plots, although similar, do not show the effects of Levy noise as what paper shows, any idea why?

The code for new paper:

µ(*Natural mortality rate of S,I,*)= {0.05, 0.05, 0.05, 0.4};
\[Beta] (*Transmission rate*)= {0.3, 0.3, 0.3, 0.8};
\[Delta](*Transmission rate*)= {0.05, 0.05, 0.05, 0.8};
\[Gamma] (*Recovered rate*)= {0.1, 0.1, 0.1, 0.3};
\[Sigma]1 (*Intensity of W1(t)*)= {0.1, 0.1, 0, 0.169};
\[Lambda]1 (*Intensity of W2(t)*)= {0.2, 0, 0, 0.15};


tmax = 801; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1], {0, tmax}];
ListStepPlot[{pL1}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;

eq1 = -s'[t] + (mu  - mu s[t] - beta s[t] i[t]) - 
   sigma1 s[t] i[t] dW1[t] + lambda1  L1[t];
eq2 = -i'[t] + (beta s[t] i[t] - (mu + delta + gamma) i[t]) + 
   sigma1 s[t] i[t] dW1[t] + lambda1  L1[t];

ic = {s[0] == 0.5, i[0] == 0.1};
rul[j_] := {beta -> \[Beta][[j]], gamma -> \[Gamma][[j]], 
   mu -> µ[[j]], delta -> \[Delta][[j]], sigma1 -> \[Sigma]1[[j]], 
   lambda1 -> \[Lambda]1[[j]]};

eqn[j_] := {eq1, eq2} /. rul[j];

sol[j_] := NDSolve[{eqn1[j] == {0, 0}, ic}, {s, i}, {t, 0, tmax - 1}];

With[{sol = sol[1], sol2 = sol[2], 
  sol3 = sol[3]}, {Plot[
   Evaluate[s[t] /. {sol, sol2, sol3}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black, Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["S(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, FrameStyle -> Black, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black, Blue}, {"With jumps", 
       "without jumps", "Deterministic"}, 
      LegendFunction -> Framed], {.85, .85}]], 
  Plot[Evaluate[i[t] /. {sol, sol2, sol3}], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> {Red, Black, Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["I(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18, FrameStyle -> Black, 
   PlotLegends -> 
    Placed[LineLegend[{Red, Black, Blue}, {"With jumps", 
       "without jumps", "Deterministic"}, 
      LegendFunction -> Framed], {.8, .8}]]}]

Thanks for the code, now you "This code yields many errors.", and... you can't be bother to find them? you don't understand them? They are inconsistent? Please do your part, show some diligence. What are these errors? — rhermans, Commented Jan 24, 2023 at 13:54

Math · Accepted Answer · 2023-01-28 05:28:21Z

With new parameter mu definition and initial condition for delayed equations we have

SeedRandom[1234];
A (*Recruitment rate*)= {0.9, 0.3, 0.6, 0.6};
µ1 (*Natural mortality rate of S*)= {0.3, 0.3, 0.4, 0.4};
\[Beta] (*Transmission rate*)= {0.07, 1.3, 0.35, 0.8};
\[Gamma] (*Recovered rate*)= {0.05, 0.05, 0.2, 0.3};
µ2 (*General mortality of I*)= {0.5, 0.5, 0.3, 0.3};
\[Eta] (*Exponentially fading memory rate*)= {0.09, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= {0.15, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= {0.25, 0.25, 0.15, 0.15};
\[Sigma]4 (*Intensity of W4(t)*)= {0.27, 0.27, 0.13, 0.13};
\[Lambda]1 (*Jump intensity of S*)= {0.2, 0.2, 0.5, 0.5};
\[Lambda]2 (*Jump intensity of I*)= {0.23, 0.23, 0.3, 0.3};
\[Lambda]4 (*Jump intensity of D*)= {0.1, 0.1, 0.7, 0.7};
\[Kappa](*----*)= {0.2, 0.1, 0.7, 0.7};
\[Tau]1(*----*)= {0.7, 0.1, 0.7, 0.7};
\[Tau]2(*----*)= {0.7, 0.1, 0.7, 0.7};
mu0 = {1, 1, 1, 1};
tmax = 301; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1]; pWe2 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW2 = 
 Interpolation[Table[{(j - 1), pWe2[[j]]}, {j, Length[pWe2]}], 
  InterpolationOrder -> 1]; pWe4 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW4 = 
 Interpolation[Table[{(j - 1), pWe4[[j]]}, {j, Length[pWe4]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1.], {0, tmax}]; pL2 = 
 RandomFunction[PoissonProcess[1.1], {0, tmax}]; pL4 = 
 RandomFunction[PoissonProcess[.9], {0, tmax}]; ListStepPlot[{pL1, 
  pL2, pL4}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences; dpL2 = 
 pL2["SliceData", Range[tmax]] // First // Differences; dpL4 = 
 pL4["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;
L2[t_] := If[t < 1, 0, dpL2[[Round[t]]]]/tmax;
L4[t_] := If[t < 1, 0, dpL4[[Round[t]]]]/tmax;

eq1 = -s'[t] + (a - (mu1 + kappa) s[t] - beta s[t] d[t] + 
     kappa s[t - tau1] Exp[mu tau1] + 
     gamma i[t - tau2] Exp[mu tau2]) + sigma1 s[t] dW1[t] + 
   lambda1 s[t] L1[t];
eq2 = -i'[t] + (beta s[t] d[t] - (mu2 + gamma) i[t]) + 
   sigma2 i[t] dW2[t] + lambda2 i[t] L2[t];
eq3 = -r'[t] + (gamma i[t] - mu3 r[t]);
eq4 = -d'[t] + eta (i[t] - d[t]) + sigma4 d[t] dW4[t] + 
   lambda4 d[t] L4[t];

rul[j_] := {a -> A[[j]], beta -> \[Beta][[j]], gamma -> \[Gamma][[j]],
    eta -> \[Eta][[j]], mu1 -> µ1[[j]], mu2 -> µ2[[j]], 
   sigma1 -> \[Sigma]1[[j]], sigma2 -> \[Sigma]2[[j]], 
   sigma4 -> \[Sigma]4[[j]], lambda1 -> \[Lambda]1[[j]], 
   lambda2 -> \[Lambda]2[[j]], lambda4 -> \[Lambda]4[[j]], 
   kappa -> \[Kappa][[j]], tau1 -> \[Tau]1[[j]], tau2 -> \[Tau]2[[j]],
    mu -> mu0[[j]]};

eqn[j_] := {eq1, eq2, eq4} /. rul[j];



ic1 = {s[t /; t <= 0] == 0.6, i[t /; t <= 0] == 0.3, 
  d[t /; t <= 0] == 0.05}; 
sol1[j_] := 
 NDSolve[{eqn[j] == {0, 0, 0}, ic1}, {s, i, d}, {t, 0, tmax - 1}];

With[{sol = sol1[1]}, {Plot[Evaluate[s[t] /. sol], {t, 0, tmax - 1}, 
   PlotRange -> All, Frame -> True, PlotStyle -> Green, 
   FrameLabel -> {"Time t (Days)", "S"}], 
  Plot[Evaluate[i[t] /. sol], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> Blue, 
   FrameLabel -> {"Time t (Days)", "I"}], 
  Plot[Evaluate[d[t] /. sol], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> Red, 
   FrameLabel -> {"Time t (Days)", "D"}]}]

Using input data from the paper, we have the following

SeedRandom[123];
A (*Recruitment rate*)= {1, 0.3, 0.6, 0.6};
mu0 (*Natural mortality rate of S*)= {0.09, 0.3, 0.4, 0.4};
\[Beta] (*Transmission rate*)= {0.18, 1.3, 0.35, 0.8};
\[Gamma] (*Recovered rate*)= {0.55, 0.05, 0.2, 0.3};
\[Delta] (*Recovered rate*)= {0.44, 0.05, 0.2, 0.3};
\[Alpha]1 (*General mortality of I*)= {0.4, 0.5, 0.3, 0.3};
\[Alpha]2 (*General mortality of I*)= {0.02, 0.5, 0.3, 0.3};
\[Epsilon] (*Exponentially fading memory rate*)= {0.6, 0.09, 0.7, 0.2};
\[Rho] (*Exponentially fading memory rate*)= {0.2, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= {0.2, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= {0.85, 0.25, 0.15, 0.15};
\[Sigma]3 (*Intensity of W3(t)*)= {0.5, 0.27, 0.13, 
  0.13}; \[Sigma]4 (*Intensity of W4(t)*)= {0.5, 0.27, 0.13, 0.13};
\[Lambda]1 (*Jump intensity of S*)= {0., 0.2, 0.5, 0.5};
\[Lambda]2 (*Jump intensity of I*)= {0., 0.23, 0.3, 0.3};
\[Lambda]3 (*Jump intensity of D*)= {0., 0.1, 0.7, 
  0.7}; \[Lambda]4(*Jump intensity of R*)= {0., 0.1, 0.7, 0.7};
\[Tau]1 (*Jump intensity of S*)= {0.5, 0.2, 0.5, 0.5};
\[Tau]2 (*Jump intensity of I*)= {0.5, 0.23, 0.3, 0.3};
\[Tau]3 (*Jump intensity of D*)= {0.5, 0.1, 0.7, 0.7};

tmax = 101; pWe1 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW1 = 
 Interpolation[Table[{(j - 1), pWe1[[j]]}, {j, Length[pWe1]}], 
  InterpolationOrder -> 1]; pWe2 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW2 = 
 Interpolation[Table[{(j - 1), pWe2[[j]]}, {j, Length[pWe2]}], 
  InterpolationOrder -> 1];
pWe3 = RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 
   1]]; dW3 = 
 Interpolation[Table[{(j - 1), pWe3[[j]]}, {j, Length[pWe3]}], 
  InterpolationOrder -> 1]; pWe4 = 
 RandomFunction[WhiteNoiseProcess[], {0, tmax}][[2, 1, 1]]; dW4 = 
 Interpolation[Table[{(j - 1), pWe4[[j]]}, {j, Length[pWe4]}], 
  InterpolationOrder -> 1];

pL1 = RandomFunction[PoissonProcess[1.], {0, tmax}];
pL2 = RandomFunction[PoissonProcess[1.1], {0, tmax}];
pL3 = RandomFunction[PoissonProcess[1.1], {0, tmax}];
pL4 = RandomFunction[PoissonProcess[.9], {0, tmax}];
ListStepPlot[{pL1, pL2, pL3, pL4}];

dpL1 = pL1["SliceData", Range[tmax]] // First // Differences;
dpL2 = pL2["SliceData", Range[tmax]] // First // Differences;
dpL3 = pL3["SliceData", Range[tmax]] // First // Differences;
dpL4 = pL4["SliceData", Range[tmax]] // First // Differences;

L1[t_] := If[t < 1, 0, dpL1[[Round[t]]]]/tmax;
L2[t_] := If[t < 1, 0, dpL2[[Round[t]]]]/tmax;
L3[t_] := If[t < 1, 0, dpL3[[Round[t]]]]/tmax;
L4[t_] := If[t < 1, 0, dpL4[[Round[t]]]]/tmax;

eq1 = -s'[t] + (a - (mu + rho) s[t] - beta s[t] i[t] + 
     rho s[t - tau1] Exp[-mu tau1] + gamma i[t - tau2] Exp[-mu tau2] +
      epsilon q[t - tau3] Exp[-mu tau3]) + sigma1 s[t] dW1[t] + 
   lambda1 s[t] L1[t];
eq2 = -i'[t] + (beta s[t] i[t] - (mu + alpha1 + delta + gamma) i[t]) +
    sigma2 i[t] dW2[t] + lambda2 i[t] L2[t];
eq3 = -q'[t] + (delta i[t] - (mu + alpha2 + epsilon) q[t]) + 
  sigma3 q[t] dW3[t] + 
  lambda3 q[t] L3[t]; eq4 = -r'[t] + (gamma i[t] + rho s[t] + 
    epsilon q[t] - mu r[t] - rho s[t - tau1] Exp[-mu tau1] - 
    gamma i[t - tau2] Exp[-mu tau2] - 
    epsilon q[t - tau3] Exp[-mu tau3]) + sigma4 d[t] dW4[t] + 
  lambda4 d[t] L4[t];

ic1 = {s[t /; t <= 0] == 5, i[t /; t <= 0] == 15, q[t /; t <= 0] == 5};
rul[j_] := {a -> A[[j]], beta -> \[Beta][[j]], rho -> \[Rho][[j]], 
   gamma -> \[Gamma][[j]], epsilon -> \[Epsilon][[j]], 
   delta -> \[Delta][[j]], alpha1 -> \[Alpha]1[[j]], 
   alpha2 -> \[Alpha]2[[j]], sigma1 -> \[Sigma]1[[j]], 
   sigma2 -> \[Sigma]2[[j]], sigma3 -> \[Sigma]3[[j]], 
   sigma4 -> \[Sigma]4[[j]], lambda1 -> \[Lambda]1[[j]], 
   lambda2 -> \[Lambda]2[[j]], lambda3 -> \[Lambda]3[[j]], 
   lambda4 -> \[Lambda]4[[j]], tau1 -> \[Tau]1[[j]], 
   tau2 -> \[Tau]2[[j]], tau3 -> \[Tau]3[[j]], mu -> mu0[[j]]};

eqn[j_] := {eq1, eq2, eq3} /. rul[j];


sol1[j_] := 
  NDSolve[{eqn[j] == {0, 0, 0}, ic1}, {s, i, q}, {t, 0, tmax - 1}, 
   MaxSteps -> 10^6];


With[{sol = sol1[1]}, {Plot[Evaluate[s[t] /. {sol}], {t, 0, tmax - 1},
    PlotRange -> All, Frame -> True, PlotStyle -> {Green}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["S(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18], 
  Plot[Evaluate[i[t] /. {sol}], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> {Blue}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["I(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18], 
  Plot[Evaluate[q[t] /. {sol}], {t, 0, tmax - 1}, PlotRange -> All, 
   Frame -> True, PlotStyle -> {Red}, 
   FrameLabel -> {Style["Time(Days)", 20, Black], 
     Style["Q(t)", 20, Black]}, ImageSize -> 500, 
   FrameTicksStyle -> 18]}]

Comments are not for extended discussion; this conversation has been moved to chat. — Kuba, Commented Jan 30, 2023 at 7:28

Stack Exchange Network

Extension of: Numerically solving a system of SDE's with Levy noise?

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged
differential-equations
modeling
stochastic-calculus
system-modeling
or ask your own question.

Linked

Hot Network Questions

Extension of: Numerically solving a system of SDE's with Levy noise?

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged differential-equationsmodelingstochastic-calculussystem-modeling or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
differential-equations
modeling
stochastic-calculus
system-modeling
or ask your own question.