Using NDSolveValue to find maximum difference between two functions

Question

I have a set of differential equations I'm trying to solve with various parameters and compare some of the functions. Here is the naive approach:

F[a1_, a2_, b1_, b2_] :=
 NDSolve[{
  f1'[t] == b1 f1[t] - a1 f1[t] f2[t],
  f2'[t] == b2 f2[t] - a2 f1[t] f2[t],
  f1[0] == 1, f2[0] == 1},
 {f1, f2},
 {t, 0, 10}
]

(Note: These are not my DEs. Just a trivial example. Additional complications of my DEs should not affect how we approach the question I have.)

So my goal is to find the difference between the functions. Proceeding naively, we can just integrate to find any $p$-norms we want, except $p=\infty$, which I've used Maximize for instead:

G[a1_, a2_, b1_, b2_] :=
 Block[{s1, s2},
  {s1, s2} = ReplaceAll[{f1, f2}, F[a1, a2, b1, b2][[1]]];
  {
   NIntegrate[Abs[s1[t] - s2[t]], {t, 1, 10}],
   NIntegrate[Abs[s1[t] - s2[t]]^2, {t, 1, 10}]^(1/2),
   Maximize[{Abs[s1[t] - s2[t]], 0 <= t <= 10}, t][[1]]
   }
  ]

However, this is quite slow. The faster version I have can only compute the norms for $p<\infty$ by making a straightforward adjustment to combine the functions F and G into one. So first, for comparison, strip out the extra stuff from G to get less output:

G2[a1_, a2_, b1_, b2_] :=
 Block[{s1, s2},
  {s1, s2} = ReplaceAll[{f1, f2}, F[a1, a2, b1, b2][[1]]];
  {
   NIntegrate[Abs[s1[t] - s2[t]], {t, 1, 10}],
   NIntegrate[Abs[s1[t] - s2[t]]^2, {t, 1, 10}]^(1/2)
   }
  ]

And here's our new G3:

G3[a1_, a2_, b1_, b2_] :=
 NDSolveValue[
  {
   f1'[t] == b1 f1[t] - a1 f1[t] f2[t], 
   f2'[t] == b2 f2[t] - a2 f1[t] f2[t],
   f1[0] == 1, f2[0] == 1,
   L1'[t] == Abs[f1[t] - f2[t]], L1[0] == 0,
   L2'[t] == Abs[f1[t] - f2[t]]^2, L2[0] == 0
   },
  {L1[10], L2[10]^(1/2)},
  {t, 0, 10}
  ]

The new G3 is about 100x faster than G2, but it's missing the $p=\infty$ norm that I'd like to have. I've run out of ideas, and I have an itch in the back of my brain that says "it's a simple trick you're overlooking." Something internal to NDSolve and NDSolveValue that lets me track the value of some quantity at each time-step and update it (i.e. Abs[s1[t]-s2[t]])? This sounds like something that may exist, but I'm not finding it in the web of documentation.

The only practical idea I came up with was just doing something like $p=1000$ as well as $p=1,2$, which could be an approximate $p=\infty$ norm. It fits nicely in the existing code for G3 and returns fairly useful data, but it's actually slower than the original G.

My intention is to use this to compare over, say, hundreds or thousands values of each of a1, b1, a2, b2, which means I do not have the luxury of saying "1 second is fine, even if I could do it in 0.05 seconds."

Any thoughts on getting the $p=\infty$ norm at similar speeds as we find the other two quantities?

Answer 1

From RepeatTiming of the various functions for the test parameters given in a comment,

G2[.05, .07, .1, .2] // RepeatedTiming
(* {0.0386, {8.42801, 3.41397}} *)
G3[.05, .07, .1, .2] // RepeatedTiming
(* {0.0010, {8.47102, 3.41433}} *)
s = Flatten@F[.05, .07, .1, .2]; // RepeatedTiming
(* {0.00083, Null} *)

we see that G3 is about 40 times faster G2 but its results differ from those of G2 by about several tenths of a percent. More importantly, F is about the same speed as G3, indicating that essentially all of the extra time required by G2 goes to NIntegrate. This suggests that G3 is about the best that can be done, apart from reducing the accuracy required of NDSolve.

Here is an alternative approach.

c = (f1 /. s)["Coordinates"] // Flatten;
v1 = (f1 /. s)["ValuesOnGrid"];
v2 = (f2 /. s)["ValuesOnGrid"];

(* p == 1 *)
(Sum[Abs[v1[[i]] - v2[[i]]] (c[[i + 1]] - c[[i - 1]])/2, {i, 2, Length[c] - 1}]
    + Abs[v1[[1]] - v2[[1]]] c[[1]]/2 + Abs[v1[[-1]] - v2[[-1]]] (c[[-1]] - c[[-2]])/2) 
    // RepeatedTiming
(* {0.00011, 8.48083} *)

(* p == 2 *)
(Sum[Abs[v1[[i]] - v2[[i]]]^2 (c[[i + 1]] - c[[i - 1]])/2, {i, 2, Length[c] - 1}] 
    + Abs[v1[[1]] - v2[[1]]]^2 c[[1]]/2 + Abs[v1[[-1]] - v2[[-1]]]^2 (c[[-1]] - c[[-2]])/2)
    ^(1/2) // RepeatedTiming
(* {0.00012, 3.42291} *)

This alternative, which simply uses the trapezoidal integration rule on the raw data produced by NDSolve is at least as fast and as accurate as G3. Moreover, it offers the advantage of separating the computation of {f1, f2} from the computation of the norms. Note that G3 and the approach just given must have the same intrinsic accuracy, because they both rely directly on the same raw data from NDSolve.

Addition

{* p == 10 *)
(* {0.00012, 2.20872} *)
(* p == 100 *)
(* {0.00012, 2.33251} *)
(* p == 1000 *)
(* {0.00017, 2.366285331305866415} *)
(* p == 10000 *)
(* {0.00024, 2.3696899749044680823} *)

The last few compare well with the p == Infinity norm for the array Abs[v1 - v2], namely

(* p == Infinity *)
Max[Abs[v1 - v2]]
(* 2.37007 *)

Answer 2

A way to get the infinity norm straight out of NDSolveValue is with

{WhenEvent[f1'[t] - f2'[t] == 0,                      (* critical point *)
  max = norm[Infinity][t];                            (* update  max and *)
  norm[Infinity][t] -> Max[max, Abs[f1[t] - f2[t]]]   (* norm[Infinity] *)
  ],
 norm[Infinity][0] == (max = Abs[f1[0] - f2[0]]) (* diff. of ICs *)}

and return the value

Max[norm[Infinity][10], Abs[f1[10] - f2[10]]]

from NDSolveValue. (Here max is a non-diff.-eq. variable that should be initialized to the absolute difference of the functions at the initial condition.)

Here's a function that computes a list of norms p:

ClearAll[G4];
G4[a1_, a2_, b1_, b2_, p_: {1, 2, Infinity}] :=
  Module[{normlist, norm, max, f1, f2, t, pode, pnorm, discvar},
   (* code to compute p-norm in NDSolve[] *)
   pode[Infinity] := {WhenEvent[f1'[t] - f2'[t] == 0,
      max = norm[Infinity][t]; 
      norm[Infinity][t] -> Max[max, Abs[f1[t] - f2[t]]]],
     norm[Infinity][0] == (max = Abs[f1[0] - f2[0]]) (* diff. of ICs *)};
   pode[pp_] := {norm[pp]'[t] == Abs[f1[t] - f2[t]]^pp, norm[pp][0] == 0};
   (* code for the return value of the p-norm in NDSolveValue[] *)
   pnorm[Infinity] := Max[norm[Infinity][10], Abs[f1[10] - f2[10]]];
   pnorm[pp_] := norm[pp][10]^(1/pp);
   (* DiscreteVariables *)
   discvar[{___, Infinity, ___}] := {norm[Infinity]};   
   discvar[_] := {Nothing};

   normlist = Flatten[{p}];
   NDSolveValue[{
     f1'[t] == b1 f1[t] - a1 f1[t] f2[t], 
     f2'[t] == b2 f2[t] - a2 f1[t] f2[t], f1[0] == 1, f2[0] == 1,
     pode /@ normlist},
    pnorm /@ normlist,
    {t, 0, 10},
    DiscreteVariables -> discvar[normlist]]
   ];

Include the infinity norm makes it about three times slower, but only three times. It's going to be faster, I think, to Reap and Sow[Abs[f1[t] - f2[t]]] at critical points and then post-process the reaped local extrema, as Carl Woll suggested.

G4[0.3, -0.2, 0.5, -0.1] // RepeatedTiming      (* default is p = {1, 2, Infinity} *)
(*  {0.0044, {11.4175, 4.74135, 2.45422}}  *)

G4[0.3, -0.2, 0.5, -0.1, {1, 2}] // RepeatedTiming
(*  {0.0016, {11.4175, 4.74135}}  *)

Answer 3

Here is what I came up with, after trying a few different implementations and combinations of the ideas proposed by Michael E2, bbgodfrey, and Carl Woll.

G5[a1_, a2_, b1_, b2_] :=
 Block[{max = 0},
  NDSolveValue[
   {
    f1'[t] == b1 f1[t] - a1 f1[t] f2[t], 
    f2'[t] == b2 f2[t] - a2 f1[t] f2[t], f1[0] == 1, f2[0] == 1,
    L1'[t] == Abs[f1[t] - f2[t]], L1[0] == 0,
    L2'[t] == Abs[f1[t] - f2[t]]^2, L2[0] == 0,
    Linf[0] == 0,
    WhenEvent[f1'[t] - f2'[t] == 0,
     max = Linf[t]; Linf[t] -> Max[max, Abs[f1[t] - f2[t]]]
     ]
    },
   {L1[10], L2[10]^(1/2), Max[Linf[10], Abs[f1[10] - f2[10]]]},
   {t, 0, 10},
   DiscreteVariables -> {Linf}
   ]
  ]

This seems very competitive with G3, taking about 50-75% more time (quite reasonable considering we are asking for more information).

I am not sure how to choose best answer, as the two answers provided valuable input which I distilled into this -- which I think represents a sort of minimal change from G3 to get what I was looking for.

Without choosing a best answer yet, I should thank the three contributors. Your input was very helpful and valuable. Perhaps I will leave this open a bit longer in case anyone else has another interesting perspective to add.