# Get a numerical solution of an equation containing an integral?

Posted 5 months ago
672 Views
|
2 Replies
|
2 Total Likes
|
 I have the following problem: given r, I have to find the {x,y} which satisfies: Integrate[1/(E^((2*x^2-4*x*d+2*d^2+2*y^2)/(8*d+4*r^2))*(Sqrt[d]*(2*d+r^2))),{d,0,Infinity}]==1 however, of all the solutions, I am interested to find the one which has the maximum y. Unfortunately, the integral can be solved only numerically. Any help? Bye
2 Replies
Sort By:
Posted 5 months ago
 If you restrict x, say to be nonnegative, then this is a straightforward optimization problem once the integral is made into a "black box" function that only evaluates for numeric input. f[x_?NumberQ, y_?NumberQ, r_?NumberQ] := NIntegrate[ 1/(E^((2*x^2 - 4*x*d + 2*d^2 + 2*y^2)/(8*d + 4*r^2))*(Sqrt[ d]*(2*d + r^2))), {d, 0, Infinity}] - 1 Example: FindMaximum[{y, f[x, y, 1/2] == 0, x >= 0}, {x, y}] (* Out[69]= {1.26438464527, {x -> 0.340246814653, y -> 1.26438464527}} *) 
 Dear Daniel Lichtblau, thank you very much for your double help. I can not restrict x to be non negative. However, thank to your help, I can now solve the complete problem, that is: ranger = {0, 0.001, 0.01, 0.1, 1, 10}; rangeRy = {10^-6, 10^-4, 10^-2, 10^0, 10^2, 10^4, 10^6}; f2[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}]; solu2 = Table[FindMaximum[{y, Ry/Sqrt[\[Pi]] f2[x, y, 0, r] == 1, y > 0}, {x, y}], {r, ranger}, {Ry, rangeRy}] As you see, I need to solve it for different r and Ry. Unfortunately, the numerical computation does not seem to converge over the all ranges defined by me, although it easily converge for Ry close to 1, say between 10^-2 - 10^2, and very small values of r, around 0. You also advised me to reformulate the problem as: f1[x_, y_, z_, r_] := 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[NMaximize[{y, Ry/Sqrt[\[Pi]] f1[x, y, 0, r] == 1, y > 0}, {x, y}], {r, ranger}, {Ry, rangeRy}] but this version has the same problems, even if it seems to converge quicker.Finally I used the following calculations to check the convergence of the solutions, but often I do not have good results. MatrixForm[Transpose[Table[rangeRy[[j]]/Sqrt[\[Pi]] f1[solu1[[i, j, 2, 1, 2]], solu1[[i, j, 2, 2, 2]], 0, ranger[[i]]] - 1, {i, 1, Length[ranger]}, {j, 1, Length[rangeRy]}]]] MatrixForm[Transpose[Table[rangeRy[[j]]/Sqrt[\[Pi]] f2[solu2[[i, j, 2, 1, 2]], solu2[[i, j, 2, 2, 2]], 0, ranger[[i]]] - 1, {i, 1, Length[ranger]}, {j, 1, Length[rangeRy]}]]] Do I need to increase somehow the accuracy of the FindMaximum or of the NMaximize?Kind regards Umberto Attachments: