# Why are there problems solving this SDE system?

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?

• A minimal working example of your problem would be easier to answer. Commented Jan 19, 2023 at 18:55
• @rhermans I dont know how to simplify it..
– Math
Commented Jan 19, 2023 at 19:05
• 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. Commented Jan 19, 2023 at 19:09

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]}]


• You are amazing Sir! This is dearly appreciated!
– Math
Commented Jan 20, 2023 at 12:15
• @Math You are welcome! Commented Jan 20, 2023 at 12:56