4
$\begingroup$

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

enter image description here

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.)

$\endgroup$

3 Answers 3

7
$\begingroup$

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!

$\endgroup$
5
  • 1
    $\begingroup$ @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}}. $\endgroup$
    – user64494
    Commented Mar 21 at 16:42
  • 1
    $\begingroup$ @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}}. $\endgroup$ Commented Mar 21 at 18:09
  • 1
    $\begingroup$ From where do you know 20 in FindRoot[(D[T[t], t] /. s) == 0, {t, 20}]? $\endgroup$
    – user64494
    Commented Mar 21 at 18:53
  • 1
    $\begingroup$ @user64494 From the plot. $\endgroup$ Commented Mar 21 at 19:17
  • 1
    $\begingroup$ @AlexeiBoubitch: Thank you. NMaximize does need any plot. $\endgroup$
    – user64494
    Commented Mar 21 at 19:29
3
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

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