# Speed up a differential equation solution + parameter integration + parameter maximization

The problem I am trying to solve can be described as a maximization of an integral of the solution of a differential equation. A MWE is

ClearAll["Global*"]
\$Assumptions =  a > 0 && n > 1;
sol[a_, n_] := NDSolve[{X'[u] == a^2/2 ((1 + n) (1 - 2 X[u]) + Sqrt[4 n + (1 - n)^2 (1 - 2 X[u])^2]), X[0] == 1/2}, X[u], {u, 0, 1}];
P[a_, n_?NumericQ, t_] := (X[u] /. sol[a, n] /. {u -> t})[[1]];
Pmax[a0_, n0_] := Module[{a = a0, n = n0},
f[m_?NumericQ] := NIntegrate[P[a, n - m, 1] E^-(m  a^2 + b^2) 2 b BesselI[0, 2 Sqrt[m] a b], {b, 0, ∞}];
maxIntegral = NMaximize[f[m], {m} ∈ Interval[{0., n}]]; {maxIntegral[[1]], maxIntegral[[2, 1, 2]]}];
Pmax[1.5, 2] // Timing (* gives me {420.376, {0.998829, 4.47458*10^-12}} *)


How can I speed up the evaluation? I am interested in Pmax ranging 0<a<2 and 1<n<20, I suspect there will be numerical issues as n increases.

In addition, is there a better practice to formulate the problem in Mathematica?

EDIT: As suggested in an answer, in this particular MWE the term P[a, n - m, 1] can be left out of NIntegrate, leading to a simplification of the evaluation (the integral gives 1). In the original problem, P also depends on b, so I am looking to more computationally-efficient ways to perform the nested sequence of NDSolve, NIntegrate and NMaximize.

• Look in the options of NIntegrate and NMaximize. Dec 4, 2020 at 12:45

The easiest speed-up is adding the option Method -> "NelderMead" to NMaximize. Then I get the following, which happens to be the exact (most precise possible) solution:

Pmax[1.5, 2] // AbsoluteTiming
(*  {4.38104, {0.998829, 0.}}  *)


If the objective function is unimodular, use FindMaximum instead of NMaxmize for what is probably the fastest way:

sol[a_, n_] :=
NDSolve[{X'[u] ==
a^2/2 ((1 + n) (1 - 2 X[u]) +
Sqrt[4 n + (1 - n)^2 (1 - 2 X[u])^2]), X[0] == 1/2},
X[u], {u, 0, 1}];
P[a_, n_?NumericQ, t_] := (X[u] /. sol[a, n] /. {u -> t})[[1]];
Pmax[a0_, n0_] :=
Module[{a = a0, n = n0},
f[m_?NumericQ] :=
NIntegrate[
P[a, n - m, 1] E^-(m a^2 + b^2) 2 b BesselI[0,
2 Sqrt[m] a b], {b, 0, \[Infinity]}];
maxIntegral = FindMaximum[f[m], {m, 1/2, 0, 1}];
{First@maxIntegral, m /. Last@maxIntegral}
];

Pmax[1.5, 2] // AbsoluteTiming
(*  FindMinimum::reged... warning message *)
(*  {0.120378, {0.998829, 0.}}  *)


It gives a FindMinimum::reged warning ("at the edge of the search region"), which is to be expected for a monotonic objective function.

In this toy example, f[m] is strictly decreasing and the maximum is at m -> 0. For some reason, the optimizers struggle to reach zero and stay within Interval[{0, n}]. There must be something about the penalty function used to enforce the constraints that causes trouble. Even FindMaximum[{f[m], 0 <= m <= 1}, {m, 1/2}] has some trouble. It takes ~0.6 sec. and obtains m -> 0.0000708571, which is at best okay. If you add PrecisionGoal -> 10, it takes longer than I would wait.

Here my solution attempt

In the definition f[m_]:=… one can leave out P[a, n - m, 1] because it doesn't depend on b. The remaining Integral

Integrate[E^-(m a^2 + b^2) 2 b BesselI[0, 2 Sqrt[m] a b], {b, 0, \[Infinity]},Assumptions -> Element[m, Reals]]
(*ConditionalExpression[1, m != 0]*)


equals 1 for real m(m!=0)!

Modified definition Pmax

Pmax[a0_, n0_] :=Module[{a = a0, n = n0}, f[m_?NumericQ] := P[a, n - m, 1] ;
maxIntegral =NMaximize[
f[m], {m} \[Element] Interval[{0., n}]]; {maxIntegral[[1]],maxIntegral[[2, 1, 2]]}]


now is evaluated

Pmax[1.5, 2] // Timing
(*{32.8125, {0.998829, -6.705*10^-26}}*)


in 32seconds! (<< 420s)

• Indeed in this MWE you can leave out P[a,n-m,1] and the integral (over a probability distribution) gives 1. This simplifies a lot the evaluation, but in the original problem P also depend on b (I didn't keep it here in order to formulate a "simple" MWE). I will edit the question to clarify that. I wonder: are there other ways to speed up the evaluation? I am studying compile and other methods but I am not familiar with them. Dec 3, 2020 at 13:07
• Unfortunately NMaximize and NIntegrate cannot be compiled. Dec 3, 2020 at 13:31

Use a detour with Plot , Reap  the plotted points and generate an interpolating function for f[m] to get maximum in half a second.

See Edit at end

Xsol[a_?NumericQ, n_?NumericQ] :=
X /. First@
NDSolve[{Derivative[1][X][u] ==
a^2/2 ((1 + n) (1 - 2 X[u]) +
Sqrt[4 n + (1 - n)^2 (1 - 2 X[u])^2]), X[0] == 1/2},
X, {u, 0, 1}]

Xsol[3, 2][.3]

nint[m_, a_, n_] :=
NIntegrate[
Xsol[a, n - m][1] E^-(m a^2 + b^2) 2 b BesselI[0,
2 Sqrt[(m)] a b], {b, 0, \[Infinity]}]


Get maximum with Plot and see what is going on.

Manipulate[{pl2 =
Plot[nint[m, a, n], {m, 0, n}, PlotRange -> {.45, 1.05},
ImageSize -> 400], {"Max=",
Max[pl2[[1, 1, 3, 2, 1]][[All, 2]]]}}, {{a, 1}, 0, 2,
Appearance -> "Labeled"}, {{n, 3}, 1, 20, Appearance -> "Labeled"},
ContinuousAction -> False]

rp[a_, n_] :=
Reap[Plot[fr = nint[m, a, n], {m, 0, n}, PlotPoints -> 30,
MaxRecursion -> 1, EvaluationMonitor :> Sow[{m, fr}]]][[2, 1]] //
Sort

f[a_?NumericQ, n_?NumericQ] := Interpolation[rp[a, n]]

f[1, 3][1]

Plot[Evaluate[f[1, 3][m]], {m, 0, 3}, PlotRange -> All]

max[a_, n_] := Maximize[{f[a, n][m], 0 <= m <= n}, m]

nmax[a_, n_] := NMaximize[{f[a, n][m], 0 <= m <= n}, m]

{max[2, 2] // AbsoluteTiming, nmax[2, 2] // AbsoluteTiming}

(*   {{0.4843818, {0.999989, {m -> 0.}}}, {0.5312506, {0.999989, {m ->
0.}}}}   *)


Maximize and NMaximize need the same time. You gain a lot of speed due to MaxRecursion -> 1  in Plot, but accuracy seems not to suffer.

Edit I just see, for high n the interpolation generates artificial peaks for f[m] at m near n, which give maxima above 1.0. Have no more time to play with InterpolationOrder or else to correct it. (As rp[a,n] shows,you can set it to 1 this case.) Leave it to you.

• In this problem this solution works very well. Comparing its solutions on some points evaluated in my original formulation I get an error lower than 10^-4. Before accepting the solution, could you elaborate more on the role of Plot in rp? Does it just a sampling of the function? Also I had to change Manipulate in Manipulate[{pl2 = Plot[nint[m, a, n], {m, 0, n}, PlotRange -> {.45, 1.05}, ImageSize -> 400], {"Max=", Max[Cases[pl2, Line[data_] :> data, Infinity][[1, All, 2]]]}}, {{a,1}, 0, 2, Appearance -> "Labeled"}, {{n, 3}, 1, 20, Appearance -> "Labeled"}, ContinuousAction -> False]` Dec 10, 2020 at 11:12
• Plot does just generate data points for the allowed m. The higher MaxRecursion, the more points, where funtion changes fast. Reap collects these points for later generating an interpolation function. Then maximizing an InterpolationFunction is an easy and fast job. Plot calls nint about 30 plus lets say 10 times, whereas NMaximize would call nint some hundred times to get closer to max. You can make error lower with higher PlotPoints and higher MaxRecursion, but this takes more time. I am using older version 8.0, that's why you had to change routine to get the sampled points from pl2. Dec 10, 2020 at 13:37