# Finding the maximum of NDSolve's result

I'm solving a set of differential equations as follows:

s = NDSolve[{p'[t] == -gamma p[t] f[t], f'[t] == -c f[t] + gamma p[t] f[t],
T'[t] == b (1 - T[t]) (r^n/(f[t]^n + r^n)) (1 - 0.00262/T[t]) - m (f[t]^n/(f[t]^n + r^n)) (1 - 0.00262/T[t]), p[0] == p0,
f[0] == 0.001, T[0] == T0} /. {gamma -> 0.5431799751931238, b -> 0.08850801132176628, m -> 0.9780307913586779,
r -> 0.8264557187463307, n -> 14.983207992804674, c -> 0.006, p0 -> 0.8, T0 -> 0.8}, {p, f, T}, {t, 0, 72}];

p2 = Plot[Evaluate[{T[t], f[t]} /. s], {t, 0, 72}, PlotStyle -> {Red, Blue}, AxesOrigin -> {-1, 0}, PlotRange -> All, PlotLegends -> {T, f}]


The blue curve corresponds to the variable $$f$$ in the equations. How can I find the coordinates of the maximum of this curve? (And the same for the red curve, its maximum and also its minimum.)

Try this:

FindRoot[(D[T[t], t] /. s) == 0, {t, 20}]

(*  {t -> 17.9042}  *)


or this:

FindMaximum[T[t] /. s, {t, 20}]

(*  {0.956916, {t -> 17.9042}}  *)


and the same with the second function.

Have fun!

• @Anovice: FindInstance[{((f /. s)[[1]][t] - 0.025) == 0, t >= 0 && t <= 72}, t, Reals,2] results in {{t -> 7.57455}, {t -> 7.57455}} and NMinimize[{((f /. s)[[1]][t] - 0.025)^2, t >= 0 && t <= 72}, t] performs {1.14078*10^-19, {t -> 7.57455}}. Commented Mar 21 at 16:42
• @A novice it works for the function f: FindRoot[((f[t] /. s) - 0.0025) == 0, {t, 1}] yielding {t -> 2.13987}. In the case of T(t), there is something weird with the solution. You can see it by warnings. There is a workaround: FindMinimum[((T[t] /. s) - 0.0025)^2, {t, 30}] returning: {1.4399*10^-8, {t -> 31.3609}}. Commented Mar 21 at 18:09
• From where do you know 20 in FindRoot[(D[T[t], t] /. s) == 0, {t, 20}]?  Commented Mar 21 at 18:53
• @user64494 From the plot. Commented Mar 21 at 19:17
• @AlexeiBoubitch: Thank you. NMaximize does need any plot. Commented Mar 21 at 19:29

This can be done as follows.

NMaximize[{(f /. s)[[1]][t], t >= 0 && t <= 72}, t]


{0.742649, {t -> 25.7451}}

NMaximize[{(T /. s)[[1]][t], t >= 0 && t <= 72}, t]


{0.956916, {t -> 17.9042}}

NMinimize[{(T /. s)[[1]][t], t >= 0 && t <= 72}, t]


{0.00261999, {t -> 33.1547}}

Here are three approaches with increasing accuracy that take advantage of the data computed by NDSolve:

(1) The maximum value of T[t] and the value of t at the maximum, among the steps computed by NDSolve:

maxPos = PositionLargest[T["ValuesOnGrid"] /. First[s]];
maxTVal = Extract[T["ValuesOnGrid"] /. First[s], maxPos]
maxt1Val =
First@Extract[T["Grid"] /. First[s], maxPos]
(*
0.9569154933291244
17.89633907171858
*)


(2) The maximum of the interpolating function returned by NDSolve, which occurs between the steps computed by NDSolve but contains but contains additional error from the approximation by the interpolation:

maxTPVal = Extract[T'["ValuesOnGrid"] /. First[s], maxPos]
maxt2Val =
First@Extract[T["Grid"] /. First[s], maxPos + Sign[maxTPVal]]
FindMaximum[T[t] /. First[s], {t, maxt1Val, maxt2Val}]
(*  {0.9569156862237768, {t -> 17.904199030539463}}  *)


One can reduce the error due to interpolation by add the option InterpolationOrder -> All to NDSolve. It affects the 6th digit of precision slightly, so the extra precision may not be needed in some cases.

Replace T by f to get the corresponding results for f[t].

(3) The following computes the maximum using NDSolve (you have to plan ahead and include the code that saves the maximum values). WhenEvent will locate the maximum as accurately as the ODE integration method allows.

maxT = {0., 0.};
{ode, ics} = {
{p'[t] == -gamma  p[t]  f[t]
, f'[t] == -c  f[t] + gamma  p[t]  f[t]
, T'[t] ==
b  (1 - T[t])  (r^n/(f[t]^n + r^n))  (1 - 0.00262/T[t]) -
m  (f[t]^n/(f[t]^n + r^n))  (1 - 0.00262/T[t])}
, {p[0] == p0, f[0] == 0.001, T[0] == T0}
} /. {gamma -> 0.5431799751931238, b -> 0.08850801132176628,
m -> 0.9780307913586779, r -> 0.8264557187463307,
n -> 14.983207992804674, c -> 0.006, p0 -> 0.8, T0 -> 0.8};
derivSol = First@Solve[ode, {T'[t], f'[t], p'[t]}];
s = NDSolve[{ode, ics,
WhenEvent[T'[t] < 0, If[T[t] > maxT[[2]], maxT = {t, T[t]}];
"RestartIntegration"] /. derivSol,
WhenEvent[f'[t] < 0, maxf = {t, f[t]}; "RestartIntegration"] /.
derivSol}
, {p, f, T}, {t, 0, 72}];

maxT
maxf
(*
{17.904137483723822, 0.9569156846482132}
{25.74518018578614, 0.7426487539346786}
*)


You could use a similar If[] statement to save the maximum of f[t], if f[t] might not have just one maximum. (T[t] has several maximum due to jitter when we have T[t] == 0 approximately, between around t == 30 and t == 40.)

Remark: In fact, the length of the interval where T[t] is approximately zero can vary depending on PrecisionGoal, AccuracyGoal, WorkingPrecision and WhenEvent[]. It makes me wonder if the jump up from T[t] == 0 in the OP's plot is an anomaly due to rounding error.