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?
WhenEvent[f1'[t]-f2'[t]==0, Sow[Abs[f1[t]-f2[t]]]]
andReap
the extrema? $\endgroup${a1, a2, b1, b2}
that you wish used for timing purposes.G3
generates error messages for{1, 2, 3, 4}
. $\endgroup$NIntegrate
, andG3
probably is faster thanG2
, because the integration is coarser. $\endgroup$0.1, 0.2, 0.05, 0.07
respectively. For randomized testing, tryRandomReal[{0.1, 10}], RandomReal[{0.1, 10}], RandomReal[{0.01, 1}], RandomReal[{0.01, 1}]
. $\endgroup$0.05, 0.07, 0.1, 0.2
respectively. For randomized testing, tryRandomReal[{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. $\endgroup$