1
$\begingroup$

Following this solution by Alex: Numerically solving a system of SDE's with Levy noise?, I tried Alex's solution for a different model:

b(*Recruitment rate*)= {0.8, 0.3, 0.6, 0.6};
µ(*Natural mortality rate of S,I,*)= {0.1, 0.3, 0.4, 0.4};
\[Beta]a(*Transmission rate*)= {0.7, 1.3, 0.35, 0.8};
\[Beta]i(*Transmission rate*)= {0.7, 1.3, 0.35, 0.8};
\[Rho] (*- rate*)= {0.2, 1.3, 0.35, 0.8};
\[Xi](*- rate*)= {0.1, 1.3, 0.35, 0.8};
\[Alpha](*DID rate*)= {0.1, 1.3, 0.35, 0.8};
\[Gamma]a(*Recovered rate*)= {0.2, 0.05, 0.2, 0.3};
\[Gamma]i(*Recovered rate*)= {0.2, 0.05, 0.2, 0.3};
\[Nu](*vac rate*)= {0.4, 0.05, 0.2, 0.3};
\[Theta] (*transmission rate of I*)= {0.3, 0.5, 0.3, 0.3};
\[Eta] (*Treatment rate*)= {0.5, 0.09, 0.7, 0.2};
\[Sigma]1 (*Intensity of W1(t)*)= {0.1, 0.15, 0.2, 0.169};
\[Sigma]2 (*Intensity of W2(t)*)= {0.2, 0.25, 0.15, 0.15};
\[Sigma]3 (*Intensity of W3(t)*)= {0.15, 0.25, 0.15, 0.15};
\[Sigma]4 (*Intensity of W4(t)*)= {0.2, 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 A*)= {0.23, 0.23, 0.3, 0.3};
\[Lambda]3 (*Jump intensity of I*)= {0.15, 0.23, 0.3, 0.3};
\[Lambda]4 (*Jump intensity of R*)= {0.1, 0.1, 0.7, 0.7};

tmax = 1001; 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], {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] + (b (1 - nu) - (betaa a[t]/(s[t] + a[t] + i[t] + r[t]) + 
        betai i[t]/(s[t] + a[t] + i[t] + r[t])) s[t] - (mu + rho) s[
       t] + xi r[t]) + sigma1 s[t] dW1[t] + lambda1 s[t] L1[t];
eq2 = -a'[
     t] + ((betaa a[t]/(s[t] + a[t] + i[t] + r[t]) + 
        betai i[t]/(s[t] + a[t] + i[t] + r[t])) s[
       t] - (gammaa + theta + mu) a[t]) + sigma2 a[t] dW2[t] + 
   lambda2 a[t] L2[t];
eq3 = -i'[t] + (theta a[t] - (gammai + eta + alpha + mu) i[t]) + 
   sigma3 i[t] dW3[t] + lambda3 i[t] L2[t];
eq4 = -r'[t] + (b nu + rho s[t] + 
     gammaa a[t] + (gammai + eta) i[t] - (xi + mu) r[t]) + 
   sigma4 r[t] dW4[t] + lambda4 r[t] L4[t];

ic = {s[0] == 0.6, a[0] == 0.3, i[0] == 0.1, r[0] == 0.05};
rul[j_] := {b -> b[[j]], nu -> \[Nu][[j]], betaa -> \[Beta]a[[j]], 
   betai -> \[Beta]i[[j]], gammaa -> \[Gamma]a[[j]], 
   gammai -> \[Gamma]i[[j]], eta -> \[Eta][[j]], mu -> µ[[j]], 
   rho -> \[Rho][[j]], xi -> \[Xi][[j]], theta -> \[Theta][[j]], 
   alpha -> \[Alpha][[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]]};

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

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

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

Ignore the parameter definitions. This gives me errors for some reason even though it is almost identical to Alex's answer.. any idea why?

$\endgroup$
3
  • $\begingroup$ A minimal working example of your problem would be easier to answer. $\endgroup$
    – rhermans
    Jan 19, 2023 at 18:55
  • $\begingroup$ @rhermans I dont know how to simplify it.. $\endgroup$
    – Math
    Jan 19, 2023 at 19:05
  • $\begingroup$ You are not even trying. You don't need a list of 4 Plot to see the error. sol[1] fails but eqn[1] doesn't. $\endgroup$
    – rhermans
    Jan 19, 2023 at 19:09

1 Answer 1

2
$\begingroup$

There are typos in the code with b definition. Introducing b0 in equations and rul we have

SeedRandom[1234];

b(*Recruitment rate*)= {0.8, 0.3, 0.6, 0.6};
µ(*Natural mortality rate of S,I,*)= {0.1, 0.3, 0.4, 0.4};
βa(*Transmission rate*)= {0.7, 1.3, 0.35, 0.8};
βi(*Transmission rate*)= {0.7, 1.3, 0.35, 0.8};
ρ (*-rate*)= {0.2, 1.3, 0.35, 0.8};
ξ(*-rate*)= {0.1, 1.3, 0.35, 0.8};
α(*DID rate*)= {0.1, 1.3, 0.35, 0.8};
γa(*Recovered rate*)= {0.2, 0.05, 0.2, 0.3};
γi(*Recovered rate*)= {0.2, 0.05, 0.2, 0.3};
ν(*vac rate*)= {0.4, 0.05, 0.2, 0.3};
θ (*transmission rate of I*)= {0.3, 0.5, 0.3, 0.3};
η (*Treatment rate*)= {0.5, 0.09, 0.7, 0.2};
σ1 (*Intensity of W1(t)*)= {0.1, 0.15, 0.2, 0.169};
σ2 (*Intensity of W2(t)*)= {0.2, 0.25, 0.15, 0.15};
σ3 (*Intensity of W3(t)*)= {0.15, 0.25, 0.15, 0.15};
σ4 (*Intensity of W4(t)*)= {0.2, 0.27, 0.13, 0.13};
λ1 (*Jump intensity of S*)= {0.2, 0.2, 0.5, 0.5};
λ2 (*Jump intensity of A*)= {0.23, 0.23, 0.3, 0.3};
λ3 (*Jump intensity of I*)= {0.15, 0.23, 0.3, 0.3};
λ4 (*Jump intensity of R*)= {0.1, 0.1, 0.7, 0.7};

tmax = 1001; 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], {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] + (b0 (1 - nu) - (betaa a[t]/(s[t] + a[t] + i[t] + r[t]) + 
        betai i[t]/(s[t] + a[t] + i[t] + r[t])) s[t] - (mu + rho) s[
       t] + xi r[t]) + sigma1 s[t] dW1[t] + lambda1 s[t] L1[t];
eq2 = -a'[
     t] + ((betaa a[t]/(s[t] + a[t] + i[t] + r[t]) + 
        betai i[t]/(s[t] + a[t] + i[t] + r[t])) s[
       t] - (gammaa + theta + mu) a[t]) + sigma2 a[t] dW2[t] + 
   lambda2 a[t] L2[t];
eq3 = -i'[t] + (theta a[t] - (gammai + eta + alpha + mu) i[t]) + 
   sigma3 i[t] dW3[t] + lambda3 i[t] L2[t];
eq4 = -r'[t] + (b0 nu + rho s[t] + 
     gammaa a[t] + (gammai + eta) i[t] - (xi + mu) r[t]) + 
   sigma4 r[t] dW4[t] + lambda4 r[t] L4[t];

ic = {s[0] == 0.6, a[0] == 0.3, i[0] == 0.1, r[0] == 0.05};
rul[j_] := {b0 -> b[[j]], nu -> ν[[j]], betaa -> βa[[j]], 
   betai -> βi[[j]], gammaa -> γa[[j]], 
   gammai -> γi[[j]], eta -> η[[j]], mu -> µ[[j]], 
   rho -> ρ[[j]], xi -> ξ[[j]], theta -> θ[[j]], 
   alpha -> α[[j]], sigma1 -> σ1[[j]], 
   sigma2 -> σ2[[j]], sigma3 -> σ3[[j]], 
   sigma4 -> σ4[[j]], lambda1 -> λ1[[j]], 
   lambda2 -> λ2[[j]], lambda3 -> λ3[[j]], 
   lambda4 -> λ4[[j]]};

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


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


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

Figure 1

$\endgroup$
2
  • $\begingroup$ You are amazing Sir! This is dearly appreciated! $\endgroup$
    – Math
    Jan 20, 2023 at 12:15
  • $\begingroup$ @Math You are welcome! $\endgroup$ Jan 20, 2023 at 12:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.