Laguerre basis

Quantum Mechanics
 K. Konishi, G. Paffuti

Use a basis of Laguerre polynomials to solve radial problems with the variational method.

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

Laguerre basis

We have run into Laguerre polynomials several times: in the study of radial eigenstates of Hydrogen atom, in an example of the  Gauss formulas etc. The main point of interest here is the completeness property of the set in the interval [0,∞]

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

The normalization used in this book is the one used in Mathematica,

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

These properties make the set as a natural candidate for the study of radial problems in Quantum Mechanics, as 0 ≤ r ≤ ∞.

While the set () is perfect for computing integrals, and for many other uses, it is not best suited for radial Schrödinger equation. It is easily seen that the momentum operator is not represented by an adjoint matrix in the basis (). For example ⟨1| p| 1⟩ ≠ 0 and purely imaginary: this is a pedestrian way for seeing that the radial momentum is not a self-adjoint operator. Related to this problem, and of more interest to us, is the fact that in this set it is not easy to implement boundary conditions. The typical radial Schrödinger equation for the reduced function f[r] = r R[r] has he form

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

The boundary condition is not satisfied by functions (), and this is a quite big problem as it is not true that an arbitrary linear combination of those functions can be a trial function for a variational approach.

Fortunately the cure to this problem is quite easy. We just multiply the set by x. Our set of basis functions which satisfy zero boundary conditions is

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

To conform with the usual conventions we start counting vectors from 1.  This is why we shifted the index of the Laguerre polynomials. The price we have to pay in using () is the loss of orthogonality, but this is not a big issue in a variational approach. We will always denote by N the matrix*

⟨(Subscript[f, i])|(Subscript[f, j])⟩ = (Subscript[N, ij]).

*) We will use the symbol matN in Mathematica  instructions to avoid conflicts with the Mathematica function N (which means the numerical value of).  See for instance Eq.(1.8) below.

The Laguerre polynomials satisfy several identities, those of interest to us are the following recursive relations

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

from the first one it easily follows that :

(d/dx)(Exp[-x/2] x L[n,x]) = Exp[-x/2](((n+1)/2)L[n+1,x] + (1/2)L[n,x] - (n/2)L[n-1,x])

These relations and the orthogonality conditions () allow us to write almost by inspection, or with an easy computation, the matrix elements of several operators. We list only non vanishing entries:

matN[k,n] = (Piecewise[{{6(n^2)+6n+2;,  if k=n }, {-4(n^2);,        if k=n-1}}])       ;     matN[n,k] = matN[k,n];

D1[k,n] = (Piecewise[{{-(n^2);,         if k=n-1 }, {n(n-1)/2;,        if k=n-2}}])       ;     D1[n,k] = D1[k,n]

D2[k,n] = (Piecewise[{{-((n^2)+n+1)/2;,         if k=n }, {n(n-1)/4;,        if k=n-2}}])   ;    D2[n,k] = D2[k,n]

D1 and D2 are the first and second derivative respectively.

Before continuing let us note that, denoting the "bad" basis () by (Subscript[b, k]) we have immediately, for any operator

⟨(Subscript[f, i])|A|(Subscript[f, j])⟩ = ⟨(Subscript[b, i])|x A x|(Subscript[b, j])⟩

This allows a very easy computation of power of position operators using the simple recursion relation (). For example the operator 1/(x^2)is just the identity matrix

⟨(Subscript[f, i])|(1/(x^2))|(Subscript[f, j])⟩ = ⟨(Subscript[b, i])|(Subscript[b, j])⟩= (Subscript[δ, ij])

While 1/x is given by

⟨(Subscript[f, i])|(1/x)|(Subscript[f, j])⟩ = ⟨(Subscript[b, i])|x|(Subscript[b, j])⟩

x[k,n] = (Piecewise[{{2n+1;,  if k=n }, {-n;,        if k=n-1}}])       ;     x[n,k] = x[k,n];

In the same way, positive powers can be computed, but it is even simpler to note that in b basis x operator has matrix elements only for elements which differ at most by one place, i.e. it is a tridiagonal matrix. If we work with matrix n×n it is sufficient to enlarge the basis to (n+δ)×(n+δ), make matrix products and keep the resulting  n×n submatrix. It is easily seen that δ is the power of x added by 2, as one must take into account the first and last x coming from passage form b basis to f basis. This feature will be used below to write a function which produces matrices for arbitrary positive powers of the radial variable.

Comment on non orthogonal bases

One must take some care in using non-orthogonal basis in translating operators products in matrix products. The main point is that the identity of a complete set of eigenstates is a bit more complicated than usual, in effect the identity operator is given by

𝕀 = ∑ |(Subscript[f, i])⟩ (Subscript[N, ij]^(-1)) ⟨(Subscript[f, j])|

This is easily seen by applying this operator to basis vectors

|(Subscript[f, k])⟩ = 𝕀 |(Subscript[f, k])⟩ =  ∑ |(Subscript[f, i])⟩ (Subscript[N, ij]^(-1)) ⟨(Subscript[f, j])|(Subscript[f, k])⟩ =  ∑ |(Subscript[f, i])⟩(Subscript[N, ij]^(-1)) (Subscript[N, jk]) = |(Subscript[f, k])⟩.

This means that an operator relation of the form C = A B translate, as matrices, in

matC = matA.(N^(-1)).matB

This is also easily understood if one works in the b basis. N matrix is the matrix of (x^2)in this basis and we have, writing explicitly the x factors to convert from the b to f basis

x C x = (x A x)(1/(x^2))(x B x).

Variational calculation

Our basis will be enriched by a parameter λ, i.e.

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

This parameter is reabsorbed in the matrix elements by a scale transformation  x → λ x. Homogenous quantities, like derivatives or homogenous potentials, change by a power of λ. The variational principle is written as

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

which, in the rescaled variables, produce the matrix equation  (we will always use natural units):

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

matV is the matrix computed with the rescaled potential V. For each ℰ[λ] a minimization with respect to λ is understood.

N.B. We choose to write x/λ as argument for the excited states, which have typically large characteristic lengths, corresponding to large values of λ. After the computation of eigenvalues and eigenstates, if one want to recover the usual wave functions one must make the inverse transformation x → λ x. 
In general we shall try to use the symbol r for the usual radial variable, while x will denote rescaled variables,  r = λ x. In particular the potential function under rescaling will be written as (Subscript[V, λ])[x] ≡ V[λ x].

Norms

Let ψ be a state definite by its spectral coefficients, computed for instance by solving an eigenvalue problem. In original non-rescaled coordinates

ψ = (∑_k) (Subscript[c, k]) (Subscript[f, k]) [ r/λ]

Its square norm is then given by (we take for simplicity all coefficients to be real):

((‖ψ‖)^2)=  (∑_(k,j)) (Subscript[c, k]) (Subscript[c, j]) ∫𝕕r (Subscript[f, k]) [ (r/λ)](Subscript[f, j]) [ (r/λ)] = λ (∑_(k,j)) (Subscript[c, k]) (Subscript[c, j]) ∫𝕕x (Subscript[f, k]) [ x](Subscript[f, j]) [ x]= λ (∑_(k,j)) (Subscript[c, k]) (Subscript[N, kj]) (Subscript[c, j])

If we want normalized states we have to divide by the square roots of this factor.

Functions

Matrices

Translating in Mathematica  language what has been said above we can immediately write the function laguerreMatrices below. In input it is necessary to specify the dimension of the basis and as an optional argument a list with power of wanted radial variables.

In output there is a list on matrices. The first two are D2 and the normalization matrix N, then the requested powers follow.

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]]

Wave functions

The basis is defined by

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

It is useful to have at our disposal a tool which computes quickly the wave functions, for plotting or for 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 the 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, possibly, the reader has to call 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 vector of coefficient by a sequence of basis functions. This is done simply with:

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

The reader has to provide, if necessary, the normalization of the vector and the possible scale factor λ. The trivial check of consistency for the function laguerreP  is left as an exercise.

Checks

It is useful to have some checks for the function. The basis is defined by

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

All integrals are of the form  (x^a) Exp[-x] and can be computed by defining a rule (it is supposed that the variable is x):

myInt[funz_]:= Block[{a},Together[Expand[funz]//.{x^a_. Exp[-x]-> a!, Exp[-x]-> 1}]];

Let us define a small dimension to make the checks

dim = 5;

Here we compute D2, the norm matrix, and potential 1/r, r, (r^2).

{D2,matN,Rm1,R1,R2}=laguerreMatrices[dim,{-1,1,2}];

Here we compute the same matrices by performing explicit integrations :

checkD2 = Table[myInt[fBase[m,x]* (∂_(x,x))fBase[n,x]],{n,1,dim},{m,1,dim}];
checkmatN =  Table[myInt[fBase[m,x]* fBase[n,x]],{n,1,dim},{m,1,dim}];
checkRm1 = Table[myInt[fBase[m,x]1/x  fBase[n,x]],{n,1,dim},{m,1,dim}];
checkR1 = Table[myInt[fBase[m,x] x  fBase[n,x]],{n,1,dim},{m,1,dim}];
checkR2 = Table[myInt[fBase[m,x] (x^2)  fBase[n,x]],{n,1,dim},{m,1,dim}];

and compare the two expressions

{D2== checkD2,matN== checkmatN,Rm1== checkRm1,R1== checkR1,R2== checkR2}

We can also check an approximate form of the relation () (the approximation is due to the missed matrix elements at the boundary of the matrices)

{MatrixForm@(R1.Inverse[matN].R1),MatrixForm@R2}

It is seen that the two matrices coincide in a subspace of dimension dim - 1.

Example

Hydrogen atom: energy levels

We take the hydrogen atom as a simple model to check the method. In this case, the potential is homogeneous and the computation is simplified. We take for the moment s-states, then V[r] = - 1/r, but we shall compute everything for a generic L.

One starts fixing the dimension

dim = 50;

computes the matrices

{D2,matN,Rm1} = N@laguerreMatrices[dim,{-1}];
Rm2 = IdentityMatrix[dim];

and writes a function which computes the eigenvalues for a fixed λ. As our matrices will be rather small we use the full version of Eigenvalues to compute the results.  It will be shown later how this procedure can be modified.

energy[state_,L_?IntegerQ,λ_?NumberQ]:=Sort[Re[Eigenvalues[{((- 1)/(2(λ^2))) D2-(1/λ)  Rm1 + (1/(2(λ^2)))L(L+1)Rm2,matN}]]][[state]];

Here we choose the state, in particular its angular momentum

L=0; state = 1;

We have chosen the first s - state, i.e. 1s.

We choose an arbitrary starting point for our parameter, (Subscript[λ, start])

λstart = 0.7;

as a check

energy[state,L,λstart]

It is already the correct result, but we shall write anyway a routine to minimize with respect to λ, first let us have a look at λ dependence of the eigenvalues

λplot=Range[0.5,50,0.5];

Show[
Table[ListPlot[Transpose[{λplot,Evaluate[(energy[i,L,#1]&/@λplot)]}],
AxesOrigin-> {0,-1},
PlotRange-> {All,{-1,0}},Joined-> True,
AxesLabel-> {"λ","eig"}],{i,1,5}]]

It is seen that for large λ, the ground state energy is not well reproduced.  A similar conclusion holds for excited levels.

((Range[15])^2)energy[Range[15],L,λstart]

As an empirical rule the value λ = dim/12 seems to give a reasonable stability for the first few states, up to dim/3.  Beyond this limit the eigenvalues are little related to the real coulombic problem. As dim → ∞ we expect that every fixed set of levels reach a stable value.

One needs a procedure which automatically selects the best value for λ. The simplest procedure is to use the Mathematica  function FindMinimum:

min=FindMinimum[energy[state,L,λvar],{λvar,0.01,100}]

N.B. With our definition of λ higher states tend to have very high λ, one must check the extrema of the interval in FindMinimum. It is not advisable to use the form of FindMinimum with a single starting point: this form uses derivatives with respect to the parameters. The derivatives are usually small in a region of stability for λ and the function give a warning message.

The energy and λ can be extracted from the variable min defined above:

ene0=min[[1]];
λ = λvar/.min[[2]]

Wave function

Let us give now an example in which also the wave function is computed. Let us consider the third state with L=1 of hydrogen atom, i.e. the state 4p. The "state" is actually the radial quantum number, the principal quantum number  of the hydrogen atom,  n = (Subscript[n, r])+L = state+L.

state =3; L=1;

dim = 50; 
{D2,matN,Rm1} = N@laguerreMatrices[dim,{-1}];
Rm2 = IdentityMatrix[dim];

ψ[x] = ∑ (Subscript[c, k])(Subscript[f, k])[x/λ];

min=FindMinimum[energy[state,L,λvar],{λvar,0.3,100}];
{ene,λ} = {min[[1]],λvar/.min[[2]]}

The value is the correct one :

eneth = - 0.5/(((state+L))^2)

To compute the coefficients c  one makes a call to Eigensystem with the value of λ just computed.

H =- (1/(2(λ^2)))D2 - (1/λ)Rm1 +((L(L+1))/(2(λ^2)))Rm2;

{eigH,eigV}=Eigensystem[{H,matN}];
order = Ordering[eigH];
eigH=eigH[[order]];
eigV = eigV[[order]];

eigH[[state]]

coeff = Re[eigV[[state]]];

To compute the eigenfunction one uses the function getPsi defined in the first section.

To normalize the coefficients we use relation ():

norm2 = λ coeff.matN.coeff;
coeff = coeff/(Sqrt[norm2]);

The wave function is

ψ[r_?NumberQ]:= getPsi[coeff,r/λ];

An alternative formula is :

ψ[r_, λ_] := Sum[coeff[[k]] fBase[k, r/λ], {k, 1, dim}]

Let us verify the normalization :

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

The known reduced functions for hydrogen atom are :

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)] ;

Plot[{Rreduced[r,state+L,L],ψ[r]},{r,0,60},PerformanceGoal-> "Speed"]

They are indistinguishable on this scale, so let us plot their difference (the scale is (10^(-14))):

Plot[Evaluate/@((10^14)Abs[Rreduced[r,state+L,L]-ψ[r]]),{r,0,60},AxesLabel-> {"r","(10^(-14))"},PlotRange-> All]

Generic potential

To work with a generic potential one of the possibility is to compute all integrals using Gauss' integration formulas. From notebook[*] we copy the function which generates abscissas and weights for:

(∫_0^∞)Exp[-r]F[r] 𝕕r ≃ (∑_k) (Subscript[w, k])F[(Subscript[x, k])]

We choose to make explicit the exponential factor, the same appearing in the definition of our basis functions.

laguerreRoots[n_]:=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]];
w=(radEig/((((n+1))^2)((LaguerreL[n+1,radEig])^2)));
{radEig,w}]

We compute usual matrices

dim = 50;

{D2,matN,Rm1} = N@laguerreMatrices[dim,{-1}];
Rm2 = IdentityMatrix[dim];

and grid points

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

The values of basis functions on grid points  must be repeatedly used, so it is better to compute them once for all. We omit the exponential term, which will be taken into account by the weights of Gauss formula.

matBaseLag=Transpose[Transpose[laguerreP[dim,xgrid]]*xgrid];

An alternative formula is

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

With this matrix the matrix element of any potential can be computed by using Gauss' formulas :

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]];

We define a potential for our problem

V[x_]:= - 1/x

and define which state we want to study :

state=3; L=1;

One must write a function which computes energy, but, in contast to the simple Coulombic case, the form of the Hamiltonian can change with λ, so one computes inside the function the matrix matV, for the potential.  The matrix elements are computed with Gauss' formula:

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

min=FindMinimum[energyGenericV[dim,λvar,V,{xgrid,weights},matBaseLag],{λvar,0.5,20}];

{{ene,λ} ={ min[[1]],λvar/.min[[2]]},"eneth=",-0.5/(((state+L))^2)}

We write a function which need as input only the dimension we want to use, the potential and the state, and produce eigenfunctions and eigenvalues. An additional input parameter is the interval in which to look for λ, the default is [0.01,5] but it is much better to have an idea of the interval by plotting eigenvalues as a function of λ. This can be done with the function getEnergy[], given below. Inside the main routine the function energy is redefined to avoid a recomputation of matrix elements once λ has been fixed.

The functions defined outside this one are:

As a by-product the function defines a wave function for the level, the name of the function is psiGenericV[r]. It depends on the original, non rescaled, variable r.

N.B. The right value of the interval in which FindMinimum look for the optimal λ value can affect the results on the wave function, the reader in each case must check the stability.

getEnergy[dim_,state_,L_,λ_,V_]:= Module[{matV1,matV2,Vk,xgrid,weights,matBaseLag,matV},
{D2,matN} = N@laguerreMatrices[dim];
{xgrid,weights} =laguerreRoots[dim];
(* matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}]; *)
matBaseLag=Transpose[Transpose[laguerreP[dim,xgrid]]*xgrid];
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
Sort[Re[Eigenvalues[{((- 1)/(2(λ^2))) D2+matV + (1/(2(λ^2)))L(L+1)IdentityMatrix[dim],matN}]]][[state]]];

Here we show a plot of the eigenvalues as functions of λ. While in this simple case (hydrogen atom) the choice of the interval search for λ is quite stable, in different potentials the situation can be worse and it is better to have an idea of the optimal interval before doing any heavy computational work.

λRange= Join[Range[0.025,1,0.025],Range[1.25,10,0.25]];

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

enePlot=getEnergy[20,state,L,#1,V]&/@λRange;

ListPlot[Table[{λRange[[k]],enePlot[[k,p]]},{p,1,Length[state]},{k,1,Length[λRange]}],
AxesLabel-> {"λ","E"},
PlotLabel-> "dim = "<>ToString[dim],Joined-> True]

Here is the general function

energyAndPsi[diminput_,stateinput_,Linput_,λInterval_:{0.1,5}]:=
 Module[{D2,matN,Rm1,Rm2,H,eigH,eigV,order},
Block[{dim,state,L,matBaseLag,xgrid,weights,min,matV,λI},
λI = Flatten[λInterval];
{dim,state,L} = {diminput,stateinput,Linput};
{D2,matN,Rm1} = N@laguerreMatrices[dim,{-1}];
Rm2 = IdentityMatrix[dim];
{xgrid,weights} =laguerreRoots[dim+4];
(* matBaseLag= Table[xgrid*LaguerreL[k-1,xgrid],{k,1,dim}]; *)
matBaseLag=Transpose[Transpose[laguerreP[dim,xgrid]]*xgrid];

energyLocal[dim_,λ_?NumberQ,V_,{xgrid_,weights_}]:= Module[{matV1,matV2,Vk},
matV1 = 0*IdentityMatrix[dim];
(* the values of the potential are computed  on the rescaled grid *)
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}];
matV=DiagonalMatrix[matV2] + matV1 + Transpose[matV1];
Sort[Re[Eigenvalues[{((- 1)/(2(λ^2))) D2+matV + (1/(2(λ^2)))L(L+1)Rm2,matN}]]][[state]]];
If[Length[λI]== 1,
ene = energyLocal[dim,λI[[1]],V,{xgrid,weights}];
λ = λI[[1]];,
min=FindMinimum[energyLocal[dim,λvar,V,{xgrid,weights}],{λvar,λI[[1]],λI[[2]]}];
{ene,λ} ={ min[[1]],λvar/.min[[2]]};]
H =- (1/(2(λ^2)))D2+ matV +((L(L+1))/(2(λ^2)))Rm2;
{eigH,eigV}=Eigensystem[{H,matN}];
order = Ordering[eigH];
eigH=eigH[[order]];
eigV = eigV[[order]];
coeffgenericV = Re[eigV[[state]]];
norm1 =Sqrt[ λ coeffgenericV.matN.coeffgenericV];
psiGenericV[r_]:=(1/norm1) Sum[coeffgenericV[[k]](r/λ)Exp[-r/(2λ)]LaguerreL[k-1,(r/λ)],{k,1,diminput}];
ene]];

Check

As a check we test again the Coulomb potential. We use default values for λInterval.

V[x_]:= -1/x;

ene = energyAndPsi[50,3,1]

Here is a comparison with respect exact eigenfunction :

sign=Sign[psiGenericV[0.01]];
Plot[Evaluate/@{Rreduced[r,4,1],sign*psiGenericV[r]},{r,0,60}]

FindMaximum[Abs[Rreduced[r,4,1]-sign*psiGenericV[r]],{r,0,60}]

Orthogonal basis

The method presented in the previous section can be considered a standard method to solve radial problems. We can bring some change in it by diagonalizing the set of the basis function. In fact there is no computation to be done, it is known that the generalized Laguerre polynomials L[n,2,x] are orthogonal with respect the measure (x^2)Exp[-x], i.e. the one generated by the prefactor x in the old form of the basis functions. These polynomials define then the orthonormal basis

f[n,x]=(1/(Sqrt[n(n+1)])) Exp[-x/2] x  LaguerreL[n-1,2,x];    n =  1 , 2

The normalization is the correct one to get norm one functions. This basis certainly simplifies the algebra, the corresponding operators - matrices in more direct. There is no need to solve a generalized eigenvalue problem to find the energy levels, i.e. the matrix matN used above is the identity. The price to pay is that matrices for derivatives are no longer tridiagonal, and sparse. Some calculations show that for second derivatives:

D2[i,j] = (Piecewise[{{- (1/(Sqrt[i(i+1)j(j+1)]))((i(i+1)(4i-1))/12);,         if i=j }, { - (1/(Sqrt[i(i+1)j(j+1)]))((i(i+1)(2i+1))/6);,        if i<j}}])   ;    D2[i,j] = D2[j,i]; i=1,2,...

This matrix can be generated with the following function

laguerreOrthoD2[dim_]:= Module[{mat2,D2},
mat2 = 0*IdentityMatrix[dim]; 
Table[mat2[[i,j]]=- (1/(Sqrt[i(i+1)j(j+1)]))((i(i+1)(2i+1))/6),{i,1,dim-1},{j,i+1,dim}];
D2 = mat2+Transpose[mat2];
Table[D2[[i,i]]= - (1/(Sqrt[i(i+1)i(i+1)]))((i(i+1)(4i-1))/12),{i,1,dim}];
D2];

The basis vectors are defined by

fBase[n_,x_]:=(1/(Sqrt[n(n+1)])) Exp[-x/2] x  LaguerreL[n-1,2,x];

Check

Let us give a simple check of the above computation

dim = 5;
D2= laguerreOrthoD2[dim];
Table[myInt[fBase[i,x]*(∂_(x,x)) fBase[j,x]],{i,1,dim},{j,1,dim}] == D2

Table[myInt[fBase[i,x]* fBase[j,x]],{i,1,dim},{j,1,dim}] == IdentityMatrix[dim]

Variational computation

Almost nothing must be changed from the previous section, except the base definitions for Laguerre polynomials on the grid, and the absence of matN in the eigenvalue problem:

dim = 20;

D2= N@laguerreOrthoD2[dim];

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

L=1; state=3;

The centrifugal potential now must be included in V, as we have not defined a matrix for it. We take as usual the Coulomb potential for a check

V[x_]:= - (1/x) + L(((L+1))/(2 (x^2)))

The base for polynomial components of the basis changes, involving LaguerreL[n, 2, x]

matBaseLag= Table[(1/(Sqrt[k(k+1)]))xgrid*LaguerreL[k-1,2,xgrid],{k,1,dim}];

The function getMatV is unchanged and gives us the matrix corresponding to our potential, for a fixed λ. For energy we drop matN from the function:

energyGenericV[dim_,λ_?NumberQ,V_,{xgrid_,weights_},matBaseLag_]:= Module[{matV},
matV=getMatV[λ,V,{xgrid,weights},matBaseLag];
Sort[Re[Eigenvalues[((- 1)/(2(λ^2))) D2+matV ]]][[state]]]

min=FindMinimum[energyGenericV[dim,λvar,V,{xgrid,weights},matBaseLag],{λvar,0.5,20}];

{ene,λ} ={ min[[1]],λvar/.min[[2]]}

One can study eigenfunctions, etc.