Laguerre basis and pseudospectral method

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the pseudospectral method in conjunction with the Laguerre polynomials.

ClearAll["Global`*"];
Off[General::spell, General::spell1];

Introduction

Review of the theory

General remarks

In notebook [4] a general overview of collocation and pseudospectral methods has been presented.  Here we specialize to the problem of the solution of radial Schrödinger equation for the reduced radial function, ψ = r R[r]. We repeat some of the ideas introduced in notebook [4] to make the notebook self-contained, and we give a detailed discussion on the boundary conditions.

In general a collocation method is defined by the choice of a grid points (Subscript[x, a]) and by the discretization of the equation on this grid:

-(1/2)(Subscript[𝒟, ab]^((2)))(Subscript[ψ, b]) + V[(Subscript[x, a])](Subscript[ψ, a]) = E (Subscript[ψ, a])

The derivative matrix 𝒟 has in any finite difference scheme elements different from zero for |a-b| ≤ k, for example in the nearest neighbors approximation of second derivative on a uniform grid

f''[(Subscript[x, a])] ≃ (1/(2(h^2)))(f[(Subscript[x, a])+h] - 2 f[(Subscript[x, a])] + f[(Subscript[x, a])-h])

we have a tridiagonal matrix, i.e. k=1. The number of points involved in the computations of derivatives can be related to the degree of a polynomial which interpolate the function (in general a local interpolation): the approximation are made as to be exact for these polynomials. For instance polynomials of second degree generate formula (). A collocation method is defined with a maximum degree polynomial, i.e. for a grid of (Subscript[n, G]) points by a polynomial of degree (Subscript[n, G])-1. In this case the matrix 𝒟 is in general a full matrix and all points concur to define the derivative at each point of the grid.

The best grid from the numerical point of view are related to the grid points used by Gaussian integration formulas, i.e. are zeros of an n-degree polynomial, orthogonal with respect to some measure on the given interval. In this case collocation method is also called pseudospectral method. The reason being that equations () are seen to be equivalent to the usual spectral equations, as those used in notebooks [6], if matrix elements are computed by Gaussian integration formulas, as we have indeed done in notebook [7].

In contrast to the spectral methods equation () is defined directly on "real" space, not in the space of coefficients of a spectral expansion. Moreover it is often much faster to solve: for an n-grid points we have to compute only n values of the potential, while we have to compute of the order of n^2/2 values of the corresponding matrix elements in the spectral method.

The grid : internal points and boundary

Let us take a grid consisting of  n internal points in the interval [0 , ∞]. The values of the wave function on these points will be our unknown quantities. We have to embed the grid and the values in some interpolation scheme which allows the use of derivatives, acting on the interpolated function. One must choose between two possibilities:

The second procedure is more flexible, as it can account for different boundary conditions, and is also more easily programmable, so we follow this second alternative. Our extended grid will contain

(Subscript[n, E])= n + 1

grid points (our boundary is only r = 0 point).

The set of cardinal functions

The particular grid chosen depends on which set of functions are used for interpolation. As discussed briefly in notebook [], convergence properties strongly suggest a choice based on set of orthogonal polynomials. These orthogonal polynomials will be in general (but not necessarily) multiplied by a "weight factor" which we call α(x). typically (α^2) will be the weight factor with respect to polynomials are orthogonal. For Laguerre polynomials α[x] = Exp[-x/2]. The set of interpolating functions will be the linear space (Subscript[𝒱, n]) spanned by functions of the type

α(x) (Subscript[p, k])[x];    deg[(Subscript[p, k])] ≤ n

Let us stress a few points :

(Subscript[f, 0])[x] = Exp[-x/2];  (Subscript[f, k])[x] = Exp[-x/2]x L[k-1,x];  k = 1,2,...n

The f' s with k ≥ 1 are exactly those used in the spectral decomposition of the problem.

This choice is somewhat arbitrary, had we chosen α (x) = x Exp[-x/2] the vanishing boundary conditions had been included automatically, but the expression of the derivative matrix given below would be different.

Once fixed our space of interpolating functions, (Subscript[𝒱, n ])it is trivial to define derivatives. The interpolating function for ψ (which will be called always ψ) is

ψ[x] = (∑_(a=0)^n)((α[x])/(Subscript[α, a]))(Subscript[ℒ, a])[x](Subscript[ψ, a])= (∑_(a=0)^n)(Subscript[C, a])[x](Subscript[ψ, a]);     (Subscript[α, a])≡ α[(Subscript[x, a])];

ℒ are Lagrange polynomials of degree n, (see notebook [3]) with explicit expression

(Subscript[ℒ, i])[x]=(∏_(j≠i)^n)(((x-(Subscript[x, 0]))(x-(Subscript[x, 1]))..(x-(Subscript[x, i-1]))(x-(Subscript[x, i+1])).. (x-(Subscript[x, n])))/(((Subscript[x, i])-(Subscript[x, 0]))((Subscript[x, i])-(Subscript[x, 1]))..((Subscript[x, i])-(Subscript[x, i-1]))((Subscript[x, i])-(Subscript[x, i+1])).. ((Subscript[x, i])-(Subscript[x, n]))))

which have the trivial but crucial property :

(Subscript[ℒ, i])[(Subscript[x, a])]= (Subscript[δ, ia])

The functions C, i.e. polynomials times exponentials in our case, will be called cardinal functions, and they also have the property

(Subscript[C, i])[(Subscript[x, a])]= (Subscript[δ, ia])

Derivative are now immediate :

Dψ [x] = (∑_(b=0)^n)(Subscript[C, b]^('))[x](Subscript[ψ, b]);   ⇒ Dψ [(Subscript[x, a])] = (∑_(b=0)^n) (Subscript[DC, b])[(Subscript[x, a])](Subscript[ψ, b])≡ (Subscript[𝒟, ab])(Subscript[ψ, b])

Similarly for any derivative :

(Subscript[𝒟, ab]^((L)))= ( (D^((L)))(Subscript[C, b]))[(Subscript[x, a])]

Since  cardinal functions are explicitly known  it is just a technical problem to compute matrix of derivative, and there will be a function for that.

We want to underline some aspects of what had been discussed up to now:

For any function in (Subscript[𝒱, n])we can write an expansion in term of basis functions, we can take the basis () as they are more convenient for later development:

∀ ψ∈(Subscript[𝒱, n]) :   ψ[x] = (∑_(k=0)^n)(Subscript[c, k])(Subscript[f, k])[x]

The sum starts from index 0 and contains n + 1 terms. From () evaluated on the grid points it follows that there exist a linear relation between function's values and expansion coefficients:

ψ[(Subscript[x, a])]≡(Subscript[ψ, a]) = (∑_(k=0)^n) (Subscript[c, k])(Subscript[f, k])[(Subscript[x, a])] ≡ (∑_(k=0)^n)(Subscript[S, ak])(Subscript[c, k]);             (Subscript[c, k])=(∑_(a=0)^n) (Subscript[S, ka]^(-1)) (Subscript[ψ, a]);          (Subscript[S, ak]) = (Subscript[f, k])[(Subscript[x, a])].

To impose the validity of Schrödinger equation and boundary conditions on ψ values means:

-(1/2)(∑_(b=1)^n)(Subscript[𝒟, ab]^((2)))(Subscript[ψ, b]) + V[(Subscript[x, a])](Subscript[ψ, a]) = E (Subscript[ψ, a]);  a= 1,...n.     and      (Subscript[ψ, 0])=0.

The collocation points for Schrödinger equation are internal points of the grid. The sum on b index in () starts from 1 as (Subscript[ψ, 0])= 0 does not contribute. This means that the matrices which appear in this equation are obtained by the larger (n+1)×(n+1) matrices () by dropping first the row and the first column. This mechanism has already been presented in notebook [4].

The relation to spectral method and the choice of the grid

Up to this point we have carefully avoided any statement on the nature of the grid. The reader has also noticed that the choice of the orthogonal polynomials practically had no role up to now.

The choice of the grid is dictated both from numerical experience and theoretical analysis on the stability properties and on coincidence of collocation method with usual spectral method in some limiting case. For the spectral method there exist in literature several criteria for stability, bounds on errors as a function of dimension of the basis, etc. All this is based on the variational formulation of the problem. It is important then to link collocation method with this more "stable" method. The grid must to be chosen as  set of zeros of orthogonal polynomials relative to the interval we are studying. The statement is a little vague as there can be several orthogonal sets, depending on the measure used on the interval. So let us express the idea in these terms: to each spectral method, based on an orthogonal set of polynomials, there corresponds a collocation method based  on the grid of zeros of these polynomial.

We have still an alternative: we can choose a Gauss basis, i.e. based on the internal points, plus the boundary point or a Gauss-Lobatto basis, i.e. with boundary points included in the grid from the beginning. The corresponding roots and weights are slightly different but almost equivalent from the numerical point of view. In this notebook we will leave the choice to the reader, pointing out the small differences in the two cases.

For radial problems a natural basis is the one based on the Laguerre polynomials, so there are two choices:

In both cases we have n internal points and 1 boundary point in the extended grid.

Let us now discuss the relation to the spectral method. Due to the relation () the equation () can be transformed into a set of equations for the coefficients (Subscript[c, k])of the approximant ():

-(1/2)(∑_(b=1)^n)(∑_j)(Subscript[𝒟, ab]^((2)))(Subscript[S, bj])(Subscript[c, j]) + (Subscript[V, a])(∑_j)(Subscript[S, aj])(Subscript[c, j]) = E (∑_j)(Subscript[S, aj])(Subscript[c, j]) ; a>0;         (Subscript[ψ, 0])= 0 =(∑_k) (Subscript[S, 0k])(Subscript[c, k])

Let us take as a basis the set (). This is the one selected, for k ≥ 1,  by the spectral methods, as the basis functions in this case must obey boundary conditions. For this basis clearly (Subscript[ψ, 0])= (Subscript[c, 0]) as all basis functions with index greater than zero vanish at the origin. This means that in this basis

(Subscript[S, 00])=1;   (Subscript[S, 0i])= 0; i ≥1

As index a in () is greater than 0, equation () implies that all indexes in () run from 1 to n, the zero being excluded.

We show now that, if matrix elements are computed by Gaussian formulas, equation () is equivalent to the equation for the spectral coefficients.  Let us call these coefficients (Subscript[u, k]), they satisfy

-(1/2)(∑_(j=1)^n)(Subscript[D, ij])(Subscript[u, j]) + (∑_(j=1)^n) (Subscript[V, ij])(Subscript[u, j]) = E (∑_(j=1)^n) (Subscript[N, ij])(Subscript[u, j])

The reader has noticed that the index j runs from 1 to n. Let us now compute the matrix elements with Gaussian formulas for the internal points of the grid: this is correct as the basis functions (Subscript[f, k])[x],k ≥ 1, vanish at the boundary. For derivatives at grid points the relation () can be used which is valid, as we have stressed before, for these functions:

(Subscript[D, ij]) = (∫_0^∞)𝕕x (Subscript[f, i])[x](∂_xx)(Subscript[f, j])[x] = (∑_a) (Subscript[W, a]) (Subscript[f, i])[(Subscript[x, a])](∑_b) (Subscript[𝒟, ab]^((2)))(Subscript[f, j])[(Subscript[x, b])] =  (∑_(a,b)) (Subscript[W, a]) (Subscript[S, ai]) (Subscript[𝒟, ab]^((2)))(Subscript[S, bj])

(Subscript[V, ij])= (∫_0^∞)𝕕x (Subscript[f, i])[x]V[x](Subscript[f, j])[x]= (∑_a) (Subscript[W, a]) (Subscript[f, i])[(Subscript[x, a])](Subscript[V, a])(Subscript[f, j])[(Subscript[x, a])]= (∑_a) (Subscript[W, a]) (Subscript[S, ai]) (Subscript[V, a])(Subscript[S, aj])

(Subscript[N, ij])= (∫_0^∞)𝕕x (Subscript[f, i])[x](Subscript[f, j])[x]= (∑_a) (Subscript[W, a]) (Subscript[f, i])[(Subscript[x, a])](Subscript[f, j])[(Subscript[x, a])]= (∑_a) (Subscript[W, a]) (Subscript[S, ai]) (Subscript[S, aj])

Here for shortness we wrote  (Subscript[W, a])= (Subscript[w, a])Exp[(Subscript[x, a])];   (Subscript[w, a]) are  the usual weights.

With (Subscript[M, ia])= (Subscript[W, a]) (Subscript[S, ai]) equation () can be written in the form

(∑_a)(Subscript[M, ia]) ( -(1/2)(Subscript[𝒟, ab]^((2)))(Subscript[S, bj])(Subscript[u, j])+ (Subscript[V, a])(Subscript[S, aj]) (Subscript[u, j]) - E (Subscript[S, aj])(Subscript[u, j])) = 0

As the matrix M is in general non-singular this equation coincides with eq.().

N.B. A much shorter "proof" could be obtained assuming for (∂_xx)ψ an expansion through matrix (Subscript[D, ij]), i.e. writing () in x-space but we prefer the proof above as it make clear where and when boundary conditions come in.

In conclusion the two approach are equivalent when Gaussian integrals are exact.

In this respect Gauss grid and Gauss-Lobatto grids are slightly different. The basis functions are polynomials of at most n-degree times an exponential. For matrix (Subscript[N, ij])this means that all elements are weighted integrals of polynomials of degree less than 2n. For all these polynomials Gauss lobatto formula is exact, even if we limit sums at internal points, as functions vanish at the boundary. Gauss formulas are correct up to polynomials of degree n-1, this means that element (Subscript[N, nn])is reproduced only approximately. For the Laguerre basis the same phenomenon occurs for the matrix (Subscript[D, ij]). In the derivatives of Exp[-x/2] x L[k-1,x] we have the derivative of the polynomial, which lower his degree, and the derivative of the exponential, which leave it unchanged. In matrix element (Subscript[D, ij]) we have weighted polynomials of maximum degree i+j.  This means that elements (Subscript[D, nn])is exactly reproduced in the Gauss Lobatto grid while only approximately reproduced in the Gauss case. As n goes large evidently the two procedures give the same results. The only difference we have been able to detect will be presented below, and concerns the levels of hydrogen atom. The potential 1/r happens to be exactly reproduced in both Gauss and Gauss Lobatto basis, as the negative power cancels the x factors in the basis functions and we end always with polynomials. In this case the collocation method must be identical to the spectral method. It will be seen that this is true for physical relevant region for both methods, but from a mathematical point of view the Gauss method has an additional instability.

The previous relations allow a number of checks. First, as already stated, relation () holds in (Subscript[𝒱, n])space and in particular for basis functions. We must have

(Subscript[S, ai]^((2)))≡ (D^((2)))(Subscript[f, i])[(Subscript[x, a])] = (Subscript[𝒟, ab]^((2)))(Subscript[f, i])[(Subscript[x, b])] =(Subscript[𝒟, ab]^((2))) (Subscript[S, bi])⇒ (Subscript[𝒟, ab]^((2)))= (Subscript[((S^((2)))(S^(-1))), ab])

As the right hand side can be easily computed, this provides a check for matrix 𝒟. This relation is valid for any grid, and has nothing to do with the Gauss grid or similar things. It depends only of our choice of the basis.

A second partial check come from the N matrix. As was already stated one must have, with at most one exception,

(Subscript[N, ij]) =  (∑_a) (Subscript[W, a]) (Subscript[S, ai]) (Subscript[S, aj]) =  (∑_a) (Subscript[W, a]) (Subscript[f, i])[(Subscript[x, a])](Subscript[f, j])[(Subscript[x, a])].

This is a check for the grid points and weights.

Interpolation and norms

A final remark about the use of future results. Suppose that one has solved problem ().  One will obtain eigenvalues and eigenfunction at grid points. As usual we will use a scale parameter λ, i.e. in terms of real radial variable r our basis functions will be of the form

(Subscript[F, k])[r] = (Subscript[f, k])[(r/λ)]=Exp[-r/λ] (r/λ) L[k-1,r/λ]

The parameter will be reabsorbed with the change of variables r = λ x. This means that the grid will be in x variable and in particular (Subscript[V, a])will be computed as V[λ (Subscript[x, a])].  In r - space the grid will be (Subscript[r, a])= λ (Subscript[x, a]). An optimal λ value will be fixed by a  search for stationary ranges for λ.

For functions defined on the same grid scalar product, will be given by

∫φ[r]ψ[r] 𝕕r = λ ∫φ[x]ψ[x] 𝕕x = λ ∑ (Subscript[W, a])(Subscript[φ, a])(Subscript[ψ, a])

In particular (we suppose all quantities real):

((|ψ|)^2) = λ ∑ (Subscript[W, a])(Subscript[ψ, a])(Subscript[ψ, a]).

For derivative operators

∫φ[r](D^((L))) ψ[r] 𝕕r =λ (1/(λ^L))∫φ[x](D^((L)))ψ[x] 𝕕x = (λ^(1-L)) ∑ (Subscript[W, a])(Subscript[φ, a])(Subscript[𝒟, ab]^((L)))(Subscript[ψ, b]).

This in principle could allow to compute any physical quantity, but for comparison with analytically known results or just for graphical reasons, sometimes the values of ψ at different points are needed. They can be obtained in a variety of direct ways, all quite consuming in computer time, such as computing spectral coefficient by the matrix (Subscript[f, i])[(Subscript[x, a])]and using analytical expression for the basis. A more logical way is to use an interpolation. Unfortunately standard Mathematica function suffers from some difficulties in  this case as the grid point are quite far apart. A better way is to use a so-called "baricentric interpolation", which will be briefly explained in the section dedicated at the functions used in the notebook. This procedure works very well for relatively small basis, up to n=40 typically.  It has some instability, which we have cured in a rather crude way, i.e. by using "bad points" Mathematica procedure.

References for the code used

To implement the method we need only two things : a function which computes grid points and a function which computes the  derivative matrices. The first one has been used many times in other notebooks, it is laguerreRoots, which computes both grid points and weights for Gaussian integration. The second one is more problematic. In principle one can write explicitly polynomials (Subscript[ℒ, i]) , make derivatives and evaluate them at the grid points (we will do this below as a check) but this has two drawbacks:

It is slow, very slow.

It is extremely unstable. The Lagrange polynomial of high degree are known to be very unstable, some examples have been presented in notebook [7], below we will show some other similar phenomena.

The way out is to write carefully a routine which minimizes the rounding errors. Mathematica has a similar routine for the simpler case of Chebyshev polynomials, as we will see, but not for grids based on the Laguerre polynomials. This problem has obviously been studied in the literature and we will adapt for Mathematica a routine written for MatLab by J.A.C. Weideman and S.C. Reddy, see ref.[1]. The main difference with the quoted work is the optional use of the Gauss Lobatto grid.

The same reference [1] has been partially followed for the interpolation problem.

Functions

Functions used for spectral method

laguerreMatrices

Compute by default matrices D2 and matN. An optional list containing powers of r grater than -2 can be given as a second argument. For instance {-1,2} will give matrices for 1/r and (r^2). These matrices are computed in exact arithmetic. In this notebook we will use the simple call sequence

{D2, matN} = laguerreMatrices[dim]

for dimension dim bases.

laguerreMatrices[dim_,listV_:{}]:= Module[{normarray,Nf,d2,D2,rm1,Rm1,δdim,listOfPotentials},
δdim = If[Length[listV]== 0,0,Max[Abs[listV]]+2];
Block[{n,j},
n = j-1;
normarray=SparseArray[{{j_,j_}->(1/2)(6 n^2+6n+2),{k_,j_}/;j-k==1->-4 n^2,
{k_,j_}/;j-k==2-> n(n-1)},{dim,dim}];
Nf =(normarray+Transpose[normarray]);
d2=Normal[SparseArray[{{k_,j_}/;j-k==0->-((n^2)+n+1)/4,
{k_,j_}/;j-k==2->n(n-1)/4},{dim,dim}]];
D2 =(d2+Transpose[d2]);
If[Length[listV]== 0,Return[{D2,Nf}]];
rm1 =  SparseArray[{{k_,j_}/;j-k==0->(2n+1)/2,
{k_,j_}/;j-k==1->-n},{dim+δdim,dim+δdim}];
Rm1 = (rm1+Transpose[rm1]);
listOfPotentials = {};
Do[mat0=Rm1;
mat1=If[listV[[sp]]== -1,mat0,Nest[Dot[#1,mat0]&,mat0,listV[[sp]]+1]];
mat1= mat1[[1;;dim,1;;dim]];
listOfPotentials = Append[listOfPotentials,mat1];
,{sp,1,Length[listV]}];
listOfPotentials = Normal/@listOfPotentials;
];
Normal/@Join[{D2,Nf},listOfPotentials]]

laguerreRoots

Compute by default n-points abscissas and weights for Gauss-Lobatto grid.  This is different form the similar routine used in other notebooks. An optional argument, ord, allow to compute the same quantities for Gauss grid. ord=0 gives 0 order Laguerre polynomials, i.e. Gauss grid, ord = 1 Gauss Lobatto grid. In both cases only internal points are given. The standard calling sequence is

{xgrid, weights} = laguerreRoots[n,ord]

A third optional argument, mostly used for checks, allow to do computation in arbitrary precision. If you use this argument you must specify the order.

laguerreRoots[n_,ord_:1,prec_:MachinePrecision]:=
Module[{J,δ,γ,radEig,w},
If [ord==1,
δ = 2Range[1,n]; γ = -N[(Sqrt[Range[n-1]Range[2,n]]),prec];,
δ = 2Range[0,n-1]+1; γ = -Range[1,n-1];];J = Normal[SparseArray[{Band[{1,1}]-> δ,Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
radEig=Sort[Eigenvalues[N[J,prec]]];
If[ord==1 ,w= (1/(n+1))(1/((LaguerreL[n+1,1,radEig])^2));,
w=(radEig/((((n+1))^2)((LaguerreL[n+1,radEig])^2)));];
{radEig,w}]

getMatV

This function computes matrix elements for the potential V. As inputs it requires: the scale λ, the symbol for the potential function, example V if the function has been defined as V[r], a list with the output of laguerreRoots, i.e. {xgrid, weights}, a matrix, matBaseLag which contains values of Laguerre polynomials on the points xgrid. The standard calling sequence will be

matV = getMatV[λ, V, {xgrid, weights}, matBaseLag]

getMatV[λ_?NumberQ,V_,{xgrid_,weights_},matBaseLag_]:=
Module[{dim,matV1,matV2,Vk},
dim = Length[matBaseLag];
matV1 = 0*IdentityMatrix[dim];
Vk = V[xgrid*λ];
Table[matV1[[i,j]]=(matBaseLag[[i]]*weights*Vk).matBaseLag[[j]];,
{i,1,dim-1},{j,i+1,dim}];
matV2 =Table[ (matBaseLag[[i]]*weights*Vk).matBaseLag[[i]],{i,1,dim}];
DiagonalMatrix[matV2] + matV1 + Transpose[matV1]];

Energy

This function simply collects computed quantities as xgrid, etc, computes matV with a call to getMatV and computes the  eigenvalues with a call to Mathematica function Eigenvalues. It is written for computing  also the higher momenta states, in our basis the operator 1/r^2 is represented by the identity matrix. The calling sequence is

energy[λ, {state, L}, V, {xgrid, weights}, matBaseLag]

energy[λ_?NumberQ,{state_,L_},V_,{xgrid_,weights_},matBaseLag_]:= Module[{matV,dim},
dim = Length[matBaseLag];
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
Sort[Re[Eigenvalues[{((- 1)/(2(λ^2))) D2+matV + (1/(2(λ^2)))L(L+1)IdentityMatrix[dim],matN}]]][[state]]]

Functions used for pseudospectral method

laguerredif

This function writes the derivative of Exp[-x/2] for the computation of (Subscript[𝒟, ab])matrices and makes a call to polidif, which does the derivative of the polynomial part, i.e. the difficult work. The function can compute in fact the derivatives of every order and the output could contain all, but to simplify the work for the reader we suppress all other matrices except the requested one, usually the second derivative. The mandatory arguments are n, the dimension of the grid, and m the order of the derivative (usually 2). An optional argument prec can be used to perform arbitrary precision calculations. The typical calling sequence is

{xgrid, D2PS} = laguerredif[dim, 2];

The output contains the grid of points (without the zero) and the matrix for second derivatives. They are a vector of dim components and a dim×dim matrix.

The optional argument ord specifies if the grid is a Gauss Lobatto grid (ord =1 , the default) or ord = 0 (a Gauss grid).

The argument prec can be used to switch to another output. If prec < 0 the extended grid (with 0 point included0 and extended matrices form first to m-th derivatives are produced. This can be of some help in checks and in particular applications. The reader must pay attention that in this case the output has the form

{xgrid, {D1PS, D2PS}} = laguerredif[dim, 2,ord, -1];

and the task of eventually reduce matrices and grid to internal points is left to the user. If this argument is used ord must be specified.

laguerredif[n_,m_,ord_:1,precIn_:MachinePrecision]:=
Block[{prec},
Module[{x,ell,beta,alpha,fullOutput},
fullOutput = False;
If[precIn<0,prec = MachinePrecision;fullOutput = True;,prec = precIn];
 x = laguerreRoots[n,ord,prec][[1]];
x = Prepend[x,0];
alpha = N[Exp[-x/2],prec];
beta = ConstantArray[0,{m,Length[x]}];
Do[beta[[ell]] = N[(((-1)/2))^ell*ConstantArray[1,n+1],prec];
,{ell,1,m}];
DM = polidif[x,alpha,beta,prec];
If[fullOutput== True,Return[{x,DM}]];
x = Take[x,{2,n+1}];
Do[
DM[[ell]] = Take[DM[[ell]],{2,n+1},{2,n+1}];,{ell,1,m}];
Return[{x,DM[[m]]}];]]

polidif

This function computes the derivative of the interpolating polinomial. For details see ref.[1]. The reader does not need to call this function, it is used only by laguerredif.

polidif[x_,alpha_,bIn_,prec_:MachinePrecision]:=Module[
{n,m,XXT,c,DX,id,k,mC,mCT,mZ,mX,mD,mY,bb,ell},
(n = Length[x];
b = N[bIn,prec];
m = Dimensions[b][[1]];
id = IdentityMatrix[n];
XXT = PadRight[{x},n,{x}];
DX = Transpose[XXT] - XXT;
Table[DX[[k,k]]=1,{k,1,n}];
c = alpha*Apply[Times,DX,1];
mCT = PadRight[{c},n,{c}];
mC = Transpose[mCT]/mCT;
mZ = N[1/DX,prec];
Table[mZ[[k,k]]=0,{k,1,n}];
bb = Drop[Flatten[mZ],{1,n^2,n+1}];
mX = Transpose[Partition[bb,n-1]];
mY = ConstantArray[N[1,prec],{n-1,n}];
(* mY = SparseArray[{},{n-1,n},1]; *)
mD = IdentityMatrix[n];
Array[pD,m];
Do[(
mY = Join[{b[[ell,All]]},ell*Take[mY,{1,n-1}]*mX,1];
Table[mY[[i]]=mY[[i-1]]+mY[[i]],{i,2,n}];
mD = ell*mZ*(mC*Transpose[PadRight[{Tr[mD,List]},n,{Tr[mD,List]}]]-mD);
Table[mD[[i,i]]=mY[[n,i]],{i,1,n}];
pD[ell] = mD;),{ell,1,m}];
Return[Table[pD[k],{k,1,m}]];)]

Functions used for interpolation

polint

This function computes a polynomial interpolation for a function, using the exponential weights Exp[-x/2]. For details on this function, which will be checked below, see ref.[1]. The critical points, i.e. the instabilities, come from the last instruction, where the denominator mD.w can be zero in the usual Mathematica precision. The needed input is:

In output one gets the values of the function on the new grid. The calling sequence is

ψnew = polint[xgrid, ψgrid, xnewgrid];

polint[xk_,fkin_,x_]:=
   Module[{n, m, mDT, p, k, fk, mD1,alpxk,alpx},
 n = Length[xk];
 m = Length[x];
alpxk = Exp[-xk/2];
alpx = Exp[-x/2];
fk = fkin/alpxk;
mD = Outer[Plus,xk,-xk];
Table[mD[[k, k]] = 1, {k, 1, n}];
w = Apply[Times, 1/mD, 1];
 mD =Outer[Plus,x,-xk];
 mD = 1/(mD + $MachineEpsilon*Map[Boole[#1== 0]&,mD,{2}]); 
p = alpx*(mD.(w*fk)/(mD.w))];

interpLaguerre

As stated above we have wrapped polint[] with a selection mechanism for bad points, where instabilities produce a Complex Infinity result in Mathematica. These point are selected and the function is interpolated at these points with usual Interpolation[] function of Mathematica . This procedure is taken for all points for dim > 45.  It needs as input:

The calling sequence is

ψR = interpLaguerre[{xgrid, eigV}, λ, r]

interpLaguerre[{xgrid_,ψgrid_},λ_,r_]:= 
Module[{temp,dim,weights,rgrid,fPoli,mask,rmask,fmask,fpo,normPsi},
dim = Length[xgrid];
rgrid=xgrid*λ;
If[dim>45,
 psi=Interpolation[Transpose[{rgrid,ψgrid}],InterpolationOrder-> 6];
Return[psi[r]]];
fPoli=polint[xgrid,ψgrid,r/λ];
mask= Boole[Map[NumberQ[#1]&,fPoli]];
rmask = Select[r*mask,#1≠0&];
If[Length[rmask]< Length[r],
fmask = Select[fPoli,NumberQ[#1]&];
fpo=Interpolation[Transpose[{rmask,fmask}],InterpolationOrder-> 1];
fPoli = fpo[r];];
fPoli];

Checks

We consider a low-dimensional case to check explicitly the relations between the spectral and pseudospectral methods, expression of interpolating polynomials etc.

Matrix 𝒟

Definition with cardinal functions

The cardinal functions are given by

(Subscript[C, i])[x] = ((α[x])/(Subscript[α, i]))(Subscript[ℒ, i])[x]= ((Exp[-x/2])/(Subscript[α, i]))(Subscript[ℒ, i])[x]

then for their derivatives :

(Subscript[C, i]^('))[x] = (-(1/2)(Subscript[ℒ, i])[x] + (Subscript[ℒ, i]^('))[x])((α[x])/(Subscript[α, i]));       (Subscript[C, i]^(''))[x] = ((1/4)(Subscript[ℒ, i])[x] - (Subscript[ℒ, i]^('))[x]+(Subscript[ℒ, i]^(''))[x])((α[x])/(Subscript[α, i]));

On grid points :

(Subscript[𝒟, ab]^((1)))= (Subscript[C, b]^('))[(Subscript[x, a])]=(-(1/2)(Subscript[ℒ, b])[(Subscript[x, a])] + (Subscript[ℒ, b]^('))[(Subscript[x, a])])((Subscript[α, a])/(Subscript[α, b]));   (Subscript[𝒟, ab]^((2)))= ((1/4)(Subscript[ℒ, b])[(Subscript[x, a])] - (Subscript[ℒ, b]^('))[(Subscript[x, a])]+(Subscript[ℒ, b]^(''))[(Subscript[x, a])])((Subscript[α, a])/(Subscript[α, b]));

The Gaussian weights are given by

(Subscript[w, i])= (∫_0^∞)𝕕x Exp[-x] (Subscript[ℒ, i])[x]

The following function computes directly these quantities.

wdFromLagrange[xgridInt_]:= Module[{ℒ,dim,w,der,expa,xgrid},
xgrid = Join[{0},xgridInt];
dim = Length[xgrid];
ℒ[x_]:=Table[InterpolatingPolynomial[Transpose[{xgrid,SparseArray[{k}-> 1,{dim}]}],x],{k,1,dim}];
SetAttributes[ℒ,Listable];
w=Table[NIntegrate[ℒ[x][[k]] Exp[-x],{x,0,∞}],{k,1,dim}];
(*expa[x_]:=Exp[-x/2]*Exp[xgrid/2]; *)
(* SetAttributes[expa,Listable]; *)
expa = Outer[Divide,Exp[-xgrid/2],Exp[-xgrid/2]];
der=Transpose[((1/4)ℒ[xgrid]-ℒ'[xgrid]+ℒ''[xgrid])]*expa;
{w,der}]

Check for Gauss Lobatto grid

ordLag = 1;

dim=3;
{xgrid,weights} = laguerreRoots[dim,ordLag]

pseudospectral derivative

{rgrid,D2PS} = laguerredif[dim,2,ordLag];
rgrid== xgrid

Here we compute weights and derivatives directly

{wL,dL} = Chop[wdFromLagrange[xgrid]];

Max[Abs[wL[[-dim;;]]- weights]]

MatrixForm/@{D2PS,dL}

Max[Abs[dL[[-dim;;,-dim;;]]-D2PS]]

Check for Gauss  grid

ordLag = 0;

dim=3;
{xgrid,weights} = laguerreRoots[dim,ordLag]

pseudospectral derivative

{rgrid,D2PS} = laguerredif[dim,2,ordLag];
rgrid== xgrid

Here we compute weights and derivatives directly

Off[NIntegrate::"slwcon",NIntegrate::"ncvb"]
{wL,dL} = Chop[wdFromLagrange[xgrid]];
On[NIntegrate::"slwcon",NIntegrate::"ncvb"]

Max[Abs[wL[[-dim;;]]- weights]]

MatrixForm/@{D2PS,dL}

Max[Abs[dL[[-dim;;,-dim;;]]-D2PS]]

Definition of 𝒟 and basis functions

We want to check relation () and approximation with Gaussian integrals.

fBase[k_,x_]:= Exp[-x/2] x LaguerreL[k-1,x];

For the Gauss Lobatto grid

ordLag = 1;
dim=3;
{xgrid,weights} = laguerreRoots[dim,ordLag];

fB = Table[fBase[i,xgrid],{i,1,dim}];
S = Transpose[fB];
S2 = Transpose[Table[Expand[(∂_(x,x))fBase[k,x]]/.x-> xgrid,{k,1,dim}]];

{rgrid,D2PS}=laguerredif[dim,2,ordLag];

MatrixForm/@{S2.Inverse[S],D2PS}

For the Gauss grid

ordLag = 0;
dim=3;
{xgrid,weights} = laguerreRoots[dim,ordLag];

fB = Table[fBase[i,xgrid],{i,1,dim}];
S = Transpose[fB];
S2 = Transpose[Table[Expand[(∂_(x,x))fBase[k,x]]/.x-> xgrid,{k,1,dim}]];

{rgrid,D2PS}=laguerredif[dim,2,ordLag];

MatrixForm/@{S2.Inverse[S],D2PS}

Gaussian integration

Let us show the correspondence between the Gaussian integration of matrix elements and the exact ones for the spectral derivative and for norm' s matrix matN.

{D2,matN} = laguerreMatrices[dim];

For the Gauss Lobatto grid

ordLag = 1;
{xgrid,weights} = laguerreRoots[dim,ordLag];

fB = Table[fBase[i,xgrid],{i,1,dim}];
S = Transpose[fB];
S2 = Transpose[Table[Expand[(∂_(x,x))fBase[k,x]]/.x-> xgrid,{k,1,dim}]];
wS = weights*Exp[xgrid]*S;  (* this is the matrix (Subscript[W, a]) (Subscript[S, ai ])*)
matNGauss = Transpose[wS].S;
matDGauss = Transpose[wS].S2;

MatrixForm/@{matNGauss,matN}

MatrixForm/@{Chop[matDGauss],N@D2}

For the Gauss  grid

ordLag = 0;
{xgrid,weights} = laguerreRoots[dim,ordLag];

fB = Table[fBase[i,xgrid],{i,1,dim}];
S = Transpose[fB];
S2 = Transpose[Table[Expand[(∂_(x,x))fBase[k,x]]/.x-> xgrid,{k,1,dim}]];
wS = weights*Exp[xgrid]*S;  (* this is the matrix (Subscript[W, a]) (Subscript[S, ai ])*)
matNGauss = Transpose[wS].S;
matDGauss = Transpose[wS].S2;

MatrixForm/@{matNGauss,matN}

MatrixForm/@{Chop[matDGauss],N@D2}

We see that in this case last element is only approximately reproduced.

Stability and precision

ordLag = 1;
Do[
{r,D2PS} = laguerredif[dim,2,ordLag];
{r,D2PShp} = laguerredif[dim,2,ordLag,30];
Print[{Norm[Flatten[D2PS-D2PShp]]/Norm[Flatten[D2PShp]],Precision[D2PShp]}],
{dim,50,-150,50}]

ordLag = 0;
Do[
{r,D2PS} = laguerredif[dim,2,ordLag];
{r,D2PShp} = laguerredif[dim,2,ordLag,30];
Print[{Norm[Flatten[D2PS-D2PShp]]/Norm[Flatten[D2PShp]],Precision[D2PShp]}],
{dim,50,-150,50}]

Example: Hydrogen atom

Clear[x,r];

We start with a relatively low order calculation, with dim = 20.

dim=20;

V[r_]:= - 1/r;

This function computes energy levels. The inputs are:

energyPS[λ_?NumberQ,{state_,L_},V_,xgrid_]:= Module[{matV,dim},
dim = Length[xgrid];
matV=DiagonalMatrix[V[λ xgrid]+(1/(2(λ^2)))L(L+1)(1/(xgrid^2))];
Sort[Re[Eigenvalues[((- 1)/(2(λ^2))) D2PS+matV ]]][[state]]]

Here we make the choice for the states and angular momentum

L=0; state=Range[1,5];

We want to show a slight difference between the Gaussian grid and Gauss - Lobatto grid. We compute for several λ energies E[λ]. In a variational scheme the correct value must be a minimum, then E[λ] cannot be lower than this value.

If the Gaussian integration is correct we must observe this kind of behavior also in our pseudospectral computation. The reader can check from figures below that this is indeed the case for the Gauss Lobatto grid, while for the Gauss grid at large λ the curve E[λ] bends toward lower values.

These values of λ are unphysical as we must look for a stability region for λ, in principle eigenvalues must be λ independent. In effect these "dangerous zones" move toward + ∞ as the dimension of the base grows, but, as in the spectral case, it is better to check that λ values used in the computation are in a reasonable range.

Energy vs λ for the Gauss  grid

ordLag=0;
{xgrid,D2PS} = laguerredif[dim,2,ordLag];

λPlot=Join[ Range[0.05,1,0.05],Range[1.25,10,0.25]];
enePlot = energyPS[#1,{state,L},V,xgrid]&/@λPlot;

ListPlot[Table[{λPlot[[k]],enePlot[[k,p]]},{p,state},{k,1,Length[λPlot]}],Joined-> True,PlotRange-> {{0,10},{-0.7,0}}]

Energy vs λ for the Gauss Lobatto grid

ordLag=1;
{xgrid,D2PS} = laguerredif[dim,2,ordLag];

λPlot=Join[ Range[0.05,1,0.05],Range[1.25,10,0.25]];
enePlot = energyPS[#1,{state,L},V,xgrid]&/@λPlot;

ListPlot[Table[{λPlot[[k]],enePlot[[k,p]]},{p,state},{k,1,Length[λPlot]}],Joined-> True,PlotRange-> {{0,10},{-0.6,0}}]

Fixing λ

Now we use FindMinimum to fix λ value and obtain an estimation of the eigenvalues. As the reader can verify the precision of the result is even embarrassing for such a low dimensional basis.

From now on we work with the Gauss Lobatto grid, the reader can check that all results are also obtained with the Gauss grid.

ordLag=1;
{xgrid,D2PS} = laguerredif[dim,2,ordLag];

eneLambda=Table[min=FindMinimum[
energyPS[λvar,{k,L},V,xgrid],{λvar,0.005, 3.0}];{min[[1]],λvar/.min[[2]]},{k,state}];

errors = Table[Abs[eneLambda[[k,1]]+1/(2 (((k+L))^2))],{k,state}];

TableForm[Table[{NumberForm[eneLambda[[k,1]],15],errors[[k]],eneLambda[[k,2]]},
{k,state}],TableHeadings-> {state,{"         E", "  |E - Th|", "   λ"}},
TableSpacing-> {3,3}]

Here we compute eigenvectors, using previous λ values

enePsi=Table[
λ=eneLambda[[k,2]];
 matV=DiagonalMatrix[V[λ xgrid]+(1/(2(λ^2)))L(L+1)(1/(xgrid^2))];
Ham = ((- 1)/(2(λ^2))) D2PS+matV;
{eigH,eigV}= Eigensystem[Ham];
order=Ordering[eigH];
eigH=eigH[[order]];
eigV=eigV[[order]];
{eneLambda[[k,1]],eigV[[k]]},{k,state}];

Let us test for example the wave function for level 5.

lev=5; 
λ=eneLambda[[lev,2]];

We extract the eigenvector from the matrix enePsi.

valPsi = enePsi[[lev,2]];

Define the "true" radial coordinates

rgrid= xgrid*λ;

and normalize the state

{temp,weights}=laguerreRoots[dim,ordLag];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi= valPsi/normPsi*Sign[valPsi[[-3]]];

As a first attempt we try to use Interpolation function of Mathematica to compute ψ at points different from the initial grid:

psi=Interpolation[Transpose[{rgrid,valPsi}]]

Plot[psi[x],{x,0.001,Max[rgrid]},PlotRange-> All]

We compare the result with the known reduced radial functions :

Rreduced[r_, n_, L_] := r*  (2/(n^2))(Sqrt[(((n-L-1)!)/((n+L)!))])*
(((((2 r)/n)))^L) * (E^(-(r/n)))* LaguerreL[n-L-1, 2 L + 1,((2 r)/n)] ;

normTh = Table[NIntegrate[((Rreduced[r,k,0])^2),{r,0,∞}],{k,state}]

Plot[Rreduced[r,lev,0],{r,0,Max[rgrid]},
PlotRange-> All]

Here is the maximum difference. We choose to compare points up Max[r-grid]/2, for higher values interpolation can be quite instable. This gives a quite large value for the fifth level:

rmax=λ Max[rgrid]/2

Let us compare now the interpolated function and exact function at the grid points and at a uniform grid up to rmax :

rsample = Range[0.1,rmax,rmax/100];

{Max[Abs[psi[rgrid]-Rreduced[rgrid,lev,0]]],Max[Abs[psi[rsample]-Rreduced[rsample,lev,0]]]}

It is seen that while the values computed on the original grid are very accurate, the extrapolated ones are quite inaccurate. This is why in the next section  the function polint, or better interpLaguerre, will be used to get a better interpolation, at least at relatively low dimensions.

Here is an example of the above estimate using interpLaguerre

Polint

state = Range[5];

{xgrid,weights}=laguerreRoots[dim,ordLag];
δψStates=Table[
λ=eneLambda[[lev,2]];
valPsi= enePsi[[lev,2]];
normPsi = Sqrt[λ Tr[weights*Exp[xgrid]*(valPsi^2)]];
valPsi = valPsi/normPsi;
rgrid = λ xgrid;
valTh=Rreduced[rgrid,lev,0]*Sign[ Rreduced[rgrid[[3]],lev,0]];
δψOnGrid = Max[Abs[valPsi*Sign[valPsi[[3]]]-valTh]];
rmax=Max[λ xgrid]/2;
r = Range[Min[xgrid]λ,rmax,rmax/100];
fPsi= interpLaguerre[{xgrid,valPsi},λ,r];
sign1 = Sign[fPsi[[3]]];
valTh = Rreduced[r,lev,0]*Sign[ Rreduced[r[[3]],lev,0]];
δψ = Max[Abs[fPsi*sign1-valTh]];
{δψOnGrid,δψ},{lev,state}];

TableForm[δψStates,TableHeadings-> {state,{" on grid"," off grid"}},TableSpacing-> {3,3}]

Systematic analysis

Here we follow the lines of notebook [7] and analyze the results varying basis's dimension. In this section we compare results directly to the known analytical results.

Off[Power::"infy",InterpolatingFunction::"dmval"]
ordLag = 1;
L=0; state=Range[5];
eneTh = -1/2/(state+L)^2;

This is the interval we choose for λ. In more complicated cases a level dependent interval can be defined.

λInt = {0.005,3.0};

dimRange =Prepend[Range[20,90,20],10];

To make clear the procedure we adopt an inefficient procedure, computing several times the same thing, as in the case of xgrid. This allow to separate the different quantities and it is intended for helping the reader.

Here we compute the couples (E, λ) for each level. This is the most time consuming part of the computation.

eneAndLambda=Table[
{xgrid,D2PS} = laguerredif[dim,2,ordLag];
Table[min=
FindMinimum[energyPS[λvar,{k,L},V,xgrid],{λvar,λInt[[1]],λInt[[2]]}];
{min[[1]],λvar/.min[[2]]},{k,state}],
{dim,dimRange}];

eneAndLambda has 3 indexes, energies and λ's are stored in the form

eneAndLambda[[ dimension, level, 1]] → energy
eneAndLambda[[ dimension, level, 2]] → λ

These are the corresponding eigenstates. The format is

eigenstates[[ dimension, level]] → eigenstate (a vector)

eigenstates=Table[
dim = dimRange[[j]];
{temp,weights}=laguerreRoots[dim,ordLag];
{xgrid,D2PS} = laguerredif[dim,2,ordLag];
Table[λ = eneAndLambda[[j,k,2]];
 matV=DiagonalMatrix[V[λ xgrid]+(1/(2(λ^2)))L(L+1)(1/(xgrid^2))];
Ham = ((- 1)/(2(λ^2))) D2PS+matV;
{eigH,eigV}= Eigensystem[Ham];
order=Ordering[eigH];
eigH=eigH[[order]];
eigV=eigV[[order]];
cf = eigV[[k]],{k,state}],
{j,1,Length[dimRange]}];

The relative errors for energies are :

deltaE= Table[Abs[(eneAndLambda[[j,All,1]]-eneTh)/eneTh],{j,1,Length[dimRange]}];

deltaψOnGrid contains the maximum error of ψ on grid points, the format is the same as eigenstates.

deltaψOnGrid = Table[
dim = dimRange[[j]];
{xgrid,weights}=laguerreRoots[dim,ordLag];
Table[
λ = eneAndLambda[[j,lev,2]];
valPsi= eigenstates[[j,lev]];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi=valPsi/normPsi*Sign[valPsi[[3]]];
valTh = Rreduced[λ xgrid,lev,0]*Sign[ Rreduced[λ xgrid[[3]],lev,L]];
Max[Abs[valPsi-valTh]]
,{lev,state}],{j,1,Length[dimRange]}];

deltaψ are errors on the wave function computed with interpLaguerre. This function needs as input the pair  {xgrid,ψgrid}, i.e. initial grid and function to be interpolated on this points; λ, the scale parameter, and r, the points  on the original non scaled grid, where we want the value of the function. Up to dimensions less than 45 a mixture of polint and Mathematica function Interpolation is used, for dim > 45 is used Interpolation alone, in this case errors grow but polint function appears to be instable. The reader can explore different kind of interpolations.

deltaψ= Table[
dim = dimRange[[j]];
{xgrid,weights}=laguerreRoots[dim,ordLag];
Table[
λ = eneAndLambda[[j,lev,2]];
valPsi= eigenstates[[j,lev]];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi=valPsi/normPsi; 
rmax = 5 (lev^2+1);
r = Range[0.01,0.9rmax,rmax/200];
fPsi= interpLaguerre[{xgrid,valPsi},λ,r];
sign1 = Sign[fPsi[[3]]];
valTh = Rreduced[r,lev,0]*Sign[ Rreduced[r[[3]],lev,L]];
Max[Abs[fPsi*sign1-valTh]]
,{lev,state}],{j,1,Length[dimRange]}];

Here we present the results

TableForm[deltaE,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}]

ScientificForm[TableForm[deltaψ,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}],3]

ScientificForm[TableForm[deltaψOnGrid,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}],3]

The general picture which emerges is that it is quite useless to go beyond dim = 30 in this case, the only effect is a growing of rounding errors. Clearly the situation can be different for higher states, and it will be different for other potentials.

Second Example

In this section we repeat step by step the previous analysis.  The aim is to show a potential for which Gaussian integration is not exact. Errors will be estimated from the computation and compared with the analytical results. The procedure follows closely notebook [7] where the same computation was made within spectral scheme.

Model

Here is the model, its eigenvalues and eigenstates.

V0 = 10;V[x_]:= - (V0/((Cosh[x])^2));

sTh =1/2*((Sqrt[1+ 8 V0])  -1);
nE = Select[Range[1,11,2],#≤ sTh &]

eneTh = -(1/8)((((Sqrt[1+ 8 V0]) - (2 nE +1)))^2)

psiLL[n_,x_]:= (((1-((Tanh[x])^2)))^((sTh-n)/2))Hypergeometric2F1[-n,2sTh -n+1,sTh-n+1,((1-Tanh[x])/2)];
normaTh = Table[NIntegrate[psiLL[nE[[i]],x]^2,{x,0,Infinity}],{i,1,Length[nE]}];

psiTh[n_,x_]:= psiLL[nE[[n]],x]/(Sqrt[normaTh[[n]]]);

Preliminary analysis

We plot energy levels as a function of scale parameter λ to have an idea of possible problems. The computation will be done for several dimensions to point out the changes in the graphs.

L=0; ordLag=1;
state=Range[Length[eneTh]];

λGrid = Join[Range[0.05,0.9,0.05],Range[1,10,0.25]];

dimRange=Range[10,60,10];

This is the computation for ground state

ene1=Table[{xgrid,D2PS} = laguerredif[dim,2,ordLag];
{energyPS[#1,{1,L},V,xgrid]&/@λGrid,energyPS[#1,{2,L},V,xgrid]&/@λGrid},
{dim,dimRange}];

Table[ListPlot[Transpose[{λGrid,ene1[[k,1]]}],Joined-> True,AxesOrigin-> {0,-6}],{k,1,Length[dimRange]}]

This is the computation for the excited state

Table[ListPlot[Transpose[{λGrid,ene1[[k,2]]}],Joined-> True,AxesOrigin-> {0,0}],{k,1,Length[dimRange]}]

It is seen that now the situation can be worse than in the Coulomb case. Fake minima appear:  they can be identified by their shift toward λ = ∞ as the dimension grows. A stable range for small λ exists and grows with the dimension: one must be sure that selected λ values are within this "safe" region.

Systematic analysis

Now we proceed as in the Coulomb case.  We just copy the above lines for computing energies and eigenstates.

dimRange =Range[20,50,10];
dimRange={10,20,30,40,45,50,60,70};

Off[Power::"infy",InterpolatingFunction::"dmval"]

We choose a wrong interval for λ to show that wrong results can be obtained

λInt = {0.005,3.0};

eneAndLambda=Table[
{xgrid,D2PS} = laguerredif[dim,2,ordLag];
Table[min=
FindMinimum[energyPS[λvar,{k,L},V,xgrid],{λvar,λInt[[1]],λInt[[2]]}];
{min[[1]],λvar/.min[[2]]},{k,state}],
{dim,dimRange}];

These are the corresponding eigenstates.

eigenstates=Table[
dim = dimRange[[j]];
{temp,weights}=laguerreRoots[dim,ordLag];
{xgrid,D2PS} = laguerredif[dim,2,ordLag];
Table[λ = eneAndLambda[[j,k,2]];
 matV=DiagonalMatrix[V[λ xgrid]+(1/(2(λ^2)))L(L+1)(1/(xgrid^2))];
Ham = ((- 1)/(2(λ^2))) D2PS+matV;
{eigH,eigV}= Eigensystem[Ham];
order=Ordering[eigH];
eigH=eigH[[order]];
eigV=eigV[[order]];
cf = eigV[[k]],{k,state}],
{j,1,Length[dimRange]}];

deltaE= Table[Abs[(eneAndLambda[[j,All,1]]-eneTh)/eneTh],{j,1,Length[dimRange]}];

deltaψOnGrid contain the maximum error of ψ on grid points, the format is the same as eigenstates.

deltaψOnGrid = Table[
dim = dimRange[[j]];
{xgrid,weights}=laguerreRoots[dim,ordLag];
Table[
λ = eneAndLambda[[j,lev,2]];
valPsi= eigenstates[[j,lev]];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi=valPsi/normPsi*Sign[valPsi[[3]]];
valTh = psiTh[lev,λ xgrid]*Sign[ psiTh[lev,λ xgrid[[3]]]];
Max[Abs[valPsi-valTh]]
,{lev,state}],{j,1,Length[dimRange]}];

deltaψ are errors on the wave function computed with interpLaguerre. This function needs as input the couple {xgrid,ψgrid}, i.e. initial grid and function to be interpolated on this points; λ, the scale parameter, and r, the points  on the original non scaled grid, where we want to compute the function. Up to dimensions less than 45 a mixture of polint and Mathematica function Interpolation is used, for dim > 45 is used Interpolation alone, in this case errors grow but polint function appears to be instable. The reader can explore different kind of interpolations.

deltaψ= Table[
dim = dimRange[[j]];
{xgrid,weights}=laguerreRoots[dim,ordLag];
Table[
λ = eneAndLambda[[j,lev,2]];
valPsi= eigenstates[[j,lev]];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi=valPsi/normPsi; 
rmax = λ Max[xgrid]/2;
r = Range[0.01,0.9rmax,rmax/80];
fPsi= interpLaguerre[{xgrid,valPsi},λ,r];
sign1 = Sign[fPsi[[3]]];
valTh =psiTh[lev,r]*Sign[ psiTh[lev,r[[3]]]];
Max[Abs[fPsi*sign1-valTh]]
,{lev,state}],{j,1,Length[dimRange]}];

Here we present the results

TableForm[Table[{deltaE[[j,1]],eneAndLambda[[j,1,2]],
deltaE[[j,2]],eneAndLambda[[j,2,2]]},{j,1,Length[dimRange]}],
TableHeadings-> {dimRange,
{" δE[1]", " λ" ," δE[2]", " λ"}},
TableSpacing-> {3,3}]

The result for dim = 45 is completely wrong, and this corresponds to the wrong λ value in the last column. The same thing happens for dim=10, but in this case one can think that the basis is too small.

TableForm[
{{ScientificForm[TableForm[deltaψ,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}],3],
ScientificForm[TableForm[deltaψOnGrid,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}],3]}},TableSpacing-> {5,5},
TableHeadings-> {None,{"       δψ", "     δψ on grid"}}]

Here is the computation with a more reasonable λ interval

λInt = {0.02,1.0};

eneAndLambda=Table[
{xgrid,D2PS} = laguerredif[dim,2,ordLag];
Table[min=
FindMinimum[energyPS[λvar,{k,L},V,xgrid],{λvar,λInt[[1]],λInt[[2]]}];
{min[[1]],λvar/.min[[2]]},{k,state}],
{dim,dimRange}];

These are the corresponding eigenstates.

eigenstates=Table[
dim = dimRange[[j]];
{temp,weights}=laguerreRoots[dim,ordLag];
{xgrid,D2PS} = laguerredif[dim,2,ordLag];
Table[λ = eneAndLambda[[j,k,2]];
 matV=DiagonalMatrix[V[λ xgrid]+(1/(2(λ^2)))L(L+1)(1/(xgrid^2))];
Ham = ((- 1)/(2(λ^2))) D2PS+matV;
{eigH,eigV}= Eigensystem[Ham];
order=Ordering[eigH];
eigH=eigH[[order]];
eigV=eigV[[order]];
cf = eigV[[k]],{k,state}],
{j,1,Length[dimRange]}];

deltaE= Table[Abs[(eneAndLambda[[j,All,1]]-eneTh)/eneTh],{j,1,Length[dimRange]}];

deltaψOnGrid contain the maximum error of ψ on grid points, the format is the same as eigenstates.

deltaψOnGrid = Table[
dim = dimRange[[j]];
{xgrid,weights}=laguerreRoots[dim,ordLag];
Table[
λ = eneAndLambda[[j,lev,2]];
valPsi= eigenstates[[j,lev]];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi=valPsi/normPsi*Sign[valPsi[[3]]];
valTh = psiTh[lev,λ xgrid]*Sign[ psiTh[lev,λ xgrid[[3]]]];
Max[Abs[valPsi-valTh]]
,{lev,state}],{j,1,Length[dimRange]}];

deltaψ are errors on the wave function computed with interpLaguerre. This function needs as input the couple {xgrid,ψgrid}, i.e. initial grid and function to be interpolated on this points; λ, the scale parameter, and r, the points  on the original non scaled grid, where we want to compute the function. Up to dimensions less than 45 a mixture of polint and Mathematica function Interpolation is used, for dim > 45 is used Interpolation alone, in this case errors grow but polint function appears to be instable. The reader can explore different kind of interpolations.

deltaψ= Table[
dim = dimRange[[j]];
{xgrid,weights}=laguerreRoots[dim,ordLag];
Table[
λ = eneAndLambda[[j,lev,2]];
valPsi= eigenstates[[j,lev]];
normPsi= (Sqrt[λ Tr[ (valPsi^2)*weights*Exp[xgrid]]]);
valPsi=valPsi/normPsi; 
rmax = λ Max[xgrid]/2;
r = Range[0.01,0.9rmax,rmax/80];
fPsi= interpLaguerre[{xgrid,valPsi},λ,r];
sign1 = Sign[fPsi[[3]]];
valTh =psiTh[lev,r]*Sign[ psiTh[lev,r[[3]]]];
Max[Abs[fPsi*sign1-valTh]]
,{lev,state}],{j,1,Length[dimRange]}];

Here we present the results

TableForm[Table[{deltaE[[j,1]],eneAndLambda[[j,1,2]],
deltaE[[j,2]],eneAndLambda[[j,2,2]]},{j,1,Length[dimRange]}],
TableHeadings-> {dimRange,
{" δE[1]", " λ" ," δE[2]", " λ"}},
TableSpacing-> {3,3}]

It is seen that the result for dim = 45 is correct now.

TableForm[
{{ScientificForm[TableForm[deltaψ,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}],3],
ScientificForm[TableForm[deltaψOnGrid,TableHeadings-> {dimRange,Table["n = "<>ToString[k],{k,state}]},
TableSpacing-> {3,3}],3]}},TableSpacing-> {5,5},
TableHeadings-> {None,{"       δψ", "     δψ on grid"}}]

Warning : for dim > 45 errors on ψ can depend quite critically on the interval in λ, in similar cases it is better to check stability varying the interval.

References

1) J.A.C. Weideman and S.C. Reddy, A MatLab Diferentiation Suite, ACM Trans. on Math. Software, Vol. 26, 465, (2000).