# How to determine for which initial conditions a system of ODE gives negative results?

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?

• 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. Feb 19 at 15:54
• But you could always also make a Manipulate and try different initial conditions using sliders and see the effect on T solution. Feb 19 at 15:56

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.

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

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}
]

• The error message is when T=0, which causes division by zero -- this is the problem. Feb 19 at 16:30
• @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. Feb 19 at 16:32
• @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. Feb 19 at 16:44
• 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. Feb 19 at 18:09

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}]


• I suspect the crude Method -> {"FixedStep", Method -> "ExplicitRungeKutta"} just hops over the singularity into negativeland. Feb 20 at 4:45
• @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. Feb 20 at 21:10
• 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! Feb 20 at 21:46
• @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. Feb 21 at 14:18