Сan this integral be solved exactly (using Integrate)?

Question

Is it possible to solve exactly the following integral?

Psi[r_, n_] := (2 E^(-(r/n)) Sqrt[n!/(-1 + n)!] Hypergeometric1F1[
      1 - n, 2, (2 r)/n])/n^2;
Px2[n1_, n2_] = 
 Integrate[
  Psi[r, n2](-(1/(3 r))Exp[-5/2*r])Psi[r, n1]r^2, {r, 0, 
   Infinity}]

n1, n2 > 0 (n1, n2 = 1, 2, 3, 4 ...)

Answer 1

Converting Hypergeometric1F1 to LaguerreL function we have: $$\int_0^{\infty } -\frac{4 e^{-\left(\left(\frac{5}{2}+\frac{1}{\text{n1}}+\frac{1}{\text{n2}}\right) r\right)} r \, _1F_1\left(1-\text{n1};2;\frac{2 r}{\text{n1}}\right) \, _1F_1\left(1-\text{n2};2;\frac{2 r}{\text{n2}}\right)}{3 \text{n1}^{3/2} \text{n2}^{3/2}} \, dr=\\\int_0^{\infty } -\frac{4 r \exp \left(-\left(\frac{5}{2}+\frac{1}{\text{n1}}+\frac{1}{\text{n2}}\right) r\right) L_{-1+\text{n1}}^1\left(\frac{2 r}{\text{n1}}\right) L_{-1+\text{n2}}^1\left(\frac{2 r}{\text{n2}}\right)}{3 \text{n1}^{5/2} \text{n2}^{5/2}} \, dr$$ and using formula from book: Integrals and Series Volume 2,Special Functions- A. Prudnikov 1986 year on page 477 formula 6 we have:

$$\int_0^{\infty } r^{\lambda } \exp (-p r) L_m^{\lambda }(b r) L_n^{\lambda }(c r) \, dr=\frac{\left(\Gamma (m+n+\lambda +1) (p-b)^m (p-c)^n\right) \, _2F_1\left(-m,-n;-m-n-\lambda ;\frac{p (p-b-c)}{(p-b) (p-c)}\right)}{m! n! p^{m+n+\lambda +1}}$$

Substituting parameters to formula from book:

ANSWER=-4/(3*n1^(5/2)*n2^(5/2))*(Gamma[m + n + \[Lambda] + 1]*(p - b)^m*(p - c)^n)/(
       m! n! p^(m + n + \[Lambda] + 1))*
       Hypergeometric2F1[-m, -n, -m - n - \[Lambda], (
       p*(p - b - c))/((p - b)*(p - c))] /. \[Lambda] -> 1 /. 
       p -> 5/2 + 1/n1 + 1/n2 /. m -> -1 + n1 /. b -> 2/n1 /. 
       n -> -1 + n2 /. c -> 2/n2 // FullSimplify;

we have the Answer:

$$\int_0^{\infty } -\frac{4 e^{-\left(\left(\frac{5}{2}+\frac{1}{\text{n1}}+\frac{1}{\text{n2}}\right) r\right)} r \, _1F_1\left(1-\text{n1};2;\frac{2 r}{\text{n1}}\right) \, _1F_1\left(1-\text{n2};2;\frac{2 r}{\text{n2}}\right)}{3 \text{n1}^{3/2} \text{n2}^{3/2}} \, dr=-\frac{4 \left(\frac{5}{2}+\frac{1}{\text{n1}}-\frac{1}{\text{n2}}\right)^{-1+\text{n2}} \left(\frac{5}{2}-\frac{1}{\text{n1}}+\frac{1}{\text{n2}}\right)^{-1+\text{n1}} \left(\frac{5}{2}+\frac{1}{\text{n1}}+\frac{1}{\text{n2}}\right)^{-\text{n1}-\text{n2}} \Gamma (\text{n1}+\text{n2}) \, _2F_1\left(1-\text{n1},1-\text{n2};1-\text{n1}-\text{n2};1-\frac{8}{5 (2 \text{n2}+\text{n1} (-2+5 \text{n2}))}-\frac{8}{5 (-2 \text{n2}+\text{n1} (2+5 \text{n2}))}\right)}{3 \text{n1}^{5/2} \text{n2}^{5/2} \Gamma (\text{n1}) \Gamma (\text{n2})}$$

HoldForm[Integrate[-(1/(3 n1^(3/2) n2^(3/2)))
  4 E^-((5/2 + 1/n1 + 1/n2) r) r Hypergeometric1F1[1 - n1, 2, (
 2 r)/n1] Hypergeometric1F1[1 - n2, 2, (2 r)/n2], {r, 0, 
Infinity}] == -((4 (5/2 + 1/n1 - 1/n2)^(-1 + n2) (5/2 - 1/n1 + 1/n2)^(-1 + 
  n1) (5/2 + 1/n1 + 1/n2)^(-n1 - n2)
  Gamma[n1 + n2] Hypergeometric2F1[1 - n1, 1 - n2, 1 - n1 - n2, 
  1 - 8/(5 (2 n2 + n1 (-2 + 5 n2))) - 8/(
  5 (-2 n2 + n1 (2 + 5 n2)))])/(
  3 n1^(5/2) n2^(5/2) Gamma[n1] Gamma[n2]))]

Verification:

ANSWER /. n1 -> 1 /. n2 -> 3
(-(5776/(279841 Sqrt[3])))

$$-\frac{5776}{279841 \sqrt{3}}$$

Integrate[-(1/(3 n1^(3/2) n2^(3/2)))
  4 E^-((5/2 + 1/n1 + 1/n2) r) r Hypergeometric1F1[1 - n1, 2, (
 2 r)/n1] Hypergeometric1F1[1 - n2, 2, (2 r)/n2] /. n1 -> 1 /. 
 n2 -> 3, {r, 0, Infinity}]
(-(5776/(279841 Sqrt[3])))

$$-\frac{5776}{279841 \sqrt{3}}$$

Answer 2

One can adapt the code in this answer:

f[a_,n_]=Integrate[r^n*Exp[-a*r],{r,0,Infinity},Assumptions->{a>0,n>0}];
coeff[n_]:=2*Sqrt[n!/(-1+n)!]/n^2;
Px[n1_,n2_]:=Px[n1,n2]=If[n1>n2,Px[n2,n1],
  coeff[n1]coeff[n2](-1/3)Total[Map[f[5/2+1/n1+1/n2,#[[1,1]]+1]#[[2]]&,
  CoefficientRules[Hypergeometric1F1[1-n1,2,2*r/n1]Hypergeometric1F1[1-n2,2,2r/n2],r]]]];

Timing:

AbsoluteTiming[Table[Px[n1,n2],{n1,1,50},{n2,1,50}];]
(* about 1.5 seconds *)

Answer 3

The answer below uses GeneratingFunction to compute the integral in terms of a DifferenceRoot (a linear recursion formula). Defining a recursive function for the linear recursion one obtains significant speed at large values of n1 or n2.

---

The integrand:

Psi[r_, n_] := (2 E^(-(r/n)) Sqrt[n!/(-1 + n)!] Hypergeometric1F1[
      1 - n, 2, (2 r)/n])/n^2;
int = Psi[r, n2](-(1/(3 r))Exp[-5/2*r])Psi[r, n1]r^2 // 
  FullSimplify

$$-\frac{4 r \, _1F_1\left(1-\text{n1};2;\frac{2 r}{\text{n1}}\right) \, _1F_1\left(1-\text{n2};2;\frac{2 r}{\text{n2}}\right) e^{-\left(r \left(\frac{1}{\text{n1}}+\frac{1}{\text{n2}}+\frac{5}{2}\right)\right)}}{3 \text{n1}^{3/2} \text{n2}^{3/2}}$$

---

We can compute the integral using GeneratingFunction

First notice :

GeneratingFunction[Hypergeometric1F1[1 - n, 2, x], n, 
  u] // FunctionExpand

$$-\frac{e^x \left(e^{\frac{x}{u-1}}-1\right)}{x}$$

We can then replace $n$ with $n1$ and $x$ with $2r/(n1)$ in the first product and a similar replacement in the second product.

We apply a generating function for each factor and obtain a function of u1 and u2. After applying Simplifythe generating function of the integrand is:

-((E^((-(5/2) + 1/n1 + 1/n2) r) (-1 + E^(-((2 r)/(n1 - n1 u1)))) (-1 +
     E^(-((2 r)/(n2 - n2 u2)))))/(3 Sqrt[n1] Sqrt[n2] r))

$$-\frac{e^{r \left(\frac{1}{\text{n1}}+\frac{1}{\text{n2}}-\frac{5}{2}\right)} \left(e^{-\frac{2 r}{\text{n1}-\text{n1} \text{u1}}}-1\right) \left(e^{-\frac{2 r}{\text{n2}-\text{n2} \text{u2}}}-1\right)}{3 \sqrt{\text{n1}} \sqrt{\text{n2}} r}$$

Assuming all the hypothesis are satisfied, the integral of the generating function is the generating function of the integral. We use this to revert back to the series coefficients after integrating:

(
Integrate[-((E^((-(5/2)+1/n1+1/n2) r) (-1+E^(-((2 r)/(n1-n1 u1)))) (-1+E^(-((2 r)/(n2-n2 u2)))))/(3 Sqrt[n1] Sqrt[n2] r)),{r,0,Infinity},GenerateConditions->False]
//SeriesCoefficient[#,{u2,0,n2}]&
//Refine[#,n2>0]&
//SeriesCoefficient[#,{u1,0,n1}]&
//Refine[#,n1>0]&
)

Important note : //Refine[#,n>0]& eliminates the case n=0. Consider checking separately the cases n1=0 and n2=0

The result is :

Verification:

Integrate[int /. n1 -> 1 /. n2 -> 3, {r, 0, Infinity}]

$$-\frac{5776}{279841 \sqrt{3}}$$

$$-\frac{5776}{279841 \sqrt{3}}$$

---

Speeding up the computation

The DifferenceRoot in the result should be reasonably quick but one might do better by extracting the recurrence formula from the DifferenceRoot using something like diffroot[[0,1]][y,x][[1]] and defining a recursive function with memoization. This is what is done below :

(* recurrence formula with y1=y[n-1] and y2=y[n-2] obtained 
 using Solve on the recurrence formula*)
g[n_, n1_, n2_, y1_, 
  y2_] := (n1^2 (-4 + 25 n2^2) (2 (-1 + n) y1 - (-2 + n) y2) - 
    4 n2^2 (2 (-1 + n) y1 + (-2 + n) y2) + 
    4 n1 n2^2 (4 y1 + 5 (-2 + n) y2))/(n (2 n2 + 
      n1 (-2 + 5 n2)) (2 n2 + n1 (2 + 5 n2)))
(* Optional compilation *)
gcompile = 
  Compile[{{n, _Integer}, {n1, _Integer}, {n2, _Integer}, {y1, 

_Real}, {y2, _Real}}, Evaluate[g[n, n1, n2, y1, y2]]];
Clear[recurrence];
recurrence[0, OrderlessPatternSequence[n1_, n2_]] := 
  recurrence[0, n1, n2] = ((-2 n1 + 2 n2 + 5 n1 n2)/(
    2 n1 + 2 n2 + 5 n1 n2))^n2;
recurrence[1, OrderlessPatternSequence[n1_, n2_]] := 
  recurrence[1, n1, n2] = (
   16 n1 n2^2 recurrence[0, n1, n2])/(-4 n1^2 + 4 n2^2 + 20 n1 n2^2 + 
    25 n1^2 n2^2);
recurrence[n_, OrderlessPatternSequence[n1_, n2_]] :=
  recurrence[n, n1, n2] = 
   g[n, n1, n2, recurrence[n - 1, n1, n2], 
    recurrence[n - 2, n1, n2]];
Clear[MemP];
MemP[n1_, n2_] := 
 MemP[n1, n2] = 
  If[n1 > n2, 
   MemP[n2, n1], -(N@recurrence[n1, n1, n2]/(3 Sqrt[n1] n2^(3/2)))] 

Benchmark

Without compilation:

Table[MemP[i, j], {i, 50}, {j, 50}]; // AbsoluteTiming
(* {0.657429, Null} *)

Increasing the maximum range of n1 to 200 and using compilation by switching g to gcompile in the function recurrence above (the output from MemP can be floating points) :

Table[MemP[i, j], {i, 200}, {j, 50}]; // AbsoluteTiming
(* {1.04846, Null} *)

I recall that one should check separately the cases n1=0 and n2=0.