# Numerical Double Integral with symbolic bounds to normalize function

I'm trying to use mathematica to normalize this trial wave function which requires finding the following integral and solving for A so that the total propability equals 1. 'a' is the constant width of the potential well and alpha is the variational parameter. I'm having a lot of dificulty with mathematica accepting non-numeric bounds such as a and while still managing to get numerical results but it also can't be done purely symbolically without encountering a ton of erf functions.

8*A^2*Integrate[((1 - x^2*y^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, {x, y, a}] + 8*A^2*Integrate[((1 - x^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, {x, a, Infinity}]

f1[(x_)?NumericQ, (y_)?NumericQ] := NIntegrate[((1 - x^2*y^2)/a^4 E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {x, y, a}]

NIntegrate[f1[y], {y, 0, a}]

NIntegrate::nlim: y = a is not a valid limit of integration.

N[8*A^2*Integrate[((1 - x^2*y^2)/a^4/
E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a},
{x, y, a}] + 8*A^2*Integrate[
((1 - x^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^
2, {y, 0, a}, {x, a, Infinity}]]

((1/(a^6*\[Alpha]^5))*0.001953125*A^2*
(8.*a^4*\[Alpha]*(-32.*\[Alpha]^2 +
a^4*(3. + 4.*\[Alpha])^2) - 10.026513098524001*
2.718281828459045^(2.*\[Alpha])*a^4*Sqrt[\[Alpha]]*
(-32.*\[Alpha]^2 + 3.*a^4*(3. + 4.*\[Alpha]))*
Erf[1.4142135623730951*Sqrt[\[Alpha]]] +
3.141592653589793*2.718281828459045^
(4.*\[Alpha])*(9.*a^8 - 32.*a^4*\[Alpha]^2 +
256.*\[Alpha]^4)*Erf[1.4142135623730951*
Sqrt[\[Alpha]]]^2))/2.718281828459045^
(4.*\[Alpha]) + ((1/(a^6*\[Alpha]^3))*0.0625*A^2*
Erf[1.4142135623730951*Sqrt[\[Alpha]]]*
(5.0132565492620005*a^2*Sqrt[\[Alpha]]*
(-8.*\[Alpha] + a^2*(3. + 4.*\[Alpha])) +
3.141592653589793*2.718281828459045^
(2.*\[Alpha])*(3.*a^4 - 8.*a^2*\[Alpha] +
16.*\[Alpha]^2)*
Erfc[1.4142135623730951*Sqrt[\[Alpha]]]))/
2.718281828459045^(2.*\[Alpha])

• I mean, you have to make a choice: either use numeric bounds, or accept the fact that A will be in terms of error functions (and what's wrong with error functions?). – march Dec 11 '19 at 5:45

Analytical integration yields good results when applying assumptions.

int = 8*A^2*
Integrate[((1 - x^2*y^2)/a^4/E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y,
0, a}, {x, y, a}, Assumptions -> a > 0 && \[Alpha] > 0] +
8*A^2*Integrate[((1 - x^2)/a^4/
E^(\[Alpha]*((x^2 + y^2)/a^2)))^2, {y, 0, a}, {x, a, Infinity},
Assumptions -> a > 0 && \[Alpha] > 0]


Result is a long expression. Solve for A:

AA[a_, \[Alpha]_] = A /. Solve[1 == int, A]


Get A for a given a = 1:

AA[1, \[Alpha]] // FullSimplify[# \[Alpha] > 0] &


Since A^2 is used, you have two solutions:

Plot3D[Evaluate@AA[a, \[Alpha]], {a, 0, 5}, {\[Alpha], 0, 4},
PlotStyle -> {Blue, Red}]