Avoid NDSolve problem with system of differential equations?

Questions
,
1

Hello,

I have a system of differential equations which I try to solve for my school project with NDSolve, but I always get an error and i don’t know why.

Here is my code:

{a[t], b[t], c[t]} /. 
 NDSolve[{Derivative[1][a][t] == -0.001 a[t] + 
     2.5*^-8 a[t] (1.2*^6 - a[t] - b[t]) - 0.0014 a[t] c[t], 
   Derivative[1][b][t] == -0.00257 b[t] + 0.0014 a[t] c[t], 
   Derivative[1][c][t] == 
    482 (-12 + t) b[t] - 0.009` c[t] - 0.8 a[t] c[t], 
   a[0] == 1.2*10^6, b[0] == 0, c[0] == 2.6*10^7}, {a, b, c}, {t, 0, 
   90}, MaxSteps -> \[Infinity]]

The program always gives me this error:

NDSolve::nderr: Error test failure at t == 23.653564657637254`; unable to continue.

I am a beginner in the program, what did I do wrong? Thank you so much in advance for all your tips.

2

I suspect that there is an error in your system of equations. The variable c disappears almost immediately.

Clear[a, b, c]
sol = NDSolve[{Derivative[1][a][t] == -0.001 a[t] + 
      2.5*^-8 a[t] (1.2*^6 - a[t] - b[t]) - 0.0014 a[t] c[t], 
    Derivative[1][b][t] == -0.00257 b[t] + 0.0014 a[t] c[t], 
    Derivative[1][c][t] == 
     482 (-12 + t) b[t] - 0.009` c[t] - 0.8 a[t] c[t], 
    a[0] == 1.2*10^6, b[0] == 0, c[0] == 2.6*10^7}, {a, b, c}, {t, 0, 
    0.0001}, MaxStepFraction -> 1/10000];
Plot[Evaluate[{a[t], b[t], c[t]} /. sol], {t, 0, 5*^-6}, 
 PlotRange -> All, PlotLegends -> {"a", "b", "c"}]

Roberta.png
Roberta.png413×201 2.75 KB

3

Hi Tim,

Thank you for your help. When I first entered the equations into the program, it rewrote it into the form I sent.

These are the original equations with initial values:

systemofODEs.jpg
systemofODEs.jpg1008×809 50.8 KB

I need the numerical solution for these equations.

Thank you for your help!

4

Hi Roberta,

The expression may have been rewritten, but you can verify that the are equivalent by using Simplify. Take the right hand side of eqn 1 and compare the two cases for example:

Simplify[-0.001*a[t] + 2.5*^-8*a[t]*(1.2*^6 - a[t] - b[t]) - 
   0.0014*a[t]*c[t] == (0.03*((1.2*^6 - (a[t] + b[t]))/1.2*^6))*a[t] - 
       0.001*a[t] - 0.0014*a[t]*c[t]]
(* True *)

True implies that the two forms are equivalent.

The problem still remains with equation 3. The main driving force for the disappearance of c is the term -0.8 a c, which is of the order -10^{13}, initially. This drives c to zero almost instantaneously. Without any c, there is nothing for equation 2 to do.

5

Thanks for your help, you helped me a lot.

6

Add the option AccuracyGoal -> 14 and the NDSolve integrates through t == 90. However, c[t] spikes negative right away and eventually climbs back up to positive values. That shouldn’t happen with a concentration.

I tested it with the options PrecisionGoal -> 10, AccuracyGoal -> 36, WorkingPrecision -> 50 and the same thing happens. So I think the result is probably an accurate integration of the system.

7

Good morning Roberta.

Are you sure that your system of equations makes sense physically?
What are you describing here?

What is the reason for (t -12) in the expression for dc/dt?

And when I rescale your amounts of material by dividing by 10^6 and let your code run I get an oscillatory behavior of concentration c ( which may well occur, see Belusov-Zhabotinsky -reaction), but it goes to negative values. But negative concentrations don’t make sense.

Clear[a, b, c]
og = 15;
a0 = 1.2;
b0 = 0;
c0 = 2.6 10^1;
sol = NDSolve[
  {Derivative[1][a][
     t] == -0.001 a[t] + .03 a[t] (a0 - a[t] - b[t])/a0 - 
     0.0014 a[t] c[t], 
   Derivative[1][b][t] == -0.00257 b[t] + 0.0014 a[t] c[t], 
   Derivative[1][c][t] == 
    482 (-12 + t) b[t] - 0.009` c[t] - 0.8 a[t] c[t],
   a[0] == a0,
   b[0] == 0,
   c[0] == 2.6*10^1}, {a, b, c}, {t, 0, og}, 
  MaxStepFraction -> 1/10000]

Plot[Evaluate[{a[t], b[t], c[t]} /. sol], {t, 0, og}, 
 PlotRange -> All]
8

Dear Hans,

Thank you for your detailed help!

These equations belong to virus growth in mammalian cells.

Since (t-12) is a time delay, it can be excluded from the code. This can be the reason why it goes to negative interval.

I have recently noticed, that the units of a and b differ from the unit of c, since the unit of a and b is [10^6 cells/ml], until the unit of c is [log HA/100 microliters]. The conversion is 2.4logHA=1.04*10^10 viruses/ml.

9

Dear Michael,

Thank you for your help. I might have the solution, which I have written to Hans, after your answer.