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

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.

• normp is complex? Commented Feb 24, 2023 at 11:17
• The actual normalization coefficient normp should be real. I've added Norm to fp. normp can also be Norm in advance. Commented Feb 24, 2023 at 11:29
• But normp is still used inside Rlend[] Commented Feb 24, 2023 at 11:33
• Thanks, I have updated the code. Commented Feb 24, 2023 at 11:39
• Do you have more information about your parameters epslon, Z, l, R0, EE? What values can they take on? There may be some preprocessing / simplifying that can be done in order to significantly accelerate the computation(s). Commented Feb 28, 2023 at 14:21

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.