# Exact solution of the integral

Can the following integral be solved exactly?

ClearAll["Global*"]

Psi[r_, z_, b_, a_] = Exp[-b*z^2]*Exp[-a*r^2];

VB[r_, z_] = 1/(2*Sqrt[r^2 + z^2])*Exp[-3/10*Sqrt[r^2 + z^2]];

PB[b1_, a1_, b2_, a2_] =
Integrate[
Psi[r, z, b2, a2]*VB[r, z]*Psi[r, z, b1, a1]*r, {r,
0, \[Infinity]}, {z, -Infinity, Infinity},
Assumptions -> {a1 > 0, b1 > 0, a2 > 0, b2 > 0}]


• You should use 3/10 rather than 0.3 if you're desiring an analytical solution. And you might consider special cases such as a1+a2 = b1+b2 = some integer and look for patterns.
– JimB
May 29 at 22:01
• If you look at the special case with n = some positive integer, then the integral will be 1/(2 n) - (3 E^((9/400)/n) Sqrt[\[Pi]] Erfc[3/(20 Sqrt[n])])/(40 n^(3/2)). That's simply from trying positive integers 1 through 6 and figuring out the pattern. I suspect that this formula will also work for any real values of n > 0.
– JimB
May 29 at 22:37
• It turns out that Integrate works directly with the above special case and resorting to looking for patterns wasn't necessary this time.
– JimB
May 29 at 23:09

The integration over r can be done, but not the integration over z. You can see this by splitting the integration

ClearAll["Global*"]
Psi[r_, z_, b_, a_] = Exp[-b*z^2]*Exp[-a*r^2];
VB[r_, z_] = 1/(2*Sqrt[r^2 + z^2])*Exp[-3/10*Sqrt[r^2 + z^2]];
integrand = Psi[r, z, b2, a2]*VB[r, z]*Psi[r, z, b1, a1]*r And now

int1 = Integrate[integrand, {r, 0, ∞},
GenerateConditions -> False,
Assumptions -> {a1 > 0, b1 > 0, a2 > 0, b2 > 0}] However, it can't do the above over z, either definite or indefinite

 Integrate[int1,z,Assumptions->{a1>0,b1>0,a2>0,b2>0},GenerateConditions->False] Integrate[int1,{z,-Infinity,Infinity},Assumptions->{a1>0,b1>0,a2>0,b2>0},GenerateConditions->False]


Same result. Also Rubi can't.

Switching the order over the integration does not help. i.e integration over z first produces no result. So the above is the best you can get for now for analytical result, at least using direct Mathematica Integrate command.

ps. use exact numbers with Integrate. Changed your .3 to 3/10

V 13.2.1

• Thanks for the explanation! May 29 at 22:17
• What do you mean by using a direct Integrate command? Could there be any other ways? May 29 at 22:20
• @MamMam I just meant using the command as is. It might help Integrate if you use some actual numerical values for these parameters you have. The Ercf inside the integrate command makes look like there will not be easy way to integrate this. May 29 at 22:25

If the special case of $$a_1+a_2=b_1+b_2=\text{some positive constant }n$$, then

integrand = (Psi[r, z, b2, a2]*VB[r, z]*Psi[r, z, b1, a1]*r // FullSimplify) /. a1 + a2 -> n /. b1 + b2 -> n
2 Integrate[integrand, {r, 0, ∞}, {z, 0, ∞}, Assumptions -> n > 0]


results in

(* 1/(2 n) - (3 E^((9/400)/n) Sqrt[π] Erfc[3/(20 Sqrt[n])])/(40 n^(3/2)) *)


There are other special cases such as $$a_1+a_2=0$$ and $$b_1+b_2=n$$ with $$n>0$$:

2 Integrate[(E^(-n z^2 - 3/10 Sqrt[r^2 + z^2]) r)/(2 Sqrt[r^2 + z^2]), {r, 0, ∞}, {z, 0, ∞},
Assumptions -> n > 0]

(* (5 E^((9/400)/n) Sqrt[π] Erfc[3/(20 Sqrt[n])])/(3 Sqrt[n]) *)


The basic integral of Exp and Erf cannot be easily obtained . Its result contains the little known OwenT function (look up the Mathematica documentation for this) :

Integrate[
Erf[a*x + b]/E^(p^2*x^2), {x,0,Infinity}] =
(Sqrt[Pi]/p)*((1/2)*Erf[(b*p)/Sqrt[a^2 + p^2]] +
2*OwenT[(Sqrt*b*p)/Sqrt[a^2 + p^2], a/p]).


From here you get contracting the a's and b's together with using elementary algebra to the result of Nasser' s integral with Exp and Erfc as :

(E^(9/400/a1)*Pi*Erfc[(3*Sqrt[-a1 + b1])/(20*Sqrt[a1]*Sqrt[b1])])/(4*
Sqrt[a1]*Sqrt[-a1 + b1])


Numeric check for b1 > a1 :

a1 = 0.3; b1 = 0.8; {NIntegrate[(Exp[9/(400*a1) + (a1 - b1)*z^2]*
Sqrt[Pi]*Erfc[(3 + 20*a1*z)/(20*Sqrt[a1])])/(4*
Sqrt[a1]), {z, -Infinity, Infinity}],
(E^(9/400/a1)*Pi*Erfc[(3*Sqrt[-a1 + b1])/(20*Sqrt[a1]*Sqrt[b1])])/(4*
Sqrt[a1]*Sqrt[-a1 + b1])}
(*  {1.66006, 1.66006}  *)


a nice simple form.

Edit:

an easier way to calculate

int = Integrate[E^(-a z^2)  Erfc[b + c z], {z, -Infinity, Infinity}]


(which resembles the OP's integral after integrating over r) is by differentiation under the integral sign :

Derivative with respect to b :

D[E^(-a z^2)  Erfc[b + c z], b]
(* -((2*E^((-a)*z^2 - (b + c*z)^2))/Sqrt[Pi]) *)


Now integrate over z:

Integrate[-((2*E^((-a)*z^2 - (b + c*z)^2))/Sqrt[Pi]),
{z, -Infinity, Infinity},
GenerateConditions -> False]
(* -(2/(E^((a*b^2)/(a + c^2))*Sqrt[a + c^2])) *)


Find the Antiderivative with respect to b to 'cancel' the derivative :

Integrate[-(2/(E^((a*b^2)/(a + c^2))*Sqrt[a + c^2])), b]
(* -((Sqrt[Pi]*Erf[(Sqrt[a]*b)/Sqrt[a + c^2]])/Sqrt[a]) *)


To get the missing constant of integration compare both expressions at b = 0 :

Integrate[Erfc[c*z]/E^(a*z^2), {z, -Infinity, Infinity},
GenerateConditions -> False]
(* Sqrt[Pi]/Sqrt[a] *)


and the final result is

int = (Sqrt[Pi]/Sqrt[a])*Erfc[(Sqrt[a]*b)/Sqrt[a + c^2]]


Check by plot:

a = 0.7; c = 0.3; Plot[{NIntegrate[
Erfc[b + c*z]/E^(a*z^2), {z, -Infinity, Infinity}],
(Sqrt[Pi]*Erfc[(Sqrt[a]*b)/Sqrt[a + c^2]])/Sqrt[a]}, {b, 0, 1},
PlotStyle -> {Blue, Dashed}] This looks suspiciously like a time-independent quantum mechanics problem, solving for the expectation value of VB. Forgive me, if I am mistaken. If it is, in fact, a quantum mechanics problem, you are approaching the problem incorrectly, and the second psi should be the complex conjugate transpose of the first psi. Further, you are only integrating over two dimensions, which would yield a nonsensical result. Perhaps converting from cylindrical to spherical coordinates would help.

• Welcome to Mathematica StackExchange! Your post does not really seem like an answer but more as a comment, which are conventionally placed directly under the question (red link "Add a comment") :) May 30 at 16:39
• The site wouldn't let me comment, or I would have... In any case, this is the information the OP needs, if I am correct, and he is trying to solve a quantum mechanics problem which was not formulated correctly. May 30 at 16:47
• Your comment is good and where expertise outside of Mathematica programming is a big plus: not just about getting something to work but if one should be doing that in the first place. I suspect a moderator can change your answer to a comment.
– JimB
May 30 at 17:51
• @ndtech, thanks for the comment! This is the question about mathematics, not about physics. However, add in 2*Pi as a multiplier and we get integration over the entire volume in cylindrical coordinates. May 30 at 20:10
• @MamMam This would be true if cylindrical coordinates were useful in quantum mechanics. We generally use spherical coordinates. That would be a good point in electrodynamics. Since you’re obviously a physics student (in my experience, mathematicians don’t tend to know that sort of thing off of the top of their heads), I find it hard to believe this is not a quantum mechanics problem. May 31 at 18:19