2
$\begingroup$

I have the following coupled differential equations:

ClearAll["Global`*"];

ode1 = p'[t] == -gamma p[t] f[t];
ode2 = f'[t] == -c f[t] + gamma p[t] f[t];
ode3 = T'[t] == b (1 - T[t]) (1 - 0.214499/T[t]) - m (1/(1 + (r/f[t])^n)) (1 - 0.214499/T[t]);

ic = {p[0] == 0.817, f[0] == 0.001, T[0] == 0.850};

param = {gamma -> 0.552, c -> 0.004, b -> 0.093, r -> 0.798, n -> 10.862, m -> 0.870};

ode1 = ode1/.param
ode2 = ode2/.param
ode3 = ode3/.param
ic = ic/.param

NDSolve[{ode1, ode2, ode3, ic}, {p, f, T}, {t, 0, 40}]

For the given initial conditions in ic, the variable $T$ is positive; however, the phase-space plot suggests that there might be other initial conditions that may result in negative $T$. Is there any systematic way to determine for what ics, the variable $T$ becomes negative?

$\endgroup$
2
  • $\begingroup$ any systematic way to determine for what ic's, the variable T becomes negative? This question is not clear to me. Negative where? at any point? at infinity? why not make initial condition for T to be negative? Will this not make T negative at time zero and this meets what you are asking for? Did you try that? Assuming this initial condition still works. $\endgroup$
    – Nasser
    Feb 19 at 15:54
  • $\begingroup$ But you could always also make a Manipulate and try different initial conditions using sliders and see the effect on T solution. $\endgroup$
    – Nasser
    Feb 19 at 15:56

2 Answers 2

4
$\begingroup$

How can I vary the initial conditions and see their effects on T solution?

I mean make a Manipulate, something like this (the plot here is just T[t] and not the other functions.

enter image description here

But notice that for some IC's, you will hit

enter image description here

And this is when you have to adjust time slider or investigate more. But this at least gives you a way to see the effect of IC's on T plot.

Code

ClearAll["Global`*"];

ode1 = p'[t] == -gamma  p[t]  f[t];
ode2 = f'[t] == -c  f[t] + gamma  p[t]  f[t];
ode3 = T'[t] == 
   b  (1 - T[t])  (1 - 0.214499/T[t]) - 
    m  (1/(1 + (r/f[t])^n))  (1 - 0.214499/T[t]);
    
param = {gamma -> 0.552, c -> 0.004, b -> 0.093, r -> 0.798, 
   n -> 10.862, m -> 0.870};

ode1 = ode1 /. param
ode2 = ode2 /. param
ode3 = ode3 /. param

Manipulate[ 
 Module[{ic, sol},
  ic = {p[0] == p0, f[0] == f0, T[0] == T0};
  sol = NDSolve[{ode1, ode2, ode3, ic}, {p, f, T}, {t, 0, maxT}, 
    Method -> {"StiffnessSwitching"}];
  Plot[Evaluate[T[t] /. sol], {t, 0, maxT}, AxesOrigin -> {0, 0}]
  ],
 {{p0, 0.817, "p[0]"}, 0.001, 2, 0.001, Appearance -> "Labeled"},
 {{f0, 0.001, "f[0]"}, 0.001, 2, 0.001, Appearance -> "Labeled"},
 {{T0, 0.850, "T[0]"}, 0.001, 2, 0.001, Appearance -> "Labeled"},
 {{maxT, 40, "time?"}, 0.001, 100, 0.001, Appearance -> "Labeled"},
 TrackedSymbols :> {p0, f0, T0, maxT}
 ]
$\endgroup$
4
  • 1
    $\begingroup$ The error message is when T=0, which causes division by zero -- this is the problem. $\endgroup$
    – Chris K
    Feb 19 at 16:30
  • 1
    $\begingroup$ @Anovice it seams T cannot become negative easily. yes, this is what it looks like. If you try to change IC to make it negative, you get singularity error from NDSolve. May be the Physics of your modeling does not allow T to be negative. If T supposed to mean something physically? Kinetic energy? which can't be negative. $\endgroup$
    – Nasser
    Feb 19 at 16:32
  • 1
    $\begingroup$ @Anovice The model shouldn't allow negative values for T Ok, so NDSolve agrees with you in this case :) I do not think there is an analytical way to show this, since you have to use numerical solver. Only way is to try and see. $\endgroup$
    – Nasser
    Feb 19 at 16:44
  • 2
    $\begingroup$ The only reason $T$ doesn't become negative is because it crashes NDSolve since $dT/dt$ is infinitely negative when $T=0$. This is not good behavior. As can be seen in the phase plane in my answer to your related question, any initial conditions with $p=0$ and $T<0.214499$ will result in $T$ diving towards 0. If $p \neq 0$ then it's more complicated but since $p\to 0$ the basic result holds. This can be seen in @Nasser's answer, where $T<0.214499$ crashes NDSolve. $\endgroup$
    – Chris K
    Feb 19 at 18:09
1
$\begingroup$

There's a singularity when $T=0$, but it seems possible that $T'$ can stay finite.

ode1 = p'[t] == -gamma  p[t]  f[t];
ode2 = f'[t] == -c  f[t] + gamma  p[t]  f[t];
ode3 = T'[t] == 
   b  (1 - T[t])  (1 - a/T[t]) - m  (1/(1 + (r/f[t])^n))  (1 - a/T[t]);

ic = {p[0] == p0, f[0] == f0, T[0] == T0};

param = {gamma -> 0.552, a -> 0.214499, c -> 0.004, b -> 0.093, 
   r -> 0.798, n -> 10.862, m -> 0.870, p0 -> 0.817, f0 -> 0.001, 
   T0 -> 0.21430976608484623`};

sol = NDSolve[{ode1, ode2, ode3, ic} /. param, {p, f, T}, {t, 0, 40}, 
  Method -> {"FixedStep", Method -> "ExplicitRungeKutta"}]


Plot[T[t] /. First@sol, {t, 0, 40}]

enter image description here

$\endgroup$
4
  • 3
    $\begingroup$ I suspect the crude Method -> {"FixedStep", Method -> "ExplicitRungeKutta"} just hops over the singularity into negativeland. $\endgroup$
    – Chris K
    Feb 20 at 4:45
  • 1
    $\begingroup$ @ChrisK I think it's a removable singularity. Both the function and derivative stay bounded. But numerical instability near the singularity drives the standard error estimators crazy, and methods like LSODA reduce the step size until they stop. $\endgroup$
    – Goofy
    Feb 20 at 21:10
  • 1
    $\begingroup$ I'm not sure, for most values of f,T' is discontinuous with an asymptote at T==0. Anyhow, whether T gets stuck at T==0 or crosses into negative values, I think the model isn't well formulated. Real systems shouldn't behave so weird! $\endgroup$
    – Chris K
    Feb 20 at 21:46
  • 1
    $\begingroup$ @ChrisK "I think the model isn't well formulated" — Probably so. The Parker solar wind model (Parker, 1958), which has a similar removable singularity, was criticized on similar grounds. $\endgroup$
    – Goofy
    Feb 21 at 14:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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