# Using NDSolveValue to find maximum difference between two functions

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 == 1, f2 == 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][]];
{
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][]
}
]


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][]];
{
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 == 1, f2 == 1,
L1'[t] == Abs[f1[t] - f2[t]], L1 == 0,
L2'[t] == Abs[f1[t] - f2[t]]^2, L2 == 0
},
{L1, L2^(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?

• Maybe you can add WhenEvent[f1'[t]-f2'[t]==0, Sow[Abs[f1[t]-f2[t]]]] and Reap the extrema? Jun 13, 2017 at 23:16
• Please provide the specific {a1, a2, b1, b2} that you wish used for timing purposes. G3 generates error messages for {1, 2, 3, 4}. Jun 13, 2017 at 23:44
• All the time goes into doing NIntegrate, and G3 probably is faster than G2, because the integration is coarser. Jun 14, 2017 at 0:03
• For fixed test values, try 0.1, 0.2, 0.05, 0.07 respectively. For randomized testing, try RandomReal[{0.1, 10}], RandomReal[{0.1, 10}], RandomReal[{0.01, 1}], RandomReal[{0.01, 1}] . Jun 14, 2017 at 0:14
• Sorry, please reverse those: For fixed test values, try 0.05, 0.07, 0.1, 0.2 respectively. For randomized testing, try RandomReal[{0.01, 1}], RandomReal[{0.01, 1}], RandomReal[{0.1, 10}], RandomReal[{0.1, 10}] . I had "a" and "b" mixed-up when I made those up, although they might also work. Again, this is not the actual set of DEs and I'm not looking to optimize beyond solving the problem mentioned. Jun 14, 2017 at 0:31

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[] - v2[]] c[]/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[] - v2[]]^2 c[]/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.

{* 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 *)

• Thank you for your input. I wonder, have you modified G3 in some way so that F[G3[ -- ]] is syntactically correct? (See: 5th line of first code block.) Jun 14, 2017 at 1:19
• @KellenMyers Typo, now fixed. Thanks. Jun 14, 2017 at 1:54

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] == (max = Abs[f1 - f2]) (* diff. of ICs *)}


and return the value

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


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] == (max = Abs[f1 - f2]) (* diff. of ICs *)};
pode[pp_] := {norm[pp]'[t] == Abs[f1[t] - f2[t]]^pp, norm[pp] == 0};
(* code for the return value of the p-norm in NDSolveValue[] *)
pnorm[Infinity] := Max[norm[Infinity], Abs[f1 - f2]];
pnorm[pp_] := norm[pp]^(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 == 1, f2 == 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}}  *)

• Is the line  norm[Infinity] == (max = Abs[1 - 1])  correct? Do you want f1-f2 instead of 1-1? Probably obvious, but I want to make sure I follow. Jun 14, 2017 at 9:12
• @KellenMyers Yes, Abs[f1 - f2] -- I didn't think it would work. I wasn't sure the values would be substituted. Thanks. Jun 14, 2017 at 13:46

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 == 1, f2 == 1,
L1'[t] == Abs[f1[t] - f2[t]], L1 == 0,
L2'[t] == Abs[f1[t] - f2[t]]^2, L2 == 0,
Linf == 0,
WhenEvent[f1'[t] - f2'[t] == 0,
max = Linf[t]; Linf[t] -> Max[max, Abs[f1[t] - f2[t]]]
]
},
{L1, L2^(1/2), Max[Linf, Abs[f1 - f2]]},
{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.

• I did some looking at it seems like it is common to pick one's own answer in a situation like this. Sorry if it is any offense to the other helpful answers, which certainly make up the foundation of my approach to this final answer I came up with. (I will add, in practice, this implementation was more useful for me because I was also working with a much more complex system, with more DEs, no fixed end time, greater complexity/time-costs, etc.) Jun 22, 2017 at 0:10