Non Linear System of Diferential Equations with periodic function

Question

I am programming a system of differential equations consisting of 3 cell populations; Normal cells N, Tumor cells T and chemotherapy treatment Q, when a function q is constant. The system is as follows

s = NDSolve[{NN'[t] == \[Alpha] NN[t] (1 - NN[t]/k1) - (\[Gamma] NN[t] T[t])/(b + T[t]) - \[Mu]1 Q[t] NN[t], T'[t] == \[Beta] T[t] (1 - T[t]/k2) - \[Mu]2 Q[t] T[t], Q'[t] == q - \[Lambda] Q[t], NN[0] == 100, T[0] == 200, Q[0] == 300}, {NN[t], T[t], Q[t]}, {t, 1, 100}]

Now I need to program the same system but the function q is variable and periodic, that is.

g[t_] := q (1 - Exp[-(nn + 1) k t])/(1 - Exp[-k t])

Where nn is an integer.

However, it is not possible with:

s = NDSolve[{NN'[t] == \[Alpha] NN[t] (1 - NN[t]/k1) - (\[Gamma] NN[t] T[t])/(b + T[t]) - \[Mu]1 Q[t] NN[t], T'[t] == \[Beta] T[t] (1 - T[t]/k2) - \[Mu]2 Q[t] T[t], Q'[t] == g[t] - \[Lambda] Q[t], NN[0] == 100, T[0] == 200, Q[0] == 300}, {NN[t], T[t], Q[t]}, {t, 1, 100}]

Changing q to g[t], I get the error message: power::infy: "Infinite expression 1/0. encountered".

Can someone help me, how can I fix my code?

Thanks a lot

Answer 1

Your main problem was not to define all paramters.

A good way to look whether you got all variables is

Variables[Level[eqs, -1]]

Two posibilities of solution:

a) The third equation can be solved separately since it depends only on Q. (I set your constant q in q[t] to g and set k to 1/4 and nn to 3 for test).

b) Complete numerial solution

q[t_] = g (1 - Exp[-(nn + 1) k t])/(1 - Exp[-k t]);

(eqs = {NN'[
      t] == \[Alpha] NN[
        t] (1 - NN[t]/k1) - (\[Gamma] NN[t] T[t])/(b + 
         T[t]) - \[Mu]1 Q[t] NN[t], 
    T'[t] == \[Beta] T[t] (1 - T[t]/k2) - \[Mu]2 Q[t] T[t], 
    Q'[t] == q[t] - \[Lambda] Q[t], NN[0] == 100, T[0] == 200, 
    Q[0] == 300}) // TableForm

Qsol[g_, k_, nn_, \[Lambda]_] = 
 Q /. Flatten@DSolve[eqs[[{3, 6}]], Q, t]

(*   Function[{t}, (1/(k*(k*nn - \[Lambda])*(k + \[Lambda])))*
   ((300*k^3*nn - 300*k^2*\[Lambda] + 300*k^2*nn*\[Lambda] - 300*k*\[Lambda]^2 - 
     E^(t*(k + \[Lambda]))*g*k^2*nn*Hypergeometric2F1[1, (k + \[Lambda])/k, 
       2 + \[Lambda]/k, E^(k*t)] + E^(t*(k + \[Lambda]))*g*k*\[Lambda]*
      Hypergeometric2F1[1, (k + \[Lambda])/k, 2 + \[Lambda]/k, E^(k*t)] - 
     E^(t*((-k)*nn + \[Lambda]))*g*k^2*Hypergeometric2F1[1, -nn + \[Lambda]/k, 
       1 - nn + \[Lambda]/k, E^(k*t)] - E^(t*((-k)*nn + \[Lambda]))*g*k*\[Lambda]*
      Hypergeometric2F1[1, -nn + \[Lambda]/k, 1 - nn + \[Lambda]/k, E^(k*t)] - 
     g*k^2*nn*PolyGamma[0, (k + \[Lambda])/k] + 
     g*k*\[Lambda]*PolyGamma[0, (k + \[Lambda])/k] - 
     g*k*nn*\[Lambda]*PolyGamma[0, (k + \[Lambda])/k] + 
     g*\[Lambda]^2*PolyGamma[0, (k + \[Lambda])/k] + 
     g*k^2*nn*PolyGamma[0, -nn + \[Lambda]/k] - 
     g*k*\[Lambda]*PolyGamma[0, -nn + \[Lambda]/k] + 
     g*k*nn*\[Lambda]*PolyGamma[0, -nn + \[Lambda]/k] - 
     g*\[Lambda]^2*PolyGamma[0, -nn + \[Lambda]/k])/E^(t*\[Lambda]))]   *)

Limit[Qsol[g, k, nn, \[Lambda]][t], t -> 0, 
 Assumptions -> 
  Element[nn, Integers] && Element[{g, k, \[Lambda]}, Reals]]

(*   300   *)

Qs[g_, k_, nn_, \[Lambda]_] = 
 Function[t, Piecewise[{{300, t == 0}}, Qsol[g, k, nn, \[Lambda]][t]]]

eqs2[b_, g_, k_, k1_, k2_, 
  nn_, \[Alpha]_, \[Beta]_, \[Gamma]_, \[Lambda]_, \[Mu]1_, \[Mu]2_] =
  eqs[[{1, 2, 4, 5}]] /. Q -> Qs[g, k, nn, \[Lambda]]

pars = {\[Alpha] -> 13.9, \[Beta] -> 37, k -> 1/4, k1 -> 1460, 
    k2 -> 167, nn -> 3, \[Gamma] -> 0.0075, 
    b -> 100, \[Mu]1 -> 0.006, \[Mu]2 -> 0.0008, \[Lambda] -> 0.9, 
    g -> 100} // Rationalize[#, 0] & // Sort

ndsol = NDSolve[
  Evaluate[eqs2[b, g, k, k1, k2, 
     nn, \[Alpha], \[Beta], \[Gamma], \[Lambda], \[Mu]1, \[Mu]2] /. 
    pars], {NN, T}, {t, 0, 100}]

Plot[Evaluate[NN[t] /. ndsol[[1]]], {t, 0, 30}, PlotRange -> All]

Plot[Evaluate[T[t] /. ndsol[[1]]], {t, 0, 30}, PlotRange -> All]

Plot[Evaluate[Qs[g, k, nn, \[Lambda]][t] /. pars], {t, 0, 30}, 
 PlotRange -> All]

b) complete numerical solution

ndsol2 = NDSolve[Evaluate[eqs /. pars], {NN, T, Q}, {t, 0, 100}]

Plot[Evaluate[NN[t] /. ndsol2[[1]]], {t, 0, 30}, PlotRange -> All]