1
$\begingroup$

Short description at the beginning. I have a function of two variables:

dσ[x_, γ_] := ((-1 + γ)^2 +  4 x^4 γ^2 (1 + γ)^2 - 2 x^2 (-1 + γ (1 + 2 γ (2 + γ))))/(
x^3 (-1 + γ)^2 (1 + γ)^(3/2) Sqrt[1 - 2/x^2 + γ])

One can see (and also check with FunctionDomain) that this functin is real for γ > 2/x^2 - 1. And I need to integrate this with some Gaussian distribution:

f1[x_?NumericQ] := Block[{γav = 10^5, σ = 100}, NIntegrate[Exp[-(γ - γav)^2/(2 σ^2)]/(
Sqrt[2 π] σ) dσ[x, γ], {γ,2/x^2 - 1, γav + 20 σ}]]

Unfortunately, this integral can't be solved symbolically, so I make numeric integration. I've substituted for some values of γ av and σ. The lower limit of integral is because of restriction on γ, upper limit is arbitrary -- one can use 10, or 100 σ, NIntegrate can't deal with . Now I can plot this function vs. x and see this has peak near x0 = Sqrt[2/(γav + 1)]. My final goal is to build a function describing the dependence of this peak's distance form x0 as function of γ and σ. So, I think I need some kind of interpolation for f1[x].

I've tried f2 = FunctionInterpolation[f1[x],{x, x0 - 10^-4, x0 + 10^-4}]. This gives many warnings about NIntegrate and FunctionInterpolation, but the plot of interpolation is similar to the plot of f1. Then FindMaximum[f2[x], {x, x0}] gives strange messages about extrapolation -- starting value x0 is inside of interpolation region and very near to maximum, why MMA goes outside region of interpolation?

NMaximize and FindMaximum on f1[x] gives many warnings about 'function has complex value near ...' and no output.

How to make robust procedure of finding maximum for arbitrary γav and σ?

$\endgroup$
0
$\begingroup$

To overcome all these difficulties, there is a simple graphical solution:

dσ[x_, γ_] = ((-1 + γ)^2 + 4 x^4 γ^2 (1 + γ)^2 - 2 x^2 (-1 + γ (1 + 2 γ (2 +γ))))/(x^3 (-1 + γ)^2 (1 + γ)^(3/2) Sqrt[1 - 2/x^2 + γ]);
f1[x_?NumericQ] := Block[{γav = 1/2*10^5, σ = 50}, 
  NIntegrate[Exp[-(γ - γav)^2/(2 σ^2)]/(Sqrt[2 π] σ) dσ[x, γ], {γ,2/x^2 - 1, γav + 100 σ}]]

p = Plot[f1[x], {x, 0.0063, 0.0064}, PlotRange -> All, GridLines -> Automatic]//Quiet

enter image description here

pnt = First@Cases[Normal@p, Line[pt__] :> pt, ∞];
max = MaximalBy[pnt, Last]

{{0.00632692, 0.45993}}

A check with Table:

tab=ParallelTable[{x, f1[x]}, {x, 0.0063, 0.0064, 10^-7}];
MaximalBy[tab, Last]

{{0.0063269, 0.459927}}

ListLinePlot[tab, GridLines -> Automatic, Epilog -> {Red, PointSize@Medium, Point@max}]

enter image description here

"... why MMA goes outside region of interpolation?" , I don't know. Perhaps a developer can answer.

Edit

Because of your objection, I have recalculated the system with γav=1/2*10^5, σ=50. Where is the problem?

dσ[x_, γ_] = ((-1 + γ)^2 + 4 x^4 γ^2 (1 + γ)^2 - 2 x^2 (-1 + γ (1 + 2 γ (2 +γ))))/(x^3 (-1 + γ)^2 (1 + γ)^(3/2) Sqrt[1 - 2/x^2 + γ]);
f1[x_?NumericQ] := Block[{γav = 1/2*10^5, σ = 50}, 
  NIntegrate[Exp[-(γ - γav)^2/(2 σ^2)]/(Sqrt[2 π] σ) dσ[x, γ], {γ,2/x^2 - 1, γav + 20 σ}]]

p = Plot[f1[x], {x, 0.0063, 0.0064}, PlotRange -> All, GridLines -> Automatic] // Quiet

enter image description here

pnt = First@Cases[Normal@p, Line[pt__] :> pt, \[Infinity]];
max = MaximalBy[pnt, Last]

{{0.00632692, 0.45993}}

tab = Table[{x, f1[x]}, {x, 0.0063, 0.0064, 10^-7}];
MaximalBy[tab, Last]

{{0.0063269, 0.459927}}

ListLinePlot[tab, GridLines -> Automatic, Epilog -> {Red, PointSize@Medium, Point@max}]

enter image description here

$\endgroup$
  • $\begingroup$ thanks for reply. But try use different values for γav and σ, say γav = 0.5 10^5, σ = 50. Plot method gives {{0.00632688, 0.459918}}. Table method Table[{x, f1[x]}, {x, 0.006, 0.0065, 10^-7}] with many warnings finally reads {{0.0063269, 0.459927}}. Table with smaller step Table[{x, f1[x]}, {x, 0.006, 0.0065, 10^-8}] takes too much time, and finally can't compute numerical integral and gives no answer. As I wrote I need some model for max vs. γav, σ, so I have to compute for many different values. $\endgroup$ – Alx Feb 12 '17 at 0:32
  • $\begingroup$ Continue comment. Plot and Table methods give not exactly same results. I read somewhere here in MSE, that such Plot method can be used as an estimation (Plot tries to show 'nice' picture, probably not too accurate in details of calculation of points). More accurate are NMaximize and the like, or NDSolve with original or interpolated function. We see, that Table (and interpolations as well) can not work for some parameters. That is the problem. $\endgroup$ – Alx Feb 12 '17 at 0:41
  • $\begingroup$ @Alx Sorry, please contact me if you found a successful method. This method provides a result, and thinks it is a graphical method, not an exact one! If you can plot your function, you can also determine the extremes. The error messages are generated during the calculation of your functions. They need a thorough examination. $\endgroup$ – user36273 Feb 12 '17 at 9:32
  • $\begingroup$ concerning your Edit. I'm not sure if 20 σ is enough in upper limit of integral. Strictly speaking, we should integrate up to infinity. So, if I change that 20 σ to, say, 50 σ I get odd messages from MMA: "... loss of precision ...", and "γ = 50000 + 50. σ is not a valid limit of integration". Why? This "Plot" method is good in the sense of it uses adaptive points' sampling and hence is faster than "Table" method, in which we have to set small enough step size for not to miss the peak. But if we think about infinity limit all goes slow and complicate. $\endgroup$ – Alx Feb 12 '17 at 11:33
  • $\begingroup$ @Alx It is not easy in life, but you can not expect anything. NIntegrate also reaches its limits at some point. Try some Methods and/or decrease the Accuray to 5. I have edited the first example in my answer with γav=1/2*10^5,σ=50 and γav+100 σ. I also used ParallelTable. With Mathematica 11.01 I have no problems. Quiet the warnings if they occur. This method yields at least a result; with NMaximize or FindMaximum you get nothing! $\endgroup$ – user36273 Feb 12 '17 at 14:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.