# Non Linear System of Diferential Equations with periodic function

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

• Please provide the values of all parameters, so that readers can run your code. Apr 25 at 22:58
• Parameters are: [Alpha] = 13.9; [Beta] = 37; k1 = 1460; k2 = 167; [Gamma] = 0.0075; b = 100; [Mu]1 = 0.006; [Mu]2 = 0.0008; [Lambda] = 0.9; q = 100; Apr 25 at 23:44
• Your problem is that (1 - Exp[-k t]) is zero at t = 0. By the way, nn remains undefined. Apr 26 at 0:36
• k also is undefined. Apr 26 at 0:47
• Please put the parameter definitions in the question formatted as code (use the edit button). Apr 26 at 0:56

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]

• Akku14, Thanks a lot! Apr 27 at 2:35