0
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ Please provide the values of all parameters, so that readers can run your code. $\endgroup$
    – bbgodfrey
    Apr 25 at 22:58
  • $\begingroup$ 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; $\endgroup$ Apr 25 at 23:44
  • $\begingroup$ Your problem is that (1 - Exp[-k t]) is zero at t = 0. By the way, nn remains undefined. $\endgroup$
    – bbgodfrey
    Apr 26 at 0:36
  • $\begingroup$ k also is undefined. $\endgroup$
    – bbgodfrey
    Apr 26 at 0:47
  • 1
    $\begingroup$ Please put the parameter definitions in the question formatted as code (use the edit button). $\endgroup$
    – Michael E2
    Apr 26 at 0:56

1 Answer 1

1
$\begingroup$

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]
$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.