Message Boards Message Boards

GROUPS:

Get a numerical solution of an equation containing an integral?

Posted 2 months ago
397 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

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}} *)
Posted 2 months ago

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:
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract