1
$\begingroup$

I want to solve a system of delay equations (DDE), that in some periods of time receive an impulse to one of its equations.

1

2

3

4

5

The equation (2) is the one receiving impulses of size $10^6$. So, I'm trying to use NDsolve of solve numerically this system with initial conditions $T=6*10^4$ ,$C_i=2*10^6$,$C_a=0$,$D=0$,$F_\beta = 0$ for $t<=0$.

dTdelay = 
  rt T[t] Log[K/T[t]] - 
   at Ca[t] T[t] (((1 - atb) etb)/(etb + Fb[t]) + atb);
dCadelay = 
  ra Exp[-μci τ] Ci[t - τ] (Den[t - τ]/(
     Den[t - τ] + θd)) + 
   re Exp[-μca τ] Den[t - τ] (Ca[t - τ]/(
     Ca[t - τ] + θa)) - μca Ca[t];
dCidelay = -ra Exp[-μci τ] Ci[t - τ] (Den[t - τ]/(
     Den[t - τ] + θd)) - μci Ci[t];
dFbdelay = rtb T[t] - Fb[t] μb;
dDendelay = (-Den[t] μd);

param3 = {rt -> 0.001060289, K -> 1*10^11, at -> 6*^-11, ht -> 5*10^8,
    atb -> 0.69, etb -> 10^4, 
   re -> 6500, θd -> 10, θa -> 212, 
   ra -> 61, μca -> 0.01925, μci -> 0.007, 
   rtb -> 5.57*10^-6, μb -> 2.7, μd -> 0.009625};

dose = 10^6; τ = 265; tf = 1000;

teraphy2 = {WhenEvent[t == 168, Den[t] -> Den[t] + dose], 
   WhenEvent[t == 336, Den[t] -> Den[t] + dose], 
   WhenEvent[t == 504, Den[t] -> Den[t] + dose]};

sys = {T'[t] == dTdelay , Ca'[t] == dCadelay , Ci'[t] == dCidelay , 
   Fb'[t] == dFbdelay , Den'[t] == dDendelay , 
   T[t /; t <= 0] ==  6*^4, Ca[t /; t <= 0] ==  0, 
   Ci[t /; t <= 0] ==  2*^6, Fb[t /; t <= 0] ==  0, 
   Den[t /; t <= 0] ==  0};

Also i'm using WhenEvent to add the impulses at $t=168$ , $t=336$ , $t=504$.But NDsolve returns:

sol = NDSolve[
   Join[sys, teraphy2] /. param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf}];

NDSolve::ndsz: At t == 432.9999999999877`, step size is effectively zero; singularity or stiff system suspected. >>

I know that the problems happen in $t=433$ when the delayed equations (3) and (4) start receiving the impulsed from $D[t-\tau]$ at $t=168$. Then suddenly (3) and (4) change. And that seems to cause problems to NDsolve.

NDsolve works fine until $t=433$ :

    sol = NDSolve[
   Join[sys, 
     teraphy2, {WhenEvent[t == τ + 168, "StopIntegration"]} ] /. 
    param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf}];
    Plot[Evaluate[{T[t], Ca[t], Ci[t], Fb[t], Den[t]} /. First[sol]], {t, 
      0, 433}, PlotRange -> All, 
     PlotLegends -> Placed[{"T", "Ca", "Ci", "Fβ", "D"}, Above],
     Epilog -> {PointSize[0.02], Red, Apply[Point, {{168, 0}} ]} ] 

Which returns:

I'm trying to overcome the problem integrating until $t=433$, save the data as interpolated functions. Using NDsolve to restart integrating at $t=433$ but im still receiving.

NDSolve::ndsz: At t == 432.9999999999877`, step size is effectively zero; singularity or stiff system suspected. >>

Any help of suggestion is welcome.

$\endgroup$
2
$\begingroup$

NDSolve will solve it if you decrease the AccuracyGoal but I am worried that your model may have some underlying problems. First, a solution up to t=500:

tf = 500;
sol = NDSolve[Join[sys, teraphy2] /. param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf},
  AccuracyGoal -> 5];
Plot[Evaluate[{T[t], Ca[t], Ci[t], Fb[t], Den[t]} /. First[sol]], {t, 0, tf},
  PlotRange -> All, PlotLegends -> Placed[{"T", "Ca", "Ci", "Fβ", "D"},
  Above]]

Mathematica graphics

That negative Ci worries me. Could be a numerical problem, but I think it might not be. Take a look at the first term in eqn 4, dCidelay. Ci decreases due to Ci τ time ago. Since Ci is decreasing with time, there is nothing preventing Ci from becoming negative, as it does. The following graphs support this argument:

Plot[Evaluate[{Ci[t], Ci[t - τ]} /. First[sol]], {t, 400, tf}, PlotRange -> {0, All}]
Plot[ra Exp[-μci τ] Ci[t - τ] (Den[t - τ]/(Den[t - τ] + θd)) /.
  param3 /. sol[[1]], {t, 400, tf}]

Mathematica graphics Mathematica graphics

If you let the next pulse hit, all hell breaks loose:

tf = 710;
sol = NDSolve[Join[sys, teraphy2] /. param3, {T, Ca, Ci, Fb, Den}, {t, 0, tf}, AccuracyGoal -> 5];
Plot[Evaluate[{T[t], Ca[t], Ci[t], Fb[t], Den[t]} /. First[sol]], {t, 400, tf},
  PlotRange -> All, PlotLegends -> Placed[{"T", "Ca", "Ci", "Fβ", "D"}, Above]]

Mathematica graphics

So, my guess is it's more of a model problem than an NDSolve problem.

$\endgroup$

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.