How to force my differential equations give a bounded solution?

Question

I have modeled the interaction of two physical quantities, $S$ and $B$, by the following differential equations (the second one is a delay differential equation):

$$S'(t) = 0.31 S(t) \Big( 1 - \frac{S(t)}{6.19} \Big)- \frac{a}{1 + B(t)} S(t),$$

$$B'(t) = c (1 - B(t)) \theta(18 - t) - d B(t) S(t) - 3 c (1 - B(t - 18)),$$

where $a$, $c$, and $d$ are some constants to be determined from the experimental data and $\theta$ is the Heaviside function. The parameters are found to be: $a = 0.72$, $c = 0.19$, and $d = 0.24$; and also, the initial conditions read: $S(0) = 1.25$ and $B(0) = 1.14$.

The solutions (plotted with Mathematica) are:

Unprotect[HeavisideTheta]; HeavisideTheta[0.] = 0.5; Protect[HeavisideTheta];
so = NDSolve[{ss'[t] + 0.72/(1 + bb[t]) ss[t] - 0.31 ss[t] (1 - ss[t]/6.19) == 0, 
bb'[t] - 0.19 (1 - bb[t]) HeavisideTheta[18 - t] + 0.24 bb[t] ss[t] + 0.57 (1 - bb[t - 18]) == 0, 
ss[t /; t < 0] == 1.25, bb[t /; t < 0] == 1.14}, {ss, bb}, {t, 0, 24}];

Plot[Evaluate[bb[t] /. so], {t, 0, 24}, PlotRange -> All, AxesOrigin -> {0, 0}, PlotLabel -> "B(t)"]
Plot[Evaluate[ss[t] /. so], {t, 0, 24}, PlotRange -> All, AxesOrigin -> {0, 0}, PlotLabel -> "S(t)"]

Image 1

Everything looks good, since I'm only interested in the time interval $0 < t < 24$ and the curves correctly predict the behavior of my system (experiment). However, if we simulate the differential equations further and for a larger time interval, let's say, till $t = 30$, the solutions look like:

Image 2

This is not desirable as $B$ is a physical quantity and cannot be negative.

My question is:

What kinds of terms should I add in my differential equations (or what assumptions need to be considered?) to force the $B$ solution do not cross the axis and to be bounded and remain zero at the end?

Any help is appreciated.

Answer 1

Changing the wheighting factor to k^alpha as can be depicted in the function fitness we can obtain a solution to the odes with the positive requirement.

fitness[a00_?NumberQ, c00_?NumberQ, d00_?NumberQ, s0_?NumberQ, b0_?NumberQ] := Module[{ss, bb, parms, odes, ODES, solode, cinits, ssbb, t, alpha = 1.1}, 
  odes = {ss'[t] + a00/(1 + bb[t]) ss[t] - 0.31 ss[t] (1 - ss[t]/6.19) == 0, bb'[t] - c00 (1 - bb[t]) + d00 bb[t] ss[t] + 3 c00 (1 - bb[t - 18]) == 0};
  cinits = {ss[t /; t <= 0] == s0, bb[t /; t <= 0] == b0};
  ODES = Join[odes, cinits];
  solode = Quiet@NDSolve[ODES, {ss, bb}, {t, tk[[1]], tk[[8]]}][[1]];
  Return[Sum[Norm[{y1[[k]], k^alpha y2[[k]]} - 
  Evaluate[{ss[t], k^alpha bb[t]} /. solode /. {t -> tk[[k]]}]], {k, 1, 8}]]]



tk = {0, 3, 6, 9, 15, 18, 21, 24};
y1 = {0.974, 0.492, 0.829, 0.554, 0.208, 0.138, 0.199, 0.0893};
y2 = {0.915094, 0.736097, 0.793694, 0.833664, 1, 0.99578, 0.897964, 0.214499};
SdatawithSD = {{0, Around[0.974, 0.289]}, {3, Around[0.492, 0.165]}, {6, Around[0.829, 0.404]}, {9, Around[0.554, 0.245]}, {15, Around[0.208, 0.191]}, {18, Around[0.138, 0.0962]}, {21, Around[0.199, 0.241]}, {24, Around[0.0893, 0.0359]}};
BdatawithSD = {{0, Around[0.915094, 0.1]}, {3, Around[0.736097, 0.091668]}, {6, Around[0.793694, 0.082575]}, {9, Around[0.833664, 0.070242]}, {15, Around[1, 0.002851]}, {18, Around[0.99578, 0.015591]}, {21,Around[0.897964, 0.04783]}, {24, Around[0.214499, 0.01231]}};


restrs = {0.7 - 0.2 < a00 < 1 + 0.2 && 0.1 < c00 < 0.5 + 0.2 && 0.1 < d00 < 0.5 + 0.2 && 0.9 - 0.2 < s0 < 1.2 + 0.2 && 0.9 - 0.2 < b0 < 1.2 + 0.2};
vars = {a00, c00, d00, s0, b0};

sol = Quiet@NMinimize[Join[{fitness[a00, c00, d00, s0, b0]}, restrs], vars, StepMonitor :> Print[{a00, c00, d00, s0, b0, fitness[a00, c00, d00, s0, b0]}], MaxIterations -> 10, WorkingPrecision -> 30, Method -> "DifferentialEvolution"]

odes = {ss'[t] + a00/(1 + bb[t]) ss[t] - 0.31 ss[t] (1 - ss[t]/6.19) == 0, bb'[t] - c00 (1 - bb[t]) + d00 bb[t] ss[t] + 3 c00 (1 - bb[t - 18]) == 0} /. sol[[2]];
cinits = {ss[t /; t <= 0] == s0, bb[t /; t <= 0] == b0} /. sol[[2]];
ODES = Join[odes, cinits];
solode = NDSolve[ODES, {ss, bb}, {t, tk[[1]], 1.5 tk[[8]]}][[1]];

g11 = ListPlot[SdatawithSD, PlotStyle -> {Orange, PointSize[0.02]}, AxesOrigin -> {-1, 0}];
gr10 = ListPlot[Transpose[{tk, y1}], PlotStyle -> Red, PlotLabel -> "S(t)"];
gr0 = Plot[Evaluate[ss[t] /. solode], {t, tk[[1]], 1.5 tk[[8]]}];
Show[g11, gr10, gr0, PlotRange -> All]

g22 = ListPlot[BdatawithSD, AxesOrigin -> {-1, 0}, PlotStyle -> {PointSize[0.02]}];
gr2 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue];
gr00 = Plot[Evaluate[bb[t] /. solode], {t, tk[[1]], 1.5 tk[[8]]}];
gr02 = ListPlot[Transpose[{tk, y2}], PlotStyle -> Blue, PlotLabel -> "B(t)"];
Show[g22, gr02, gr00, PlotRange -> All]