When p ranges from 0 to infinity, S is obtained by numerical integration

Question

The parameters and functions are loaded as follows:

Clear["Global`*"]
epslon = 8.9;
Z = 1.0/epslon; l = 1.0; R0 = 2;
EE = 2.419601;
β0 = Sqrt[-Z^2/(2*EE)];
Rnl[r_] := Hypergeometric1F1[l + 1 - β0, 2 l + 2, (2 Z r)/β0] Exp[(-2 Z r)/(2 β0)] ( (2 Z r)/β0)^l;
Rnlh[p_?NumericQ] := -Sqrt[(2/Pi)] NIntegrate[r^2 Rnl[r] (Sin[p r]/(p r)^2 - Cos[p r]/(p r)), {r, 0, R0}];
normp = Norm[Sqrt[NIntegrate[(Rnlh[p]^2) p^2, {p, 0, Infinity}]]];
Rnlend[p_?NumericQ] := Rnlh[p]/normp;
fp[p_?NumericQ, θp_?NumericQ] := (3 Cos[θp]^2)/(4 π) Norm[Rnlend[p]]^2;

What I want to calculate is: $S=-2\pi \int_0^{\infty}{\int_0^{\pi}{fp\cdot \ln \left( fp \right) \cdot p^2}}\sin \theta dpd\theta$

Because the range of p is 0 to infinity when calculating S, no result can be obtained after the calculation for a long time. I give the range of p as 0-10 for preliminary test.

Method 1：

S1 = -2 Pi  NIntegrate[fp[p, θ] Log[fp[p, θ]] p^2 Sin[θ], {p, 0, 10}, {θ, 0, Pi}] // AbsoluteTiming

{1155.65, 4.42307}

Method 2： the idea of changing variables

$$ fp=\frac{3}{4\pi}\cos ^2\theta \cdot \left| \text{R}nlend \right|^2 $$ $$ S=-2\pi \int_0^{\infty}{\int_0^{\pi}{fp\cdot \ln \left( fp \right) \cdot p^2}}\sin \theta dpd\theta $$ $$ Let\ \cos \theta =y $$ $$ S2=S=-2\pi \int_0^{\infty}{\int_{-1}^1{\frac{3y^2}{4\pi}}}\left| Rnlend \right|^2\ln \left( \frac{3}{4\pi}y^2\left| Rnlend \right|^2 \right) p^2dpdy $$

S2 = -2 Pi  NIntegrate[(3 y^2)/(4 Pi)
      Norm[Rnlend[p]]^2 Log[3/(4 Pi) y^2 Norm[Rnlend[p]]^2] p^2, {p, 
     0, 10}, {y, -1, 1}] // AbsoluteTiming

{234.691, 4.42307}

To sum up, we can find that the speed of Method 2 is much faster than that of Method 1, but it is still unable to solve the solution when the p range is 0 to infinity. So ask if there is a better solution.

Answer 1

We can try to reduce the double integral to a single integral:

$$ \begin{array}{l} fp=\frac{3\cos ^2\theta}{4\pi}\left| Rnlend \right|^2\\ S=-2\pi \int_0^{\infty}{\int_0^{\pi}{fp\cdot \ln \left( fp \right) p^2\sin \theta dpd}}\theta\\ =\int_0^{\infty}{\int_0^{\pi}{\frac{3\cos ^2\theta}{4\pi}\left| Rnlend \right|^2\cdot \left( \ln \left( \frac{3\cos ^2\theta}{4\pi} \right) +\ln \left( \left| Rnlend \right|^2 \right) \right) p^2\sin \theta dpd}}\theta\\ =\int_0^{\infty}{\int_0^{\pi}{\frac{3\cos ^2\theta}{4\pi}\left| Rnlend \right|^2\ln \left( \frac{3\cos ^2\theta}{4\pi} \right) p^2\sin \theta dpd}}\theta +\int_0^{\infty}{\int_0^{\pi}{\frac{3\cos ^2\theta}{4\pi}\left| Rnlend \right|^2\ln \left( \left| Rnlend \right|^2 \right) p^2\sin \theta dpd}}\theta\\ \end{array} $$

Let: $$ \begin{array}{l} C_1=\int_0^{\pi}{\frac{3\cos ^2\theta}{4\pi}\ln \left( \frac{3\cos ^2\theta}{4\pi} \right) \sin \theta}d\theta\\ C_2=\int_0^{\pi}{\frac{3\cos ^2\theta}{4\pi}\sin \theta}d\theta\\ \end{array} $$

So the final simplification result is as follows: $$ S=-2\pi \left( C_1\int_0^{\infty}{\left| Rnlend \right|^2p^2dp}+C_2\int_0^{\infty}{\left| Rnlend \right|^2\ln \left( \left| Rnlend \right|^2 \right) p^2dp} \right) $$

Comparing the above two methods, when p ranges from 0 to 10:

Clear["Global`*"]
epslon = 8.9;
Z = 1.0/epslon; l = 1.0; R0 = 2;
EE = 2.419601;
\[Beta]0 = Sqrt[-Z^2/(2*EE)];
Rnl[r_] := 
  Hypergeometric1F1[l + 1 - \[Beta]0, 
    2 l + 2, (2 Z r)/\[Beta]0] Exp[(-2 Z r)/(2 \[Beta]0)] ((2 Z r)/\
\[Beta]0)^l;
Rnlh[p_?NumericQ] := -Sqrt[(2/Pi)] NIntegrate[
    r^2 Rnl[r] (Sin[p r]/(p r)^2 - Cos[p r]/(p r)), {r, 0, R0}];
normp = Norm[Sqrt[NIntegrate[(Rnlh[p]^2) p^2, {p, 0, Infinity}]]];
Rnlend[p_?NumericQ] := Rnlh[p]/normp;
C1 = NIntegrate[(3 Cos[\[Theta]]^2)/(4 Pi)
     Log[(3 Cos[\[Theta]]^2)/(4 Pi)] Sin[\[Theta]], {\[Theta], 0, Pi}];
C2 = NIntegrate[(3 Cos[\[Theta]]^2)/(4 Pi) Sin[\[Theta]], {\[Theta], 
    0, Pi}];
S = -2 Pi (C1 NIntegrate[Norm[Rnlend[p]]^2 p^2, {p, 0, 10}] + 
     C2 NIntegrate[
       Norm[Rnlend[p]]^2 Log[Norm[Rnlend[p]]^2] p^2, {p, 0, 10}]) // 
  AbsoluteTiming

The result:

{47.8886, 4.42307}

When p ranges from 0 to infinity:

S = -2 Pi (C1 NIntegrate[Norm[Rnlend[p]]^2 p^2, {p, 0, Infinity}] + 
     C2 NIntegrate[
       Norm[Rnlend[p]]^2 Log[Norm[Rnlend[p]]^2] p^2, {p, 0, 
        Infinity}]) // AbsoluteTiming

The result:

{780.926, 4.43211}

This method is optimized to meet my needs in terms of computational efficiency and accuracy.