# Maximize the solution of an equation containing an integral

I have to find {x,y} which makes the integral

Integrate[(1/(E^((x^2 - 2*x*d + d^2 + y^2 + z^2 +
(r^2/d)*z^2)/(2*(d + r^2)))*
(Sqrt[d]*(d + r^2)))), {d, 0, Infinity}]


equal to Pi^0.5/Ry, where Ry is a given constant. Among all the possible solutions, I am interested in the one which maximises y, with the constraint y>0. I have also a good starting point for y. The problem has to be solved for different values of r, say from 0 to 20, and Ry, say from 10^-7 to 10^7. I have set the problem in this way:

ranger = {0,0.001,0.01,0.1,1,10};
rangeRy = {10^-6,10^-4,10^-2,10^-1,10^0,10^1,10^2,10^4,10^6};
f[x_?NumberQ, y_?NumberQ, z_?NumberQ, r_?NumberQ] :=
NIntegrate[(1/(E^((x^2 - 2*x*d + d^2 + y^2 + z^2 +
(r^2/d)*z^2)/(2*(d + r^2)))*
(Sqrt[d]*(d + r^2)))), {d, 0, Infinity}];
solu1 = Table[
FindMaximum[{y, f1[x, y, 0, r] - Sqrt[\[Pi]]/Ry == 0,
y > 0}, {x, {y, Sqrt[(2 Ry)/E]}}], {r, ranger}, {Ry, rangeRy}]


Unfortunately, NIntegrate fails to converge to the solution for all the values of r and Ry. Any help?

• f1 is the same as f? Commented Feb 12, 2019 at 11:48
• yes, it si, it was just a typo.
– umby
Commented Feb 12, 2019 at 16:53

First, it is necessary to limit the variables so that there is no division by 0 and so that the integral converges. Secondly, we must make sure that the problem has a solution. Use ContourPlot to build solution of equation f[x, y, 0, r] - Sqrt[\[Pi]]/Ry == 0:

d0 = 10^-6; r0 = 10^-6; dmax = 20;
ranger = {r0, 0.001, 0.01, 0.1, 1, 10};
rangeRy = {10^-6, 10^-4, 10^-2, 10^-1, 10^0, 10^1, 10^2, 10^4, 10^6};
f[x_?NumberQ, y_?NumberQ, z_?NumberQ, r_?NumberQ] :=
NIntegrate[(1/(E^((x^2 - 2*x*d + d^2 + y^2 +
z^2 + (r^2/d)*z^2)/(2*(d + r^2)))*(Sqrt[
d]*(d + r^2)))), {d, d0, dmax}];
Table[With[{r = ranger[[i]], Ry = rangeRy[[5]]},
ContourPlot[
f[x, y, 0, r] - Sqrt[\[Pi]]/Ry == 0, {x, -1, 2}, {y, 0, 1},
PlotLabel -> Grid[{{"r =", r*1.}, {"Ry =", Ry*1.}}]]], {i, 1,
Length[ranger]}]
Table[With[{r = ranger[[i]], Ry = rangeRy[[6]]},
ContourPlot[
f[x, y, 0, r] - Sqrt[\[Pi]]/Ry == 0, {x, -2, 15}, {y, 0, 4},
PlotLabel -> Grid[{{"r =", r*1.}, {"Ry =", Ry*1.}}]]], {i, 1,
Length[ranger]}]


And so we see that the solution exists for some data. To find the maximum value y, use the code:

With[{r = 1, Ry = 1},
FindMaximum[
y /. FindRoot[
f[x, y, 0, r] - Sqrt[\[Pi]]/Ry == 0, {y, .5}], {x, .5}]] // Quiet
(*{0.537561, {x -> 0.432659}}*)


and

With[{r = 1, Ry = 10},
FindMaximum[
y /. FindRoot[f[x, y, 0, r] - Sqrt[\[Pi]]/Ry == 0, {y, 3}], {x,
5}]] // Quiet

(*{3.27002, {x -> 4.02939}}*)

• Many thanks!! It works fine, even if I change the integral limit to {d,0,Infinity}, as I need. By the way, why you changed the limits to {d,d0,dmax}? However, if I delete Quiet from the FindRoot, Mathematica returns a lot of warning messages although the result is seemed to be correctly computed in the end. Can I really discard all these messages?
– umby
Commented Feb 13, 2019 at 10:14
• Dear Alex Trounev, it would be great if you could help me for two additional questions: 1) I want to express x and y in the ContourPlot in unit of Ry, then on the two axes instead of x and y, I'd like to have x/Ry and y/Ry. 2) it is possible to extract the datapoints of the plotted contour? I mean a matrix containing the x and y, or better the x/Ry and y/Ry, whci represent the contour for each couple of r and Ry?
– umby
Commented Feb 13, 2019 at 10:15
• @umby 1) Do not use d=0, because in the definition of f there is r^2/d;2) Do not use Infinity in numerical calculations; 3) Quiet used to cut off meaningless messages arising during the transfer of parameters; 4) What do you want to use ContourPlot for? Commented Feb 13, 2019 at 12:38
• I want express the x and y of the ContourPlot in units of Ry, this means I do not want to report the value of x or y but x/Ry and y/Ry. The ContourPlot reports the {x,y} which satisfies the condition f[x, y, 0, r] - Sqrt[[Pi]]/Ry == 0, I want a ContourPlot which instead reports x/Ry and y/Ry. It is like a rescaling of the coordinate. Furthermore, I would like to have the couples {x,y} plotted into the ContourPlot as a matrix or a table, as numerical values, let's say.I mean I want to extract the numerical values of the {x,y} plotted in the ContourPlot. Thank you so much...
– umby
Commented Feb 13, 2019 at 15:02