Mathematica Asked by Nicola on January 23, 2021
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.
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)
Answered by Ulrich Neumann on January 23, 2021
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.
Answered by Akku14 on January 23, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP