Integration of harmonic function

I am trying to calculate the average distance of electron above the surface. However when I did the Integral, for the first couple of state(n from 1 to around 15) the results are pretty good. However, for higher value of n, it starts giving the very weird value. Some people told me that it is because of the method of integration and suggested that I should use other numerical method other than just use the default Integral of Mathematica. I would really appreciate it you you guys can help me with this problem.

$$\text{\chi n}(\text{n\_},\text{z\_})\text{:=}-\frac{(2 z) e^{-\frac{z}{a_0 n}} L_{n-1}^1\left(\frac{2 z}{a_0 n}\right)}{\sqrt{a_0 n^3} \left(a_0 n\right)}$$

Here I attach the image from my Mathematica notebook

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Young Jul 17 '16 at 23:19
• People here generally like users to post code as Mathematica code instead of images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful – Michael E2 Jul 18 '16 at 0:13
• @Young, an integral equation does not mean what you think it means. – J. M.'s technical difficulties Jul 18 '16 at 2:39
• Ah, thanks @J.M. ... I admit I didn't read the tag description ... I will try to avoid that mistake in the future. – Young Jul 18 '16 at 2:42

The issue is the lack of precision specified for a0 (Reference: precw). Change to the following:

a0 = 36/10*10^-10;


and the results seem appropriate.

xn[n_, z_] := -1/Sqrt[n^3 a0] (2 z)/(n a0) Exp[-z/(n a0)] LaguerreL[n - 1, 1, (2 z)/(n a0)]

data = N[Table[{n, Integrate[xn[n, z] z xn[n, z], {z, 0, Infinity}]/a0}, {n, 1, 30}]]


(Note: For some versions of Mathematica, including V10 Student Edition, the Integrate function needs to include PrincipalValue -> True for proper evaluation.

{{1., 1.5}, {2., 6.}, {3., 13.5}, {4., 24.}, {5., 37.5}, {6., 54.}, {7., 73.5}, {8., 96.}, {9., 121.5}, {10., 150.}, {11., 181.5}, {12., 216.}, {13., 253.5}, {14., 294.}, {15., 337.5}, {16., 384.}, {17., 433.5}, {18., 486.}, {19., 541.5}, {20., 600.}, {21., 661.5}, {22., 726.}, {23., 793.5}, {24., 864.}, {25., 937.5}, {26., 1014.}, {27., 1093.5}, {28., 1176.}, {29., 1261.5}, {30., 1350.}}

As an interesting aside this reduces to a fairly simple relationship:

$\chi n = (3/2) n^2$

nlm = NonlinearModelFit[data, a x ^2, {a}, x]
Plot[nlm[x], {x, 0, 30}, Epilog -> Point[data]]


Another way to identify this relationship was suggested by Bob Hanlon:

FindSequenceFunction[data[[All, 2]] // Rationalize, n]

• Thank you very much for your help. I really appreciate it. However, when I try the same thing, some of the result turn to be zero. Can you please explain a little bit more about the working precision. I look it up but it confuses me alot – Quang Phan Jul 17 '16 at 23:44
• When I try your code, it happens like this: 1. 1.5 2. 6. 3. 13.5 4. 24. 5. 37.5 6. 54. 7. 73.5 8. 0. 9. 0. 10. 0. 11. 0. 12. 0. 13. 0. 14. 0. 15. 0. 16. 0. 17. 0. 18. 0. 19. 739.850168393288 20. 0. 21. 891.9849790911463 22. 900.1037074992012 23. 1049.7176922474973 24. 664.0116898665082 25. 1116.1338233731344 26. 952.0162224510896 27. 1042.9886917686958 28. 1114.5273086839672 29. 1400.6626376639572 30. 1372.4169189319132 – Quang Phan Jul 17 '16 at 23:54
• I use Version 10-Student Edition – Quang Phan Jul 17 '16 at 23:57
• Hi, I did try it but it keeps giving me a bunch of zeros as mentioned above. – Quang Phan Jul 18 '16 at 0:22
• It says something like this "Integral of E^(-6250000000 z/9) z^3 (301327047-732392128125000000 z+<<6>>+46938657760620117187500000000000000000000000000000000000000 z^6-582076609134674072265625000000000000000000000000000000000000000000 z^7)^2 does not converge on {0,[Infinity]}. >>" – Quang Phan Jul 18 '16 at 0:24