Revisiting the Thomas–Fermi equation: Accelerating rational Chebyshev series through coordinate transformations

We revisit the spectral solution of the Thomas–Fermi problem for neutral atoms, u r r − ( 1 / r ) u 3 / 2 = 0 on r ∈ [ 0 , ∞ ] with u ( 0 ) = 1 and u ( ∞ ) = 0 to illustrate some themes in solving differential equations when there are complications, and also to make improvements in our earlier treatment. By “complications” we mean features of the problem that either destroy the exponential accuracy of a standard Chebyshev series, or render the classic Chebyshev approach inapplicable. The Thomas–Fermi problem has four complications: (i) a semi-infinite domain r ∈ [ 0 , ∞ ] (ii) a square root singularity in u ( r ) at the origin (iii) a fractional power nonlinearity and (iv) asymptotic decay as r → ∞ that includes negative powers of r with fractional exponents. Our earlier treatment determined the slope at the origin to twenty-five decimal places, but no fewer than 600 basis functions were required to approximate a univariate solution that is everywhere monotonic, and all of the earlier tricks failed to recover an exponential rate of convergence in the truncation of the spectral series N, but only a high order convergence in negative powers of N. Here, using the coordinate z ≡ r to neutralize the square root singularity as before, we show that accuracy and rate of convergence are significantly improved by solving for the original unknown u ( r ) instead of the modified unknown v ( r ) = u ( r ) used previously. Without a further change of coordinate, a rational Chebyshev basis T L n ( z ; L ) yields twelve decimal place accuracy for the slope at the origin, u r ( 0 ) , with 70 basis functions and twenty-four decimal places with a truncation N = 100 . True exponential accuracy can be restored by using an appropriate change of coordinate, z = G ( Z ) , where G is some species of exponential. However, the various Chebyshev and Fourier series for the Thomas–Fermi function have “plural asymptotics”, that is, a n ∼ a intermediate ( n ) for 1 ≪ n ≪ n I but a n ∼ a far ( n ) for n ≫ n I for some positive constant n I . The “far-asymptotics” as n → ∞ are often of no practical significance. For this problem, the TL-with-sinh method, theoretically the best for huge n, is bedeviled with numerical ill-conditioning and its asymptotic superiority is realized only when the goal is at least forty decimal places of accuracy. This is absurd for engineering, but useful perhaps for benchmarking. To sixty decimal places, u r ( 0 ) = − 1.588071022611375312718684509423950109452746621674825616765677 . To nineteen digits after the decimal point, the constant in the Coulson–March asymptotic series is improved to F = 13.2709738480269351535 .