# Sum of all the n-th numerical evaluation of an integral and its cumulative sum of the square of the n-th value

I wish to compute the expansion coefficient of a wavefunction(i.e. in quantum mechanics) which in itself is an integral, given by $$b_{n00} = \int^{\infty}_{u = 0} \frac{1}{(n)^{3/2}} L^{1}_{n-1}\left(\frac{2u}{n}\right)e^{-\frac{u}{n}}e^{-\frac{u^2}{4}}u^2 du$$

I've been computing it in Mathematica as follows:

Monitor[For[n = 1, n < 10000, n++,Print[N[Integrate[(n^(-3/2))*LaguerreL[n - 1, 1, u/n*2]*Exp[-u/n]*Exp[(-u^2)/4]*u^2, {u, 0, \[Infinity]}]]]], n]


I used Monitor because I let this integral run overnight and I want to know what n-th value it evaluated to. Now in quantum mechanics it is know that the sum of square of the expansion coefficient must equal to one, such that: $$\sum_{n=0}^{\infty}|b_{n00}|^2 = 1$$ and I would like to check that with Mathematica. I tried using NSum(integral^2) in the for loop but I have no idea how to make it work. So how can I write my code so that my loop prints out the n-th value, the evaluated integral value at that particular n, and also the cumulative sum of the square of all the n-th value of the integral at a particular n?

• Are you sure there is no error? Where does the formula for bn00 come from? Nov 8, 2018 at 1:16
• There is definitely an error somewhere. $b_{n00}\sim n^{-1/2}$, and so $\sum^N|b_{n00}|^2\sim \log N$, which diverges. Nov 8, 2018 at 1:52
• It came from a series of intensive derivation steps that I've done myself though I can't guarantee that there's no error, which is why I'm doing all sorts of check now. I've been trying to expand a Gaussian wavepacket as a linear combination of the eigenfunctions of the Hydrogen atom. The original expansion is $\psi(\textbf{r},0) = \sum_{n=1}^{\infty}\sum_{l=0}^{n-1}\sum_{m=-l}^{l} b_{\textit{nlm}}\varphi_\textit{nlm}(\textbf{r})$ Nov 8, 2018 at 9:08
• where $\psi(\textbf{r},0) = \frac{1}{(2\pi)^{3/4}\sigma^{3/2}}\text{exp}\left[ -\frac{(\textbf{r} - \textbf{r}_0)^2}{4\sigma^2} + i\textbf{k}_0 \cdot \textbf{r}\right]$, now I've been solving for the case where the Gaussian wavepacket is stationary at the origin, such that $\textbf{r}_0=\textbf{k}_0 = 0$. Throughout my steps it had been shown that for those conditions, the quantum number $l$ and $m$ have to be 0 otherwise the integral would be zero, hence my integral becomes $b_{n00}$ instead Nov 8, 2018 at 9:13
• @AccidentalFourierTransform but how can you tell? Nov 8, 2018 at 9:16

Let's start from the beginning. Radial wave-function of the hydrogen atom:

\[Psi][n_, l_, r_] :=
Sqrt[ (n - l - 1)!/(n + l)!] E^(-(r/n)) ((2 r)/n)^l 2/
n^2 LaguerreL[n - l - 1, 2 l + 1, (2 r)/n]


Check that these functions are orthonormal.

Table[{n, m,
N[NIntegrate[\[Psi][n, 0, r]*\[Psi][m, 0, r]*r^2, {r, 0, Infinity},
WorkingPrecision -> 20, AccuracyGoal -> 10], 9]}, {n, 1, 5}, {m,
1, 5}]
Out[]= {{{1, 1, 1.00000000}, {1, 2, -3.36248205*10^-17}, {1,
3, -2.20279187*10^-15}, {1, 4, -7.94851825*10^-15}, {1,
5, -4.60327343*10^-15}}, {{2, 1, -3.36248205*10^-17}, {2, 2,
1.00000000}, {2, 3, -2.60214679*10^-17}, {2,
4, -6.53519260*10^-20}, {2, 5, -2.11210560*10^-19}}, {{3,
1, -2.20279187*10^-15}, {3, 2, -2.60214679*10^-17}, {3, 3,
1.00000000}, {3, 4, -9.91435888*10^-20}, {3,
5, -4.97630888*10^-19}}, {{4, 1, -7.94851825*10^-15}, {4,
2, -6.53519260*10^-20}, {4, 3, -9.91435888*10^-20}, {4, 4,
1.00000000}, {4, 5, -7.91052765*10^-20}}, {{5,
1, -4.60327343*10^-15}, {5, 2, -2.11210560*10^-19}, {5,
3, -4.97630888*10^-19}, {5, 4, -7.91052765*10^-20}, {5, 5,
1.00000000}}}


Put C0=(2/Pi)^(1/4); Check that the function C0*Exp[-r^2/4] is normalized to 1

 NIntegrate[(C0*Exp[-r^2/4])^2*r^2, {r, 0, Infinity}]

Out[]= 1.


Define the function

fn00[r_, n_] := C0*Exp[-r^2/4]*\[Psi][n, 0, r]


Draw up a table of coefficients

bn00 =
Table[{n, NIntegrate[fn00[r, n]*r^2, {r, 0, Infinity}]}, {n, 1, 50}]

Out[]= {{1, 0.977625}, {2, 0.0589538}, {3, 0.0173227}, {4,
0.00834252}, {5, 0.00505751}, {6, 0.00348058}, {7, 0.0025892}, {8,
0.00202827}, {9, 0.00164787}, {10, 0.0013754}, {11,
0.00117199}, {12, 0.00101514}, {13, 0.000891033}, {14,
0.000790725}, {15, 0.000708215}, {16, 0.000639328}, {17,
0.000581077}, {18, 0.000531275}, {19, 0.000488284}, {20,
0.000450856}, {21, 0.000418025}, {22, 0.00038903}, {23,
0.000363268}, {24, 0.000340251}, {25, 0.000319585}, {26,
0.000300944}, {27, 0.000284059}, {28, 0.000268706}, {29,
0.000254696}, {30, 0.000241869}, {31, 0.00023009}, {32,
0.000219241}, {33, 0.000209222}, {34, 0.000199948}, {35,
0.000191341}, {36, 0.000183338}, {37, 0.000175879}, {38,
0.000168915}, {39, 0.000162399}, {40, 0.000156294}, {41,
0.000150563}, {42, 0.000145175}, {43, 0.000140101}, {44,
0.000135317}, {45, 0.000130801}, {46, 0.00012653}, {47,
0.000122488}, {48, 0.000118656}, {49, 0.000115021}, {50,
0.000111568}}


Let's make the table of squares of factors

bn002 =
Table[{n,
NIntegrate[fn00[r, n]*r^2, {r, 0, Infinity},
WorkingPrecision -> 20]^2}, {n, 1, 50}]

Out[]= {{1, 0.9557504483674802397}, {2,
0.0034755549726118017898}, {3, 0.00030007673534634702992}, {4,
0.00006959766078190917028}, {5, 0.000025578435774571826133}, {6,
0.000012114445921016939855}, {7, 6.703961625639443923*10^-6}, {8,
4.1138947135635687093*10^-6}, {9, 2.7154720730026987933*10^-6}, {10,
1.8917220176826471449*10^-6}, {11,
1.3735559591569759483*10^-6}, {12,
1.0305155739510539128*10^-6}, {13, 7.939396470412556727*10^-7}, {14,
6.252464840331396825*10^-7}, {15, 5.015688323372287368*10^-7}, {16,
4.0874015210321704819*10^-7}, {17, 3.3765074397536246974*10^-7}, {18,
2.8225324836621724412*10^-7}, {19, 2.3842135219686209084*10^-7}, {20,
2.0327124175675325438*10^-7}, {21, 1.7474451601792588441*10^-7}, {22,
1.5134420591983166023*10^-7}, {23, 1.3196331669341879725*10^-7}, {24,
1.1577086645606575635*10^-7}, {25, 1.0213455739860918387*10^-7}, {26, 9.056731609002574339*10^-8}, {27,
8.068970488644740015*10^-8}, {28, 7.220308241309944513*10^-8}, {29,
6.487016735587959574*10^-8}, {30, 5.850078021393612247*10^-8}, {31,
5.294125813808744772*10^-8}, {32, 4.8066509546483806644*10^-8}, {33,
4.3773988935428924374*10^-8}, {34, 3.9979084192386115555*10^-8}, {35,
3.6611553842324869665*10^-8}, {36, 3.3612752356121418748*10^-8}, {37,
3.0933452373374781182*10^-8}, {38, 2.8532122930031263955*10^-8}, {39,
2.6373558848336880024*10^-8}, {40, 2.4427782599485812935*10^-8}, {41,
2.2669159091800755773*10^-8}, {42, 2.1075677972861880836*10^-8}, {43,
1.9628368559712378482*10^-8}, {44, 1.8310820411034173024*10^-8}, {45,
1.7108788528574499147*10^-8}, {46, 1.6009866723975407308*10^-8}, {47,
1.5003216174614578738*10^-8}, {48, 1.4079338882937737508*10^-8}, {49,
1.3229887842692302985*10^-8}, {50, 1.2447507346560968732*10^-8}}


Calculate the sum of the squares of the coefficients

Total[bn002[[All, 2]]]

Out[]= 0.9596561684347894350


You can check that the balance of the sum from 51 to infinity is about 3.03454*10^-7. Therefore, the sum of the squares bn00 is not equal to 1. This is because the system of functions \[Psi][n, 0, r] is not complete, but only the system \[Psi][n,l,r]is complete.

Maybe you should use List:

bn00 = (N[
Integrate[(u^2*LaguerreL[#1 - 1, 1, (u*2)/#1])/(#1^(3/2)*
E^(u^2/4)*E^(u/#1)), {u, 0, Infinity}]] & )[
Range[1, 3]]

output:
{0.547233, 0.0659997, 0.0290896}


and then:

Total[bn00^2]

• TIP: NIntegrate is much faster than N@Integrate. Nov 8, 2018 at 14:49