MathGroup Archive 2009

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Delay Differential Equations

  • To: mathgroup at smc.vnet.net
  • Subject: [mg105272] Re: [mg105232] Re: [mg105198] Delay Differential Equations
  • From: DrMajorBob <btreat1 at austin.rr.com>
  • Date: Wed, 25 Nov 2009 23:01:28 -0500 (EST)
  • References: <200911241048.FAA00661@smc.vnet.net>
  • Reply-to: drmajorbob at yahoo.com

Using Daniel Lichtblau's methods, we can find out when the type A  
population reaches zero:

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

Clear@p
p = p /. First@
     NDSolve[{p'[t] == -0.0367*p[t - 37]*UnitStep[t - 37] -
         0.0123*p[t], p[t /; t <= 0] == 30000}, p, {t, 0, 120},
      Method -> {"EventLocator", "Event" -> p[t]}];
end = InterpolatingFunctionDomain[p][[1, -1]]
Plot[p@t, {t, 0, end}, AxesOrigin -> {0, 0}]

54.2857

But this allows dead cells to mutate (and reduce the number of type A  
cells), so it's probably not what you really want.

Regardless of cell type, this gives the number ALIVE at time t:

Clear[alive, p]
p = p /. First@NDSolve[{
       p'[t] == -0.0123*p@t, p[0] == 30000}, p, {t, 0, 120}];
alive[t_?NonNegative] := p[t];
alive[t_?Negative] = 0;
Plot[alive@t, {t, -37, 120}, AxesOrigin -> {0, 0}]

and this purports to determine the proportion of living (or dead) cells  
that are type A, at time T:

Clear[typeA]
typeA = typeA /. First@NDSolve[
      {typeA'[t] == -0.0367*typeA[t - 37],
       typeA[t /; t <= 0] == UnitStep@t}, typeA, {t, 0, 120}];
FindRoot[typeA@t, {t, 0, 60}]
Plot[typeA@t, {t, 0, 120}, AxesOrigin -> {0, 0}]

{t -> 64.248}

But that function goes negative as well, and not at the same time as  
before!

So... Daniel Huber's simpler take is probably the right one:

Clear[p]
p = p /. NDSolve[{p'[t] == -0.0123 p[t] - 0.0367 Boole[t > 37] p[t],
       p[0] == 30000}, p, {t, 0, 120}][[1]];
Plot[p@t, {t, 0, 120}]

We can get a symbolic solution in this case, too:

Clear[p]
p = p /. First@
     DSolve[{p'[t] == -0.0123 p[t] - 0.0367 Boole[t > 37] p[t], p[0] ==  
30000},
       p, t];
p[t] // InputForm

Piecewise[{{30000/E^(0.0123*t), t <= 37}},
  116640.59642871171/E^(0.049*t)]

(Accurate for t>=0.)

Bobby

On Wed, 25 Nov 2009 01:31:53 -0600, Daniel Lichtblau <danl at wolfram.com>  
wrote:

> Yun Zhao wrote:
>> Hi Daniel,
>>
>>   Thanks for your reply.  But my question is I was trying to solve an
>> equation, and the solution just doesnt make much sense, and I looked in
>> the help files, but they didn't really help me much.  The problem is  
>> this:
>>
>> I have some cells of type A, and at some point in time, cells of type A
>> are going to slowly change into cells of type B.  So, population size of
>> cell type A will decrease over time due to two factors: (1) becoming
>> cells of type B and (2) dying over time.  So for a differential
>> equation, i have two rates, one is the transformation rate (0.0367), one
>> is the death rate (0.0123).  But I want to delay the transformation by
>> 37 hours.  So from time = 0 to time = 37, all that is changing in cells
>> of type A is that they are dying at the death rate, then after 37 hours,
>> they are also disappearing at the transformation rate.
>>
>> So my differential equation is like this:
>>
>> *solrandom1=NDSolve[{p'[t]*
>>
>> *==0.0367*p[t-37]-0.0123*p[t],p[t/;t<0]=30000},p,{t,0,120}]*
>>
>> What I expect when I plot p(t) over t is a gradual decrease from t=0 to
>> t=37.  Then more decrease from t=37 to t=120.  p(t) should never be less
>> than zero.  But when I ran that code in Mathematica, and plotted it
>> using this command
>>
>> *Plot[Evaluate[p[t]/.solrandom1],{t,0,120},PlotRange**ïâ??®Automatic]*
>>
>> p(t) was actually less than zero.  I am not able to figure out why this
>> is.
>>
>> Please tell me what I did wrong.  Thanks.
>
>
> First is that the delay term you have is not negative. So I'd expect an
> increase. When I do the following I indeed see a rising population.
>
> solrandom1 =
>   NDSolve[{p'[t] == 0.0367*p[t - 37] - 0.0123*p[t],
>     p[t /; t <= 0] == 30000}, p, {t, 0, 120}]
>
> Plot[Evaluate[p[t] /. solrandom1], {t, 0, 120},
>   PlotRange -> Automatic]
>
> Second is that your claim of no transformations for 37 hours is wrong.
> You are modeling that by 0.0367*p[t - 37], and p[t] is 30000 for t<=0.
> So this contribution is not going to be zero (whether you use a positive
> or negative thereof). You can repair this by using the system below.
>
> solrandom1 =
>   NDSolve[{p'[t] == -0.0367*p[t - 37]*UnitStep[t - 37] - 0.0123*p[t],
>     p[t /; t <= 0] == 30000}, p, {t, 0, 120}]
>
> Next is that there is nothing to stop this solution from going negative.
> For that you'd need an EventHandler. Could do as below.
>
> solrandom1 =
>   NDSolve[{p'[t] == -0.0367*p[t - 37]*UnitStep[t - 37] - 0.0123*p[t],
>     p[t /; t <= 0] == 30000}, p, {t, 0, 120},
>    Method -> {"EventLocator", "Event" -> p[t]}]
>
> Last is that these methods are things in the documentation I had noted
> in my last reply.
>
> Daniel Lichtblau
> Wolfram Research
>
>> On Tue, Nov 24, 2009 at 1:28 PM, Daniel Lichtblau <danl at wolfram.com
>> <mailto:danl at wolfram.com>> wrote:
>>
>>     Yun Zhao wrote:
>>
>>         Hi,
>>
>>         Does anyone have experience working with delay differential
>>         equations in
>>         Mathematica 7? I found some help files on wolfram and
>>         Mathematica help, but
>>         the information available there were very limited. If anyone can
>>         refer me to
>>         other online or textual sources, I would really appreciate it.
>>         Thank you
>>         very much.
>>
>>         Mike
>>
>>
>>     Not sure what you saw in the Documentation Center. These are perhaps
>>     the most useful, worth checking if you have not encountered them
>>     already.
>>
>>     howto/SolveDelayDifferentialEquations
>>     tutorial/NDSolveDelayDifferentialEquations
>>
>>     Daniel Lichtblau
>>     Wolfram Research
>>
>>
>
>


-- 
DrMajorBob at yahoo.com


  • Prev by Date: Re: part assigned sequence behavior puzzling
  • Next by Date: Re: Kernel crashes in ReplaceAll - with or without Maximize
  • Previous by thread: Re: Delay Differential Equations
  • Next by thread: Re: Delay Differential Equations