0
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ f1 is the same as f? $\endgroup$ – Alex Trounev Feb 12 at 11:48
  • $\begingroup$ yes, it si, it was just a typo. $\endgroup$ – umby Feb 12 at 16:53
2
$\begingroup$

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]}]

fig1 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}}*)
$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – umby Feb 13 at 10:14
  • $\begingroup$ 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? $\endgroup$ – umby Feb 13 at 10:15
  • 1
    $\begingroup$ @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? $\endgroup$ – Alex Trounev Feb 13 at 12:38
  • $\begingroup$ 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... $\endgroup$ – umby Feb 13 at 15:02

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.