Laguerre basis: examples

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the stability and the precision of the variational results using a basis of Laguerre polynomials. In particular compute numerically with 20 digit precision the ground state energy for a radial potential  - (Subscript[V, 0])/ Cosh[α ((r])^2).

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

Introduction

The general problem

In this notebook we shall study in some detail the variational calculation of eigenvalues and eigenfunctions for radial problems using Laguerre polynomials as the basis functions. The typical problem is

-((ℏ^2)/(2m)) (∂_(r,r))f[r] + (V[r]+ (((ℏ^2)L(L+1))/(2m (r^2))) )f[r] = E f[r];   f[0] = 0.

i.e. the radial Schrödinger equation for the reduced wave function, the radial part of the wave function being R = f/r. We shill limit ourselves to bound state problems. The solution of equation () can be expressed through a linear combination of the complete set (for f[r] vanishing at the origin):

(Subscript[b, k])[r]=Exp[-r/2]r L[k-1,r]; k=1,2...

L[k,r] stands for the k-th Laguerre polynomial. The normalization used in this book is the one used in Mathematica,

(∫_0^∞)Exp[-r] L[k,r]L[j,r] = (Subscript[δ, ij])

From now on the natural units will be used, i.e.,  ℏ = m = 1.

Variational calculations

A simple  approach to the solution of problem () is based on a variational calculation: a finite subset of the basis () is equipped with a scale parameter λ

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

and used to approximate the variational principle

δ ((⟨ψ|H|ψ⟩)/(⟨ψ|ψ⟩)) = 0

Matrix elements of kinetic and potential terms are computed on the basis (). The parameter λ can be reabsorbed by the change of variables

r = λ x

Homogeneous functions, like derivatives, are related by a simple power of λ to the same quantities in the r variable.

N.B. The symbol r will be used for the original space variable, the symbol x for the rescaled variable.

The matrix equation  in the rescaled variables takes the form:

-(1/(2(λ^2))) D2 .ψ+ matV[λ x ].ψ = ℰ N.ψ

The matrix N in the right hand side of () is due to non-orthogonality of the basis (), and will be denoted by matN in the computations. The matrix D2 will denote matrix elements of the second derivative

(Subscript[D2, ij]) = (∫_0^∞)𝕕x (Subscript[f, i])[x] (∂_(x.x))(Subscript[f, j])[x].

The matrix matV is given in the rescaled variables, by

(Subscript[matV, ij]) = (∫_0^∞)𝕕x (Subscript[f, i])[x]V[λ x] (Subscript[f, j])[x].

Computation of matrix elements

Matrix elements for D2 and N are computed by the function laguerreMatrices using the recurrence relations of Laguerre polynomials, see notebook [].

Matrix elements () for the potential are approximated by gaussian integration formulas, see notebook [NB-10.3-Gauss.nb]:

(∫_0^∞)Exp[-x]F[x] 𝕕x ≃ (∑_(a=1)^n)(Subscript[W, a]) F[(Subscript[x, a])]

The n abscissas (Subscript[x, a])are the roots of the n-th Laguerre polynomial, (Subscript[W, a]) the corresponding weights. Both quantities are computed by the function laguerreRoots.  Note the following point:

The weight function Exp[-x] has been explicitly factorized.  If we apply this formula to the matrix elements () one  must write only the polynomial part of the basis functions. Explicitly:

(Subscript[matV, ij]) ≃ (∑_(a=1)^n)(Subscript[W, a]) ( L[i-1,(Subscript[x, a])]L[j-1,(Subscript[x, a])] V[λ (Subscript[x, a])]).

In the following section we collect some common functions that will be used throughout the notebook with a brief description of their functionality.

Functions

laguerreMatrices

Compute by default matrices D2 and matN. An optional list containing powers of r greater 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 gaussian integration. An optional argument, prec, allow to compute the same quantities in precision prec arithmetic. The standard calling sequence will be

{xgrid, weights} = laguerreRoots[n]

laguerreRoots[n_,prec_:MachinePrecision]:=Module[{J,δ,γ,radEig,w},
δ = 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]]];
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 needs: 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 compute eigenvalues with a call to Mathematica function Eigenvalues. It is written in order to compute also 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]]]

Wave functions

The basis is defined by

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

It useful to have at disposal a tool which compute quickly wave functions, for plotting  or other needs. The analytical definition requires the computation of a Laguerre polynomial for each x, or if used in a vectorial form for a grid of points, a call to several Laguerre polynomials for the grid. It is simpler to define Laguerre polynomials via their recursion relation

x L[n,x] = - (n+1) L[n+1,x] + (2n+1) L[n,x] - n L[n-1,x];

This function, written as a compiled function for efficiency, does this work :

recursiveLaguerrePolynomials= Compile[{{n,_Integer},{x,_Real,1}},
Module[{y1,y0,p1,ris,ris0,c},
y0 = 1.0 + 0.*x;
y1 = 1.-x;
ris0={y0,y1};
ris=Table[
c=1./( 1.k+1.);
p1=c(2 k+1. - x) y1 - c*k y0;
y0=y1;
y1=p1;
p1,{k,1,n-2}];
ris=Join[ris0,ris];
ris]];

It computes the first n polynomials, from degree 0 to n - 1. The following is simply an interface which simplifies the output when x is a single coordinate instead of a grid of points. This is the only function which the reader might possibly need for Laguerre polynomials.

laguerreP[dim_,xIn_]:= Module[{x,ris},
x = Flatten[{xIn}];
ris= recursiveLaguerrePolynomials[dim,x];
If[Length[xIn]== 0,Flatten[ris],ris]]

For computing a wave function one has to multiply the coefficient' s vector by a sequence of basis functions. This is done simply with:

getPsi[vec_,r_]:= (vec.laguerreP[Length[vec],r])Exp[-r/2]*r

One must provide, if necessary, the normalization of the vector and the possible scale factor λ.

Example

Model

As a case study we will take the potential V[r] = - (Subscript[V, 0])/Cosh[α ((r])^2). After a change of variables the constant α can be reabsorbed in r and we consider as a standard form

V[r_]:= - (V0/((Cosh[ r])^2));

Plot[V[r]/.{V0-> 10,α-> 1},{r,0,3},AxesLabel-> {"r",None}]

The potential can be resolved in the one-dimensional case, its interpretation as radial potential amounts to select odd one-dimensional states, i.e. those with reduced wave function vanishing at r=0.  We will take V0 = 10 for example.

Eigenvalues and eigenstates are found in the book of Landau and Lifshitz [] and are given by (we select odd values of the principal quantum number, i.e. odd states). For (Subscript[V, 0])=10 two bound states exist

V0 = 10;
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_,r_]:= psiLL[nE[[n]],r]/(Sqrt[normaTh[[n]]]);

Plot[{psiTh[1,r],psiTh[2,r]},{r,0,6}]

Computation

In this section  a "blind" computation for the eigenvalues and eigenfunctions is done.  This procedure is not to be considered as a definitive procedure because several potential errors can be inadvertently introduced.  The next two sections will be devoted to an analysis of the results.

This is our set of basis functions

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

To start with, we fix a dimension,  say 30

dim=30;

and compute derivatives and the matrix of scalar products, matN :

{D2,matN}=N@laguerreMatrices[dim];

For possible studies of higher angular momenta we remember that in the basis we are using the operator 1/(r^2) is represented by the identity matrix.

Here are the abscissas and weights for Gauss formulas

extraDim = 4;

{xgrid,weights} =laguerreRoots[dim+extraDim];

N.B. The argument of laguerreRoots is independent of the dimension used, even if a number of points greater than dim is strongly encouraged. To have some extra precision we took as number of abscissas dim + 4, but the reader can make freely another choice, keeping attention to possible instabilities for too large number of abscissas.

This matrix contains the values of the basis functions, without the exponential part,  on the points xgrid

matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];

Eigenvalues

Let us choose quantum numbers for the states we want to study

state = {1,2}; L=0;

and do a simple check of the output of the function energy, defined above, with an arbitrary parameter λ = 1 (the natural scale of the problem)

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

We try a variational determination for λ, in the region of the natural scale.

N.B. This is a very dangerous assumption as will be discussed in the next section.

eneAndLambda=Table[min=FindMinimum[energy[λvar,{k,L},V,{xgrid,weights},matBaseLag],{λvar,0.01,1.2}];{min[[1]],λvar/.min[[2]]},{k,1,2}]

NumberForm[eneAndLambda[[All,1]],20]

errorsThdim30=eneAndLambda[[All,1]]-eneTh

errorsThdim30=eneAndLambda[[All,1]]-eneTh

Eigenfunctions

Let us now have a look at the eigenfunctions.

state=1;

We store for later uses  the value of λ computed before :

λ = eneAndLambda[[state,2]];
λGSdim30 = λ

The eigenfunctions are expressed as linear combination of basis functions. It must be remembered that r - variable (i.e. usual space variable) is obtained from rescaled x value by r = λ x, then

ψ[r] = (∑_(k=1)^dim)(Subscript[c, k]) (Subscript[f, k])[λ x]

The vectors of different possible (Subscript[c, k])coefficients are  the eigenvectors of the Hamiltonian matrix, i.e. the solutions of

-(1/(2(λ^2)))(Subscript[D2, ij]) (Subscript[c, j]) + (Subscript[matV, ij]) (Subscript[c, j])= E (Subscript[N, ij])(Subscript[c, j])

We use Eigensystem to compute these eigenvectors (for L≠0 states on must add the 1/r^2 potential term)

matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
Ham = ((- 1)/(2(λ^2))) D2+matV;
{eigH,eigV}= Eigensystem[{Ham,matN}];
order = Ordering[Re[eigH]];
eigH = eigH[[order]];
eigV = eigV[[order]];

ene = eigH[[state]];
cf = eigV[[state]];

The norm of the states is given by

((‖ψ‖)^2) = (∫_0^∞)ψ[r]dr =λ (∫_0^∞)ψ[x]dx = λ(∑_ij) (Subscript[c, i])(Subscript[c, j]) (∫_0^∞)(Subscript[f, i])[x](Subscript[f, j])[x]dx = λ (Subscript[c, i])(Subscript[N, ij])(Subscript[c, j])

normψ = Sqrt[λ*cf.matN.cf];

Then the normalized eigenfunctions are given by

ψ[r_]:=( cf.fBase[Range[dim],r/λ])/normψ

As a check :

NIntegrate[((ψ[r])^2),{r,0,∞}]

One can use the alternative function getPsi

The two functions are indistinct on a graph. In other cases some change in sign can be needed.

Plot[{ψ[r],psiTh[state,r]},{r,0,3}]

The functions coincide up to O((10^(-8))), which is a good precision for this kind of computations.

LogPlot[Abs[ψ[r]-psiTh[state,r]],{r,0,4},PlotRange-> All]

Error estimates

In this section we want to give some (partial) answers to the following questions:

Is our interval for searching λ optimal value reasonable or completely wrong?

In general, without knowing the analytical answers, how can we estimate the  errors in the results?

λ range

Let us consider the eigenvalues variation with λ in the system just considered

state={1,2}; L=0;

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

ListPlot[Table[{λPlot[[k]],ene[[k,p]]},{p,1,2},{k,1,Length[λPlot]}]]

We see where the problem lies : if the functions E[λ] are not convex a blind choice of λmin ≤ λ ≤ λmax in FindMinimum can easily select a spurious value for λ.

A point must be clear to the reader, a scale like λ is quite special as a variational parameter. In effect we know in advance that in the infinite dimensional Hilbert space eigenvalues cannot depend on λ, its introduction is just a change of variables. This means that "good" regions for λ are those in which E[λ] is nearly constant. In our parametrization wrong values of λ tend to be at higher scales, had we adopted the opposite choice, i.e. r = x/α, all this nasty stuff would be concentrated in the small region α<1 and would be easily overlooked.

The situation clearly goes better for higher dimensions, but always spurious minima prevent us to make a priori choice of λ interval. Let us make a quick scan for 10 < dim < 50, for each eigenvalue

λPlot = Join[Range[0.05,0.9,0.05],Range[1,8,0.25]];
dimRange=Range[10,60,10];

ene1=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
energy[#1,{state,L},V,{xgrid,weights},matBaseLag]&/@λPlot,
{dim,dimRange}];

Table[ListPlot[Transpose[{λPlot,ene1[[k,All,1]]}],Joined-> True,AxesOrigin-> {0,-5},PlotLabel-> Row[{"eig: ",1,";  dim= ",dimRange[[k]]}]],{k,1,Length[dimRange]}]

Table[ListPlot[Transpose[{λPlot,ene1[[k,All,2]]}],Joined-> True,AxesOrigin-> {0,-1},PlotLabel-> Row[{"eig: ",2,";  dim= ",dimRange[[k]]}]],{k,1,Length[dimRange]}]

We see that the stability region grows toward large λ as dim increases, but, especially for excited state a wrong minimum is present, at around λ = 4 for dim = 60.

Another danger for too high λ is due to the fact that in this limit the eigenvalues become completely different, they refer essentially to a free particle, this is true particularly for low dimensions, as the graph for dim=10 explicitly shows.

We have thus learned our first rule: a graphical inspection of λ dependence, even at low dimensions, can help us avoiding misinterpretation of the results, going to large basis and checking for stability is absolutely mandatory, low-dimensional results can be very poor.

A second rule is to avoid possible numerical instabilities. These can be of  different origins: a low dimensional grid for gaussian integration can cause poor results;  a too large grid can be unstable as the weights rapidly grow. This is particularly true for eigenfunctions, even the numerical evaluation of high degree Laguerre polynomials can give problems.

N.B. In the framework of the variational principle the function E[λ] cannot assume values less than the "true" minimum of the problem. The fake minima are due to the approximation of integrals with gaussian integration.

Error estimates

Now let us come to an estimation of the errors in the eigenvalues. The best thingto do is to compare the results as the dimension of the base is varied. The difference can give both an error estimate and a check of stability. This can be a time-consuming procedure. But if the result is to be reliable, it must be done.

We will use a simplified version of the Temple method. For the ground state the method tells us that an estimation of the error is

δE = (σ/((Subscript[E, 2])-(Subscript[E, 1])));           σ = ⟨ψ|(H^2)|ψ⟩ - (((⟨ψ|H|ψ⟩))^2)

E is the energy of the first excited state, σ is the quadratic fluctuation of energy on the ground state. Numerically we cannot compute σ on the same system where we compute (Subscript[E, 1]),  since  ψ in a finite dimensional system has no energy fluctuations, being an eigenstate of H.  One can give an estimate by growing dim, for instance  dim → dinNew = dim+4,  by identifying the old eigenvector with the first dim component of the state and computing on this state the fluctuation. An example will clarify the procedure. Let us do  the computation for dim = 30 from scratch.

This compute eigenvalues and eigenvectors.

Eigenvalues

state = {1,2}; L=0;
dim=30;
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
eneAndLambda=Table[min=FindMinimum[energy[λvar,{k,L},V,{xgrid,weights},matBaseLag],{λvar,0.01,1.2}];{min[[1]],λvar/.min[[2]]},{k,1,2}]

Eigenvectors

λ = eneAndLambda[[1,2]];
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
Ham = ((- 1)/(2(λ^2))) D2+matV;
{eigH,eigV}= Eigensystem[{Ham,matN}];
order = Ordering[Re[eigH]];
eigH = eigH[[order]];
eigV = eigV[[order]];
cf = eigV[[1]];

Immersion in a larger dimension

dimNew = dim+4;
{D2new,matNnew}=N@laguerreMatrices[dimNew];
{xgridnew,weightsnew} =laguerreRoots[dimNew+extraDim];
matBaseLagnew= Table[xgridnew*LaguerreL[k-1,xgridnew],{k,1,dimNew}];
matVnew=getMatV[λ,V,{xgridnew,weightsnew},matBaseLagnew];
Hamnew = ((- 1)/(2(λ^2))) D2new+matVnew;

This is the eigenvector immersed in the larger space

cfNew = PadRight[cf,dimNew];

To compute the quadratic fluctuation we have to remember that for a non-orthogonal basis operator products correspond to:

Op[a]Op[b] = matA.(matN^(-1)).matB

For σ and the estimation of the error we have

σ = cfNew.Hamnew.Inverse[matNnew].Hamnew.cfNew/cfNew.matNnew.cfNew - ((eneAndLambda[[1,1]])^2);

δETemple = σ/(eneAndLambda[[2,1]]-eneAndLambda[[1,1]])

To be compared with the true error

Abs[eneAndLambda[[1,1]]-eneTh[[1]]]

The agreement is excellent (we are speaking of orders of magnitude,  δETemple in this sense is a good error estimate).

Precision vs dimension

In this subsection we collect all previous observations and perform a systematic analysis of errors.

We select a set of dimensions, the states to be studied and a fiducial starting interval for λ search.

dimRange= Prepend[Range[20,120,20],10];
λInt = {0.01,1.8};
state={1,2};

We compute energies and relative λ, to check that λ is not in one of the dangerous zones shown in previous graphs.

eneAndLambda=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
Table[min=FindMinimum[energy[λvar,{k,L},V,{xgrid,weights},matBaseLag],{λ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]] → λ

From computed energies we can estimate errors comparing an eigenvalue at successive dimensions. The first entry for the matrix is absent evidently.

δE =Transpose[Table[ Drop[eneAndLambda[[All,k,1]],1]-Drop[eneAndLambda[[All,k,1]],-1],
{k,state}]];
δE = Prepend[Abs[δE],ConstantArray["*",Length[state]]];

δE has 2 indexes, errors are stored in the form

δE[[ dimension, level]] → δE

Finally we compute eigenvectors and estimate at the same time errors with Temple method.

δETemple is a matrix with the same format as δE.

eigenstates is a 2 indexes matrix, which store normalized eigenvectors, i.e. each eigenvectors is divided by the square root of v.matN.v, this is the form needed to compute normalized eigenfunctions. The format of eigenstates is

eigenstates[[ dimension, level]] → normalized coefficients of the eigenvectors

δETemple = ConstantArray[0,{Length[dimRange],Length[state]}];

eigenstates=Table[
dim = dimRange[[j]];
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k1-1,xgrid],{k1,1,dim}];
Table[λ = eneAndLambda[[j,k,2]];
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
Ham = ((- 1)/(2(λ^2))) D2+matV;
{eigH,eigV}= Eigensystem[{Ham,matN}];
order = Ordering[Re[eigH]];
eigH = eigH[[order]];
eigV = eigV[[order]];
cf = eigV[[k]];
normψ = Sqrt[λ*cf.matN.cf];
cfNormalized = cf/normψ;
(* immersion in a larger dimension *)
dimNew = dim+4;
{D2new,matNnew}=N@laguerreMatrices[dimNew];
{xgridnew,weightsnew} =laguerreRoots[dimNew+extraDim];
matBaseLagnew= Table[xgridnew*LaguerreL[k-1,xgridnew],{k,1,dimNew}];
matVnew=getMatV[λ,V,{xgridnew,weightsnew},matBaseLagnew];
Hamnew = ((- 1)/(2(λ^2))) D2new+matVnew;
cfNew = PadRight[cf,dimNew];
σ = cfNew.Hamnew.Inverse[matNnew].Hamnew.cfNew/cfNew.matNnew.cfNew - ((eneAndLambda[[j,k,1]])^2);
eneUp = If[k== Max[state],0,eneAndLambda[[j,k+1,1]]];
δETemple[[j,k]]= Abs[σ/(eneAndLambda[[j,k,1]]-eneUp)];
cfNormalized
,{k,state}],
{j,1,Length[dimRange]}];

Here are the results. We do not analyze the stability of wave functions, this will be done below by comparison with analytical results.

The results for the energy are printed with 15 decimal digits, this allows a quick comparison with δE in this simple case (we know that eigenvalues are -1/2, -9/2).

TableForm[
Table[TableForm[Table[Flatten[{NumberForm[eneAndLambda[[j,lev,2]],3],NumberForm[eneAndLambda[[j,lev,1]],16],
ScientificForm[#1,3]&/@ {δE[[j,lev]],δETemple[[j,lev]]} 
}],{j,1,Length[dimRange]}],
TableHeadings-> {dimRange,{ "  λ","       E["<>ToString[lev]<>"]","  δE ","δETemple "}},
TableSpacing-> {4,4}],{lev,state}],TableSpacing-> {10,5}]

Comments

We learn from previous table that errors δE and δETemple are compatible, even if the latter seems to be a better indicator (this will be seen again below). More importantly we see that precision do not seem to improve significantly for ground state after dim = 40 and for excited state after dim=60-80. This means that rounding errors come in the game, and if we want better precision we have to face this problem (see next section).

As a general lesson we have learned that the highest precision to be expected is of the order of 10^-13 and that after this precision has been attained rounding errors start to downgrade the results.

Analytical results

In this particular case we have analytical results and we can compare them with numerical ones. Here we use computed eigenvectors to estimate the precision of the wave function and write explicitly the difference between numerical and analytical eigenvalues. The errors on the wave functions have been computed on a grid spaced 0.1 between r=0 and r=30.

δEδψTh = Table[
Table[
dim = dimRange[[j]];
λ = eneAndLambda[[j,k,2]];
rgrid = Range[0,30,0.1];
ψ[r_]:=( eigenstates[[j,k]].fBase[Range[dim],r/λ]);
sign = Sign[ψ[0.5]]*Sign[psiTh[k,0.5]];
δψ = Max[Abs[ψ[#1]sign-psiTh[k,#1]]&/@rgrid];
δEth = Abs[eneAndLambda[[j,k,1]]-eneTh[[k]]];
{δEth,δψ},{k,state}],{j,1,Length[dimRange]}];

δEδψTh has 3 indexes and store the results in the form

δEδψTh[[dimension, level,1]] → δETh;
δEδψTh[[dimension, level,2]] → δψ = Max[|ψ[r]-ψTh[r]|];

We can now add to more columns to previous tables and have a direct comparison with known results

TableForm[
Table[TableForm[Table[Flatten[{NumberForm[eneAndLambda[[j,lev,2]],3],NumberForm[eneAndLambda[[j,lev,1]],16],
ScientificForm[#1,3]&/@ {δE[[j,lev]],δETemple[[j,lev]],δEδψTh[[j,lev,1]],δEδψTh[[j,lev,2]]}
}],{j,1,Length[dimRange]}],
TableHeadings-> {dimRange,{ "  λ","       E["<>ToString[lev]<>"]","  δE ","δETemple ","δEth", "  δψ"}},
TableSpacing-> {4,4}],{lev,state}],TableSpacing-> {10,5}]

Comments

The general features outlined above are confirmed, in particular δETemple appears to be a reasonable estimator of the error, δE is a bit pessimistic but this can be a good feature in more complicated cases.

Stability and precision of eigenvalues

The goal of this section is to show how we can overcome the limitations noted above. This is not a merely academic problem, as it can be a useful technique if one has to test new methods and their stability.

The first point to be stressed is that if we want high precision numbers in Mathematica  result, we must use high precision numbers also for the input parameters. In our case the external potential.  For the particular case we are studying V0 = 10 has infinite precision, but we write things in such a way that procedure can work for arbitrary V0.

V0=10;
V0 = SetPrecision[V0,100];
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) ;

We have written our functions, as laguerreRoots,  etc. in such a way that either the results are exact, as for laguerreMatrices, or their precision is controlled by the optional parameter prec which can be passed to them.

This simple procedure does not work for our function energy[] as it used Eigenvalues to solve a generalized eigenvalue problem, with matN ≠ 1. Mathematica routine in this case accepts only machine precision numbers, a precision between 15 and 16. We have then to invert matrix matN and bring it to left hand side of the eigenvalue equation. This can be a little dangerous because for large basis the inversion of the non sparse matrix matN can produce a rapid downgrade in the accuracy of the results, i.e. rounding errors can be quite big.

This is the new version for the computation of energy. The reader will notice that all functions are used with the additional argument prec.

computeEneHighPrecision[state_,λIn_,dim_,prec_:MachinePrecision]:= Module[{λ,D2,matN,matBaseLag,matV,invMatN,Ham,precEne},
λ=SetPrecision[λIn,prec];
{D2,matN}=N[laguerreMatrices[dim],prec];
{xgrid,weights} =laguerreRoots[dim+extraDim,prec];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
invMatN = Inverse[matN];
Ham=invMatN.(((- 1)/(2(λ^2))) D2+matV );
ene0=Sort[Eigenvalues[Ham]][[state]]];

It would be a useless waste of computer (and reader's) time to execute computations when downgrading due to rounding errors is greater than precision and also to go to higher dimensions when the target precision has been reached.  We do computations inside a "While" loop, which stops execution when one of these two conditions is met. The function starts by augmenting our precisionGoal by 10,  to give some room to rounding errors.

To save time we compute λ only for dimensions less than 40:  beyond that limit the λ is kept fixed.

The final error is the one estimated by δE at the last step

energyEstimate[state_,λInput_,precGoal_,dimMax_]:=
Module[{dim },
dim = 10;
shiftdim = 10;
prec = precGoal+10;
precEne = prec;
ene = computeEneHighPrecision[state,λGSdim30,dim,prec];
enesafe = ene;
deltaE = 1;
λ = λInput;
While[ 
dim ≤ dimMax -shiftdim&& precEne ≥ precGoal && deltaE ≥ 10^(-precGoal),
dim = dim+shiftdim;
If[dim< 40,
{D2,matN}=N[laguerreMatrices[dim]];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
min=FindMinimum[energy[λvar,{state,0},V,{xgrid,weights},matBaseLag],{λvar,0.01,2.}];
λ=SetPrecision[λvar/.min[[2]],prec];];
enenew = computeEneHighPrecision[state,λ,dim,prec];
deltaEold =  Abs[ene - enesafe];
precEneold = Precision[ene];
deltaE = Abs[ene - enenew];
precEne = Precision[enenew];
enesafe = ene;
ene = enenew;
risOld = {enesafe,dim-shiftdim,precEneold,deltaEold};
Print[TableForm[
{{"dim=  " <> ToString[dim],
NumberForm[enenew,Ceiling[precGoal]+1],
" prec.= "<>ToString[precEne],"δE:",ScientificForm[deltaE,4],
"δETh:",ScientificForm[ Abs[eneTh[[state]] - enenew],4]
}}]];
];
If[dim ≤ dimMax && deltaE < deltaEold,{ene,dim,precEne,deltaE},risOld]
]

We consider the ground state and fix a precision of 20 as our goal. The results are printed with 20 decimal digits in this case, in general with prec+1 digits for the complete number.

state=1;
precGoal=20;
maxDim = 120;

energyEstimate[state,0.1,precGoal,maxDim]//Timing

In this case dim = 70 was sufficient to get a 20 decimal digits precision result.

The same procedure can be used for the excited state:

energyEstimate[2,0.1,precGoal,maxDim]//Timing

Comments on  gaussian integration and the dimension of the basis

We finish this notebook with some words on the approximation used and point out some problems which can arise in more complicated and realistic cases.

First of all let us show that effectively the "fake" minima on E[λ] eigenvalues are produced by gaussian approximation. The simplest thing would be a comparison with more accurate computation of integrals, i.e., with NIntegrate. This is really time-consuming for large enough dimensions.  We just write here a couple of functions which can be used for this purpose:  one has to replace the function energy with energyNint written below.

getMatVNint[λ_?NumberQ, V_, dim_] := Module[{mat1, mat2},
  mat2 = ConstantArray[0, {dim, dim}];
  mat1 = mat2;
   Table[mat2[[i, j]] = NIntegrate[V[λ x] fBase[i, x] fBase[j, x], {x, 0, ∞}], {i, 1, dim}, {j, i + 1, dim}];
  Table[mat1[[i, i]] = NIntegrate[V[λ x] fBase[i, x] fBase[i, x], {x, 0, ∞}], {i, 1, dim}];
  mat1 + mat2 + Transpose[mat2]]

energyNint[λ_?NumberQ,{state_,L_},V_,dim_]:= Module[{matV},
matV=getMatVNint[λ,V,dim];
Sort[Re[Eigenvalues[{((- 1)/(2(λ^2))) D2+matV + (1/(2(λ^2)))L(L+1)IdentityMatrix[dim],matN}]]][[state]]]

A run would require several minutes and we choose instead to use an increasing number of the grid points.

Here we reproduce the "critical" region of the graph for the second eigenvalue varying the dimension of the grid. Our parameter extraDim being 4 the grids used for integrals are of dimension dim+4 (the usual one) dim+8 and 2×dim

state=2;L=0;

λPlot = Join[Range[0.05,0.9,0.05],Range[1,3,0.1]];
λPlot=Range[0.5,2.3,0.05];
dimRange=Range[60,60,10];

ene1=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
energy[#1,{state,L},V,{xgrid,weights},matBaseLag]&/@λPlot,
{dim,dimRange}];

ene2=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+2extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
energy[#1,{state,L},V,{xgrid,weights},matBaseLag]&/@λPlot,
{dim,dimRange}];

ene3=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+dim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
energy[#1,{state,L},V,{xgrid,weights},matBaseLag]&/@λPlot,
{dim,dimRange}];

ListPlot[{Transpose[{λPlot,ene1[[1,All]]}],
Transpose[{λPlot,ene2[[1,All]]}],
Transpose[{λPlot,ene3[[1,All]]}]},PlotStyle-> {Dashing[0.02],Automatic,Dashing[0.01]}
,Joined-> True,PlotLabel-> Row[{"eig: ",2,";  dim= ",dimRange[[1]]}]]

We see that for the highest number of gaussian points (120 in this case) the function is convex and no ambiguity in the identification of stationary range for λ appears.

In problems like this one the situation does not cause real troubles: the stable zone in λ is readily identified and to use very large grids would be just a waste of  time. Things are more complicated when the model is nearly critical, i.e. if we have an eigenvalue very small. In this case bound states can have a large characteristic length and convergence is in general more difficult.

From analytical results it is easy to see that (Subscript[V, 0])=6 is a critical value for this potential, for lower values the excited state is absent. Let us look at the variation of energy as a function of λ and for different dimensions

V1=6.02;

eneVal=-(1/8)((((Sqrt[1+ 8 V1]) - (2 3 +1)))^2)

V[r_]:= - (V1/((Cosh[ r])^2));

λPlot = Join[Range[0.05,0.9,0.05],Range[1,3,0.1]];
λPlot=Range[0.5,3.5,0.05];
dimRange=Range[40,80,20];

Here is the usual computation

ene1=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim+extraDim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
energy[#1,{state,L},V,{xgrid,weights},matBaseLag]&/@λPlot,
{dim,dimRange}];

ListPlot[Table[Transpose[{λPlot,ene1[[k,All]]}],{k,1,Length[dimRange]}]
,Joined-> True,PlotLabel-> Row[{"eig: ",2}],
PlotRange-> {{0,2.5},{2eneVal,0}},PlotStyle-> {Red,Blue,Black},
Epilog-> {Dashing[0.02],Line[{{First[λPlot],eneVal},{Last[λPlot],eneVal}}]}]

We see that in this case with a low number of dimensions it is almost impossible to detect a stable range for λ.  The relative minima can be in fact an artifact of the approximation. In similar case a computation with a large basis is mandatory, low dimensional results can be completely misleading.

Here is the computation with an higher number of gaussian points, this computation is a bit slow.

ene2=Table[
{D2,matN}=N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[2 dim];
matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}];
energy[#1,{state,L},V,{xgrid,weights},matBaseLag]&/@λPlot,
{dim,dimRange}];

ListPlot[Table[Transpose[{λPlot,ene2[[k,All]]}],{k,1,Length[dimRange]}]
,Joined-> True,PlotLabel-> Row[{"eig: ",2}],
PlotRange-> {{0,2.5},{2eneVal,0}},PlotStyle-> {Red,Blue,Black},
Epilog-> {Dashing[0.02],Line[{{First[λPlot],eneVal},{Last[λPlot],eneVal}}]}]

The situation is more stable and it is easier to identify a stable zone, but also in this case a large basis is required.

The message is that in the usual conditions the necessity of a large basis is not due to inadequacy of gaussian integration procedure but to a large length scale in the problem. Use of large sets of gaussian points can help in the identification of stable ranges for λ, but does not help avoiding the physical problem of large scales.