# Find the missing initial conditions for system of ODES

I want to solve the following system of ODEs:

   s1 (t) = -(k1*s1 *(k2*s2 + s3))/m,
s2 (t) = ((k1*s1 *(k2*s2 + s3))/m) - k3*s2,
s3' (t) = k3*s2 - k4*s3
s4' (t) = k4*s3


where k1=0.4029;k2=0.7;k3=0.41;k4=0.182; m=s1+s2+s3+s4;

I used:

    k1 = 0.4029; k2 = 0.78; k3 = 0.41; k4 = 0.182; m = s1+s2+s3+s4;

u3[s1, s2, s3,
s4] = {s1, s2, s3, s4} /.
First@DSolve[{s1'[t] == -(k1*s1[t]*(k2*s2[t] + s3[t]))/m,
s2'[t] == ((k1*s1[t]*(k2*s2[t] + s3[t]))/m) - k3*s2[t],
s3'[t] == k3*s2[t] - k4*s3[t], s4'[t] == k4*s3[t]}, {s1, s2, s3,
s4}, t] /. {C[1] -> c1, C[2] -> c2, C[3] -> c3, C[4] -> c4}


But it is not working, any help?

1) I need numerical solution for the system.

2) If possible can we get analytic solution or approximate solutions for each unknown functions si, i=1..4, given in a polynomials? Thanks

• You try to reinvent a SEIR-model? The initial conditions aren't unique, you need more information . – Ulrich Neumann Mar 27 '20 at 10:56
• Dear Ulrich, what type of information should I have? – user62716 Mar 27 '20 at 11:01
• The initial conditions! Or if you know the solution (measurement) you could try to find the optimal initial conditions . – Ulrich Neumann Mar 27 '20 at 11:18
• If Mathematica could solve your ode without i.c. the solution would depend on several free parameters! These parameters usually are adjusted by additional conditions. – Ulrich Neumann Mar 27 '20 at 11:31
• @user62716 Then just make up some plausible initial conditions. Maybe s1[0]==n-0.001, s2[0]==0.001, s3[0]==s4[0]==0 is interesting. BTW, don't use N, it's a built-in function. – Chris K Mar 27 '20 at 19:48

## 2 Answers

Let me show how to get an asymptotic polynomial solution at least for small t.

k1 = 0.4029; k2 = 0.78; k3 = 0.41; k4 = 0.182;
m = s1[t] + s2[t] + s3[t] + s4[t];

eqs = {s1'[t] == -(k1*s1[t]*(k2*s2[t] + s3[t]))/m,
s2'[t] == ((k1*s1[t]*(k2*s2[t] + s3[t]))/m) - k3*s2[t],
s3'[t] == k3*s2[t] - k4*s3[t], s4'[t] == k4*s3[t],
s1[0] == 3.9, s2[0] == 8.5, s3[0] == 3.7, s4[0] == 0.8} //
Simplify


Now define the si as general polynoms, here up to order 6 and then solve for all the prefactors of powers of t.

{s1[t_] = Sum[a[1][s] t^s, {s, 0, 6}] + O[t]^7,
s2[t_] = Sum[a[2][s] t^s, {s, 0, 6}] + O[t]^7,
s3[t_] = Sum[a[3][s] t^s, {s, 0, 6}] + O[t]^7,
s4[t_] = Sum[a[4][s] t^s, {s, 0, 6}] + O[t]^7};

le = LogicalExpand[#] & /@ eqs;

sol1 = NSolve[And @@ le,
Flatten[Table[a[j][s], {s, 0, 6}, {j, 1, 4}], 1]];

{ss1[t_], ss2[t_], ss3[t_], ss4[t_]} =
Normal[{s1[t], s2[t], s3[t], s4[t]} /. First@sol1] // Simplify;

{ss1[t], ss2[t], ss3[t], ss4[t]} // TableForm

(*   {3.9\[VeryThinSpace]- 0.960452 t + 0.0791007 t^2 + 0.013307 t^3 - 0.00484174 t^4 + 0.000626304 t^5 - 0.0000156141 t^6,

8.5\[VeryThinSpace]- 2.52455 t + 0.438432 t^2 - 0.073226 t^3 + 0.0123474 t^4 - 0.00163879 t^5 + 0.000127598 t^6,

3.7\[VeryThinSpace]+ 2.8116 t - 0.773388 t^2 + 0.106838 t^3 - 0.0123668 t^4 + 0.00146264 t^5 - 0.000156351 t^6,
0.8\[VeryThinSpace]+ 0.6734 t + 0.255856 t^2 - 0.0469189 t^3 + 0.00486112 t^4 - 0.000450151 t^5 + 0.0000443667 t^6}   *)


For comparison., i take the point data @Nasser provided. It shows, this polynomial aproximation is good for small t. Solving for higher order of polynoms however quickly needs very long calculation time of NSolve. Solving for arbitrary parameters with Solve to higher orders is nearly not practicable.

(*   tp=Transpose[  point data of @Nasser];   *)

tp2 = Transpose[{Range[0, 20, 1], #}] & /@ tp;

lp = ListPlot[tp2, PlotStyle -> {Black, Cyan, Red, Blue},
PlotRange -> {{0, 6}, {0, 9}}];

pl = Plot[Evaluate[{ss1[t], ss2[t], ss3[t], ss4[t]}], {t, 0, 6},
PlotStyle -> {Black, Cyan, Red, Blue}];

Show[pl, lp]


Edit How to determine any kind of error

Since in most cases you do not know the exact solution, you can test, how close to zero is left minus right side of all equations for all t. ( lhs-rhs is zero with exact solution)

First clear the polynom definitions of si in order to go back to the original equation and then substitute the found solutions i called ssi.

{ss1[t_], ss2[t_], ss3[t_], ss4[t_]} =
Normal[{s1[t], s2[t], s3[t], s4[t]} /. First@sol5] // Simplify

Clear[s1, s2, s3, s4]

plorder5 =
Plot[Evaluate[
eqs[[All, 1]] - eqs[[All, 2]] /. {s1 -> ss1, s2 -> ss2, s3 -> ss3,
s4 -> ss4}], {t, 0, 6}, PlotRange -> .2]


At order 7 error is smaller than with order 5. But be aware, this is not the error of found solutions (which is normally smaller as you see), but the error of the differential equations.

{ss1[t_], ss2[t_], ss3[t_], ss4[t_]} =
Normal[{s1[t], s2[t], s3[t], s4[t]} /. First@sol7] // Simplify

Clear[s1, s2, s3, s4]

plorder7 =
Plot[Evaluate[
eqs[[All, 1]] - eqs[[All, 2]] /. {s1 -> ss1, s2 -> ss2, s3 -> ss3,
s4 -> ss4}], {t, 0, 6}, PlotRange -> .2]


• Dear Akku14, many thanks for the above approximate solution...Best regards – user62716 Mar 29 '20 at 13:02
• Dear Akku14, can we evaluate the error (say relative error, or any type of error) for different order say 2,3,4..6 and plot it to see if the errors decrease when we increase the order?of course with fixed Nasser results(we consider Nasser results as the exact solution)?..Best regards. – user62716 Mar 29 '20 at 16:21
• Dear Akku14, the above code should be add at the end of previous code? or shall we re-arrange it order 5 and 7? as the original code at order 6? – user62716 Mar 29 '20 at 19:35
• You took t=0 to 6?for both order 5 and 7? is it correct? – user62716 Mar 29 '20 at 19:37

Can you send me the NDSolve line code

Updated for the new constraint just added to the question which was not there before, which is m=s1+s2+s3+s4. Where m here is total population.

I tried to use DAE to add this to the system of ode's. But it turned out this is not needed.

All what is needed is to set m to be the the sum of s1+s2+s3+s4 at time zero, then this constraint is automatically satisfied for all time.

So modified the code a little to remove m as variable using slider, and set it internally.

Manipulate[
Module[{ode1, ode2, ode3, ode4, s1, s2, s3, s4, t, ic, sys, m,data},

m = s10 + s20 + s30 + s40;
ode1 = s1'[t] == -(k1*s1[t]*(k2*s2[t] + s3[t]))/m;
ode2 = s2'[t] == ((k1*s1[t]*(k2*s2[t] + s3[t]))/m) - k3*s2[t];
ode3 = s3'[t] == k3*s2[t] - k4*s3[t];
ode4 = s4'[t] == k4*s3[t];

sys = {ode1, ode2, ode3, ode4};

ic = {s1[0] == s10, s2[0] == s20, s3[0] == s30, s4[0] == s40};
{s1, s2, s3, s4} =
Quiet@NDSolveValue[{ode1, ode2, ode3, ode4, ic}, {s1, s2, s3,
s4}, {t, 0, tMax}];

data = Table[{s1[i], s2[i], s3[i], s4[i]}, {i, 0, tMax, 1}];
Print[data];

Plot[
{s1[t], s2[t], s3[t], s4[t], m[t]}, {t, 0, tMax},
PlotLegends -> {"s1  susceptible",
"s2 exposed but not yet infectious", "s3 Infected",
"s4 Recovered"},
PlotRange -> {Automatic, {0, 20}}, GridLines -> Automatic,
GridLinesStyle -> LightGray,
ImageSize -> 500, ImagePadding -> 20,
PlotLabel -> Style[Row[{"SEIR infection model"}], Bold, 14],
PlotStyle -> {Black, Cyan, Red, Blue}
]

],
{{s10, 3.9, "susceptible s1(0)"}, 0, 10, .1, Appearance -> "Labeled"},
{{s20, 8.5, "exposed but not yet infectious s2(0)"}, 0, 10, .1,
Appearance -> "Labeled"},
{{s30, 3.7, "Infected s3(0)"}, 0, 10, .1, Appearance -> "Labeled"},
{{s40, 0.8, "Recovered s4(0)"}, 0, 10, .1, Appearance -> "Labeled"},
{{k1, 0.4029, "k1"}, 0, 1, .0001, Appearance -> "Labeled"},
{{k2, 0.78, "k2"}, 0, 1, .01, Appearance -> "Labeled"},
{{k3, 0.41, "k3"}, 0, 1, .01, Appearance -> "Labeled"},
{{k4, 0.182, "k4"}, 0, 1, .001, Appearance -> "Labeled"},
{{tMax, 20, "time?"}, 1, 100, 0.1, Appearance -> "Labeled"},
TrackedSymbols :> {s10, s20, s30, s40, k1, k2, k3, k4, tMax}
]


Prints s1,s2,s3,s4 for the time used at intervals selected in code above

{{3.9, 8.5, 3.7, 0.8}, {3.02772, 6.35149, 5.834, 1.68679}, {2.34279,
4.77191, 6.92468, 2.86063}, {1.82802, 3.58489, 7.3216,
4.16549}, {1.44877, 2.68696, 7.2659, 5.49837}, {1.17086, 2.00878,
6.92715, 6.79321}, {0.966539, 1.49902, 6.42446, 8.00998}, {0.815065,
1.11785, 5.84026, 9.12682}, {0.701542, 0.834097, 5.23006,
10.1343}, {0.615442, 0.623534, 4.62983, 11.0312}, {0.549355,
0.467548, 4.06158, 11.8215}, {0.498045, 0.352027, 3.53761,
12.5123}, {0.457782, 0.266386, 3.06359, 13.1122}, {0.425883,
0.202753, 2.64079, 13.6306}, {0.400392, 0.155313, 2.26778,
14.0765}, {0.379865, 0.119791, 1.94147, 14.4589}, {0.363226,
0.0930537, 1.65792, 14.7858}, {0.349659, 0.0728083, 1.41284,
15.0647}, {0.338543, 0.0573777, 1.20193, 15.3022}, {0.329394,
0.0455346, 1.02105, 15.504}, {0.321838, 0.0363782, 0.86638,
15.6754}}

• Great expert, please can we get approximate solution, not numerical, i.e. for example s1= function,...etc – user62716 Mar 28 '20 at 20:34
• Dear Nasser, also m in fact equal to s1+s2+s3+s4 – user62716 Mar 28 '20 at 20:38
• @user62716 please add this to your question as it was not there. This makes the system a DAE and I need to update my answer in this case. – Nasser Mar 28 '20 at 20:54
• Many thanks for helping me Dr.Nasser, I have updated my questions – user62716 Mar 28 '20 at 21:02
• @user62716 updated for m=s1+s2+s3+s4 For your question: DSolve can not solve this analytically, so only numerical solution is possible. – Nasser Mar 28 '20 at 21:40