Hermite basis

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the variational method using the Hermite polynomials.

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

Introduction

In this notebook we introduce a basis for Hilbert space built on Hermite polynomials,  i.e. the eigenstates of the Harmonic oscillator. This basis is best suited for the study of one-dimensional problems in the interval [-∞,+∞]. We will present both the full variational method (or spectral method)  and its approximate version, based on the computation of the matrix elements with Gaussian integration formulas. This last version is more easily applied to generic potentials.

Here we present only general methods and their simplest application. More specialized topics like error estimations, systematic analysis of scale dependence, special systems etc. will be covered in other notebooks

In all notebooks we will use the natural units: ℏ = m = 1.

Variational (spectral) method

The eigenfunctions of a harmonic oscillator (with ω = 1) are a complete orthonormal set for the Hilbert space (𝕃^2):

(1/(π^(1/4))) (1/(Sqrt[(2^n)n!]))(Subscript[H, n])[x]Exp[-x^2/2];  n = 0,1,...

To conform to usual convention in Mathematica , where vector indexes start from 1, the basis will be often written as

(Subscript[f, k])[x] = (1/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))(Subscript[H, k-1])[x]Exp[-x^2/2];   k = 1,2,...

The basis is orthonormal:

(∫_(-∞)^(+∞))(Subscript[f, k])[x] (Subscript[f, j])[x]𝕕x = (Subscript[δ, kj])

A truncation of this set, for k ≤ dim will be taken as a basis set for variational calculations. The functions will be usually equipped with a free scale parameter, α, and our trial wave functions will be looked for in the form

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

The variation principle states that

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

Inserting expansion () it is easily seen that stationarity with respect c's imply the usual truncated form of matrix Schrödinger equation:

(Subscript[H, ij])(Subscript[c, j]) = E (Subscript[c, i])

The eigenvalues E will be in general α dependent, and one must minimize with respect to this parameter.

We note that the orthonormality of the basis implies that in the right hand side of equation () only identity matrix appears.

The scale parameter α can be conveniently reabsorbed in a change of variables:

r = x/α;        x = α r.

Correspondingly matrix elements of () will be fixed by the Hamiltonian (in x -variables)

(Subscript[H, x]) = (1/2)(α^2)(p^2) + V[x/α]

Explicitly

(Subscript[p, kj]^2) = (∫_(-∞)^(+∞))((∂_x)(Subscript[f, k])[x])((∂_x) (Subscript[f, j])[x])𝕕x;    (Subscript[V, kj])= (∫_(-∞)^(+∞))(Subscript[f, k])[x]V[x/α] (Subscript[f, j])[x]𝕕x.

Clearly for homogeneous potentials the α factor in matrix elements of V can be factorized. The overall 1/α factor of the jacobian is not included in these formulas as it is common to all matrices and cancel in the equation ().

Scale and transformation factors

N.B. Please note that this choice is the opposite of that made for Laguerre basis, where λ=1/α was used. This convention turn out to be simpler in this case. When there is a mixing between original non rescaled coordinates and scaled coordinates we will try to adhere to the convention of calling r the original coordinates and x the rescaled ones. So, once an eigenvector of coefficient (Subscript[c, k])will be found, the eigenfunction will be written with a code like:

ψ[r_] := Sum[c[[k]] f[k, α r],{k,1,dim}]

To pass from matrix elements in x space to matrix elements in r coordinates on must keep track of the jacobian. So for scalar products, norms etc we have the table ( D stands for derivative or momentum operator and we assume all coefficients real for simplicity):

ψ[r] = (∑_i)(Subscript[c, i])(Subscript[f, i])[α r];     φ[r]= (∑_i)(Subscript[d, i])(Subscript[f, i])[α r];

⟨φ|ψ⟩ = (∑_ij)(Subscript[d, i]) (Subscript[c, j])∫𝕕r (Subscript[f, i])[α r](Subscript[f, j])[α r]= (1/α)(∑_i) (Subscript[d, i])(Subscript[c, i]);

⟨φ|V[r]|ψ⟩= (∑_ij)(Subscript[d, i]) (Subscript[c, j])∫𝕕r (Subscript[f, i])[α r]V[r](Subscript[f, j])[α r]= 
                    = (1/α)(∑_ij)(Subscript[d, i]) (Subscript[c, j])∫𝕕x (Subscript[f, i])[x]V[x/α](Subscript[f, j])[x]= (1/α)(∑_ij)(Subscript[d, i]) (Subscript[c, j]) (Subscript[V, ij])

⟨φ|(D^((L)))|ψ⟩= (∑_ij)(Subscript[d, i]) (Subscript[c, j])∫𝕕r (Subscript[f, i])[α r](∂_r^((L)))(Subscript[f, j])[α r]= 
                     = ((α^L)/α)(∑_ij)(Subscript[d, i]) (Subscript[c, j])∫𝕕x (Subscript[f, i])[x](∂_x^((L)))(Subscript[f, j])[x]= ((α^L)/α)(∑_ij)(Subscript[d, i]) (Subscript[c, j]) (Subscript[D, ij]^((L)))

N.B. On normalized states the overall 1/α factor cancels and matrix elements can be computed directly in x-representation, with α retained in the definition of V[x/α] or (α (∂_x)) obviously:

((⟨φ|ψ⟩)/(‖φ‖‖ψ‖)) = (∑_i) (Subscript[d, i])(Subscript[c, i]);       ((⟨φ|V[r]|ψ⟩)/(‖φ‖‖ψ‖))=(∑_ij)(Subscript[d, i]) (Subscript[c, j]) (Subscript[V, ij]);
((⟨φ|(Subscript[D, r]^((L)))|ψ⟩)/(‖φ‖‖ψ‖))=  (∑_ij)(Subscript[d, i]) (Subscript[c, j]) (Subscript[D, ij]^((L)))(α^L)

Matrix elements

Matrix elements for kinetic term and for simple powers of x - coordinate are well known from oscillator harmonic theory and have been already used in other notebooks. In the following formulas we will put ℏ = m = ω =1. These formulas are written in standard quantum mechanics convention, with n=0,1,... . All elements omitted are zero.

<n| x | n-1> = <n-1| x |n> = (Sqrt[(ℏ/(mω))])(Sqrt[(n/2)]);
<n-1|p|n> = -i (Sqrt[m ω ℏ])(Sqrt[(n/2)])  = (SuperStar[(<n |p|n-1>)])

<n| (x^2) | n> =  (ℏ/(mω))((2n+1)/2);
<n| (x^2) | n+2> = <n+2| (x^2) |n> = (ℏ/(mω))(Sqrt[(((n+1)(n+2))/4)]);
<n| (x^4) | n> =  ((((ℏ/(mω))))^2)(3/4)(2(n^2)+2n+1);
<n| (x^4) | n+2> = <n+2| (x^4) |n> = ((((ℏ/(mω))))^2)((2n+3)/2)(Sqrt[(n+1)(n+2)]);
<n| (x^4) | n+4> = <n+4| (x^4) |n> = ((((ℏ/(mω))))^2)(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)]);
<n| (p^2) | n> =  ℏ mω((2n+1)/2);
<n| (p^2) | n+2> = <n+2| (p^2) |n> = ℏ mω(Sqrt[(((n+1)(n+2))/4)]);

The main feature of these matrices is their band structure : only few diagonal are different from zero, and this allow an efficient use of sparse arrays for large matrices, at the same time numerical algorithms for solving eigenvalue problems are more stable.

Parity and higher powers of x

If problem is parity invariant a  simplification can be used, only even or odd states can be used to construct the matrices. It is easily checked that for k = 1,2,3 .. even and odd states refers to levels:

n = 2k-2 + parity;       parity  =0:  even ;   parity = 1:    odd

For instance for parity = 0, as k grows we have levels n=0,2,4 ... while for parity = 1,  n = 1,3,5, ...  . This allow an easy identification of even and odd states.

For polynomial interactions these matrix elements are sufficient. If higher power of x are required one can compute matrices for higher dimensions and using the band structure extract relevant matrix elements. For example if (x^6)operator is required one can use dim+3 matrices for fixed parity or dim+6 matrices in the case full matrices, for both parity states, are required.

Wave functions computation

Once we have computed spectral coefficients (Subscript[c, k]) we could use the analytical expression () for wave functions. This turn out not to be very practical from the point of view of computer time. A command like Plot[ψ[r] ... ] would take some time to execute,mainly because a lot of Hermite polynomials have to be computed, and for each of them we use a different Mathematica instruction. It is better to use ListPlot on a predefined grid, in such a way each polynomial is compute on the whole grid, using vector structure of Mathematica functions.

It is even better to avoid the use of the functions HermiteH. This is easily done using recursion relation for Hermite polynomials

(Subscript[H, n+1])[x] = 2 x (Subscript[H, n])[x] - 2 n (Subscript[H, n-1])[x];  (Subscript[H, 0])[x]=1;  (Subscript[H, 1])[x] = 2x;

In this way simple multiplications are a fast substitute to function's calls. The relation () can be used on grids of points obtaining a substantial speeding-up in the computation, without loss of precision. We will write below two utilities, hermiteP and getPsi which can help the reader to manage wave functions in an efficient way. In the section [] some tests on the procedure are performed.

Approximate computation of matrix elements

Except for polynomial potentials analytical computation of matrix elements of the potential is generally not available. One has to resort to approximate methods. A possibility is to compute all matrix elements with NIntegrate, this is usually a quite long procedure and for pedagogical purposes we follow a simpler route: the approximation of an integral by Gaussian sums.

The procedure is briefly reviewed in notebook [3], and applied to the Laguerre polynomials in notebook [7].  Here we construct a grid of n points, given by the roots of the n-th Hermite polynomial. With this grid the integrals can be approximated by

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

where (Subscript[w, a]) are the weights. The routine hermiteRoots[n] computes the  roots and weights, the reader is referred to notebook  [3] for more details.  The formula is exact for F[x] polynomial of degree less or equal 2n-1.

The basis functions in the present scheme have the form () ( (Subscript[P, k]) is a polynomial of degree k-1)

(Subscript[f, k])[x]= Exp[-(x^2)/2] (Subscript[P, k])

The matrix elements of the potentials can be approximated by

(Subscript[V, ij])= (∑_(a=1)^n) (Subscript[w, a]) (Subscript[P, i])[(Subscript[x, a])] (Subscript[P, j])[(Subscript[x, a])]V[(Subscript[x, a])/α]

where a possible scale α has been included. In this formula, as in the ones for kinetic term, the additional 1/α factor coming from jacobian from initial variable and the scaled variable x,  has not been included, as it is common to all matrices and cancels in the Schrödinger equation.

Functions

We need in fact only one function, which is just a simplification of a matrix used in notebook [9.1]. The calling sequence is

{p2, x2, x4} = hermiteMatrices[dim, parity]

It produces matrix elements for operators p^2 , x^2, x^4. Parity is 0 or 1 for even and odd states. dim is the dimension of the basis. If parity has any other value full matrices are computed, for both even and odd states, and the output is changed adding matrix elements of x in the last place. The default is parity = -1, i.e. full matrices

{p2, x2, x4,x1} = hermiteMatrices[dim, -1],  or  hermiteMatrices[dim].

hermiteMatrices[dim_,parity_Integer:-1]:=
Module[{x2,p2,x4,x1}, 
Block[{i,n,shift},
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;];x2=SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];x2 =(x2+Transpose[x2]);
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
x4=SparseArray[{{i_,i_}->(3/8)(2(n^2)+2n+1),{i_,k_}/;k-i==shift->((4n+6)/4)(Sqrt[(n+1)(n+2)]),
{i_,k_}/;k-i==2shift->(1/4)(Sqrt[(n+1)(n+2)(n+3)(n+4)])},{dim,dim}];
x4 =  (x4+Transpose[x4]);
If[parity == 0 || parity == 1,Return[Normal/@{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
Normal/@{p2,x2,x4,x1}]];

We give also separately a function which compute only kinetic term, this will be used for non polynomial interactions

hermiteP2[dim_,parity_Integer:-1]:=
Module[{p2}, 
Block[{i,n,shift},
If[parity == 0 || parity == 1,n =2i-2+parity; shift = 1;,n=i-1;shift =2;];
p2 =SparseArray[{{i_,i_}->(1/4)(2n+1),{i_,k_}/;k-i==shift->-(1/2)(Sqrt[(n+1)(n+2)])},{dim,dim}];
p2 = (p2+Transpose[p2]);
Normal[p2]]];

Basis' s functions are given in general by

fBase[k_,x_]:= (1/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,x]Exp[-x^2/2]

For fixed parity we can use the decomposition () and write

fBase[n_,x_,par_]:=(1/(π^(1/4)))(1/(Sqrt[(2^(2n-2+par))(2n-2+par)!]))HermiteH[2n-2+par,x]Exp[-x^2/2]

Wave function computation

Here we just write recursion equation (), the reader usually will not need a direct call to this function. It is compiled for speed and the integer n represents the number of polynomials, i.e. all polynomials from degree 0 to degree n-1.

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

This function computes normalized Hermite polynomials, in the form needed for oscillator wave functions. We write two versions of it, one with arguments [dim, r] gives all polynomials up to degree dim-1.  The second with three arguments, [dim, r, par], gives the first dim polynomials of parity fixed by par, (0 for even polynomials, 1 for odd polynomials).

hermiteP[dim_,xIn_]:= Module[{cf,x,ris,coeff},
x = Flatten[{xIn}];
 coeff=Table[(1/((N[π])^(1/4))) (1/(Sqrt[(2.^(k-1))(k-1)!])),{k,1,dim}]; 
ris= recursiveHermitePolynomials[dim,x]*coeff;
If[Length[xIn]== 0,Flatten[ris],ris]];

hermiteP[dim_,xIn_,par_]:= Module[{cf,dimUp,index,x,ris},
x = Flatten[{xIn}];
dimUp=2 dim - 2 + par+1;
index=Range[par+1,2 dim-2+par+1,2];
ris=hermiteP[dimUp,x][[index]];
ris];

A more efficient version would use recurrence only for even or odd polynomials, we leave to the reader the task of writing such a routine.

Finally getPsi give the wave function for a coefficient list eigv at points r. Parity can be specified, if not all polynomials are considered. This is the default.

getPsi[eigv_,r_,par_: -1]:=
If[par== 1|| par==0, 
(eigv.hermiteP[Length[eigv],r,par])*Exp[-(r^2)/2],
(eigv.hermiteP[Length[eigv],r])*Exp[-(r^2)/2]];

Functions for approximate calculations

hermiteRoots[n] gives roots and weights for the n-points Gaussian integration. The calling sequence is

{xgrid, weights} = hermiteRoots[n]

hermiteRoots[n_]:= Module[{J,δ,γ,radEig,w},
  γ = Sqrt[Range[1,n-1]/2];J = Normal[SparseArray[{Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
radEig =Sort[ Eigenvalues[N@J]];
w=(((2^(n-1))(Sqrt[π]))/(n^2))((n!)/((HermiteH[n-1,radEig])^2));
{radEig,w}]

The function getMatV is an utility routine which compute matrix element for a potential V, once Gaussian points and the rectangular matrix (Subscript[P, k])[(Subscript[x, a])]has been computed. The calling procedure is:

matV = getMatV[α, V, {xgrid, weights}, matBaseHerm]

where

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

Checks

Checks for matrices

We do some simple checks by comparing matrices to explicitly computed matrix elements. We do the check for even and odd states and for full matrices.

All integrals need for checks are of the form

(∫_(-∞)^(+∞))(x^k) Exp[-(x^2)] = (Piecewise[{{Γ(((k+1)/2)), k even}, {0, k odd}}])

We can write a simple substitution rule which avoid the use on Integrate

myInt[funz_]:= Together[Expand[funz]//.{x^a_. Exp[-x^2]->  Gamma[((1+a)/2)](1/2)(1+(((-1))^a)), Exp[-(x^2)]-> Gamma[(1/2)]}];

even

dim=6;parity=0;

{p2,x2,x4}=hermiteMatrices[dim,parity];

p2I = Table[myInt[(∂_x)fBase[i,x,parity](∂_x)fBase[j,x,parity]],{i,1,dim},{j,1,dim}];
x2I = Table[myInt[(x^2)fBase[i,x,parity]fBase[j,x,parity]],{i,1,dim},{j,1,dim}];
x4I = Table[myInt[(x^4)fBase[i,x,parity]fBase[j,x,parity]],{i,1,dim},{j,1,dim}];

{x2== x2I,x4== x4I,p2== p2I}

odd

dim=6;parity=1;

{p2,x2,x4}=hermiteMatrices[dim,parity];

p2I = Table[myInt[(∂_x)fBase[i,x,parity](∂_x)fBase[j,x,parity]],{i,1,dim},{j,1,dim}];
x2I = Table[myInt[(x^2)fBase[i,x,parity]fBase[j,x,parity]],{i,1,dim},{j,1,dim}];
x4I = Table[myInt[(x^4)fBase[i,x,parity]fBase[j,x,parity]],{i,1,dim},{j,1,dim}];

{x2== x2I,x4== x4I,p2== p2I}

full matrices

dim=6;parity=-1;

{p2,x2,x4,x1}=hermiteMatrices[dim,parity];

p2I = Table[myInt[(∂_x)fBase[i,x](∂_x)fBase[j,x]],{i,1,dim},{j,1,dim}];
x2I = Table[myInt[(x^2)fBase[i,x]fBase[j,x]],{i,1,dim},{j,1,dim}];
x4I = Table[myInt[(x^4)fBase[i,x]fBase[j,x]],{i,1,dim},{j,1,dim}];
x1I = Table[myInt[x fBase[i,x]fBase[j,x]],{i,1,dim},{j,1,dim}];

{x2== x2I,x4== x4I,p2== p2I,x1== x1I}

Checks for wave functions

To do a check on a realistic example we will use a coefficient' s vector generate by the solution of an eigenvalue problem. The commands used below will be explained in the next section, what is of interest now is a generation of a coefficient vector and of a grid of points.

In the example we take dim = 100, i.e. a 100 component vector, for a fixed parity, parity = 0. The actual wave function as a multiplicative factor (Sqrt[α]) but this is of no importance here.

α = 2.0; g=1.0;
 dim=100; parity = 0;

Here we generate the vector of coefficients, vec

{p2,x2,x4} =N/@ hermiteMatrices[dim,parity];
Ham=((α^2)/2)p2 + (1/(2 (α^2)))x2 + (g/(2(α^4)))x4;
{eigH,eigV} = Eigensystem[Ham];
order = Ordering[eigH];
eigV=eigV[[order]];
vec = eigV[[1]];

The coefficients describe the ground state of an anharmonic oscillator. To compute and plot the wave function we could simply use the functions fBase. We remember that to pass from rescaled coordinates x to real coordinates, r, we need a the scale factor α, then:

ψA[r_]:= vec.fBase[Range[dim],α r,parity]

We used the fact that fBase can act directly on the range of polynomials. Let us consider the time for computing a single value and a list of values.

rgrid = Range[-5,5,1/10];

{Do[ψA[0.1],{j,1,20}]//Timing,(ψA[#1]&/@rgrid;)//Timing}

The reader can check that adding the attribute Listable at ψA does not help. With these execution times doing a plot can take some time. Much better is a ListPlot

{Plot[ψA[x],{x,-5,5},PerformanceGoal->"Speed"]//Timing,ListPlot[Transpose[{rgrid,ψA[#1]&/@rgrid}],Joined-> True]//Timing}

Without the option PerformanceGoal -> "Speed", Plot would take almost 7 times longer to execute.

Compiled function

Now let us use the function getPsi which use a recursive definition of Hermite polynomials.

ψR[r_]:= getPsi[vec,α r,parity]

{Do[ψR[0.1],{j,1,20}]//Timing,(ψR[rgrid];)//Timing}

We see that it is  faster by an order of magnitude on single values and two order of magnitudes on lists. This has some consequences on  plotting times:

{Plot[ψR[x],{x,-5,5},PerformanceGoal-> "Speed"]//Timing,ListPlot[Transpose[{rgrid,ψR[rgrid]}],Joined-> True]//Timing}

The two functions are equal within machine precision :

Max[Abs[ψR[rgrid] -( ψA[#1]&/@rgrid)]]

The reader can check the same relations for odd parity and for full matrices.

Polynomial 's matrices

In the computation with the Gauss sums one has to compute the matrix H[k,rgrid] i.e. the values of (normalized) Hermite polynomials of a grid, the grid of Gauss points. This computation has to be done only once for each dimension, then in principle one can use also the analytical form of the polynomials without much wasting of time, but for consistency we present here the computation with the recursive procedure.

dimGauss = 30;
{xgrid,weights}= hermiteRoots[dimGauss];

matBaseHermA= N/@Table[((((-1))^(k-1))/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,xgrid],{k,1,dim}];//Timing

matBaseHerm = hermiteP[dim,xgrid];//Timing

Max[Abs[Flatten[matBaseHermA-matBaseHerm]]/Abs[Flatten[matBaseHerm]]]

Dimensions[matBaseHermA]

Clear["Global`ψ*", rgrid, xgrid, α, g, p2, x2, x4, Ham,dim,dimGauss,matBaseHermA,weights,matBaseHerm,eigH,eigV,vec,order,parity,p2I,x2I,x4I]

Example: harmonic and anharmonic oscillators

Harmonic oscillator

The first example is a bit trivial but can be useful as a first application. Take a harmonic oscillator with ω ≠ 1. We are using a basis with functions built on a harmonic oscillator with ω = 1. In scaled variables the Hamiltonian takes the form

(Subscript[H, ω]) = ((α^2)/2)(p^2) + ((ω^2)/(2(α^2)))(x^2)

It is clear that for (α^2)= ω we get a Hamiltonian (Subscript[H, ω]) = ω (Subscript[H, 1])and in this case (Subscript[H, ω]) is diagonal, i.e. we found the eigenvalues and eigenvectors. Let us see how this works in a variational approach. We will work with a low dimensional basis to allow the reader, if he wants, to print explicit forms for matrices.

We choose only even states and put dim = 10;

dim = 10; parity = 0;

Then we compute our matrices. In this simple case the potential is homogeneous. We obtain its matrix elements directly from computed powers of x

{p2,x2,x4} = hermiteMatrices[dim,parity];

The effective hamiltonian, in our scaled variables, takes the form

H = ((α^2)/2)p2 + ((ω^2)/(2 (α^2)))x2

We choose ω = 4 and plot the eigenvalues of H as a function of α.

ω = 4;

energy[α_?NumberQ,state_Integer]:= Module[{},
Sort[Eigenvalues[((α^2)/2)p2 + ((ω^2)/(2 (α^2)))x2]][[state]]]

Plot[energy[α,1],{α,0.1,5},AxesLabel-> {"α","E"}]

An extremely flat region at an eigenvalue of 2 is seen, which is correct, as ω/2 = 2 in our example.

We can use FindMinimum for locating a best value for α.

min=FindMinimum[energy[α,1],{α,1,3.5}]

We get an idea of the precision, in this case, from the conversion with NumberForm :

NumberForm[min[[1]],20]

The error is about 4 10^-15, and this is a typical precision for this kind of approach.

In this case, as in almost all other cases in the variational approach, it is better to use FindMinimum in the form

FindMinimum[f[x], {x, (Subscript[x, 1]),(Subscript[x, 2])}]

In this form the use of derivatives is avoided, and this is necessary both because eigenvalues are not perfectly smooth functions of α, due to rounding errors, and because the derivatives are almost zero if α is in a "safe" region, i.e. where eigenvalues asymptotically do not depend on α. It can be interesting to look into FindMinimum to see how many function call was made and how the search explores the interesting zone. This can be done in the following way (we use a procedure described in Mathematica help ):

This loads a package which contains some utility functions:

<<Optimization`UnconstrainedProblems`

FindMinimumPlot[energy[α,1],{α,1,3.5}]

Anharmonic oscillator

The anharmonic oscillator can be considered with almost no changes.

dim = 10; parity = 0;

{p2,x2,x4} =N/@ hermiteMatrices[dim,parity];

We use the potential

V[r] = (1/2)(r^2) + (1/2)g (r^4)

g=1;

We redefine our energy function to take into account the anharmonic term

energy[α_?NumberQ,g_,state_Integer]:= Module[{},
Sort[Eigenvalues[((α^2)/2)p2 + (1/(2 (α^2)))x2 + (g/(2(α^4)))x4]][[state]]]

and proceed as before. This is the α dependence of ground state energy

Plot[energy[α,g,1],{α,0.1,5},AxesLabel-> {"α","E"}]

This is a more refined estimation of the eigenvalue

min=FindMinimum[energy[αvar,g,1],{αvar,1,3.5}]

{enex4,α} = {min[[1]],αvar/.min[[2]]}

NumberForm[enex4,20]

An in - depth discussion of the model, the errors and high precision results is given in following notebooks.

Wave functions

Here we give just an example of how the wave functions can be computed. To make a comparison with a known case we use harmonic oscillator and choose the scale factor α = 2.

dim=50; parity=0; α = 2.0;
{p2,x2,x4} = hermiteMatrices[dim,0];
Ham =((α^2)/2)p2 + (1/(2 (α^2)))x2 ;
{eigH,eigV} = Eigensystem[Ham];
order = Ordering[eigH];
eigV=eigV[[order]];
vec = eigV[[1]];

To write the numerical wave function we have to remember two points

The theoretical wave function of the ground state is simply the first function of our basis fBase, let us write this explicitly:

ψTh[r_]:= fBase[1,r]

The numerical wave function is obtained using the function getPsi, at rescaled coordinates:

ψ[r_]:= (Sqrt[α]) getPsi[vec,α r,0]

Plot[{ψTh[x],ψ[x]},{x,-5,5},PerformanceGoal-> "Speed"]

The two functions are indistinct within the resolution of the graph.

Generic potential

For all (but simple polynomial) potentials the matrix elements must be evaluated numerically. As stated in the introductory section we use Gaussian sums for the approximation. The function hermiteRoots and getMatV given in the section Functions, computes the Gaussian grid and the matrix elements. We need only the matrix for kinetic term, to this purpose we use the simplified function hermiteP2, which produces the matrix elements of operator (p^2).

Here is the modification for the function energy which computes the eigenvalues. Its input is more complicated than before. It needs the variables:

energy[α_?NumberQ,state_,V_,{xgrid_,weights_},matBaseHerm_]:= Module[{matV,dim},
dim = Length[matBaseHerm];
matV=getMatV[α,V,{xgrid,weights},matBaseHerm];
Sort[Re[Eigenvalues[((α^2)/2) p2+matV ]]][[state]]]

Let us choose the potential and the dimension we want to use

g=1;
V[x_]:= (1/2)(x^2)+(1/2)g (x^4)

dim=20;

Here is kinetic term

p2= N/@hermiteP2[dim];

Here we specify how big  the Gauss grid must be. We will use n = dim + shiftDim

shiftDim = 10;

Here are the roots and weights :

{xgrid,weights}= hermiteRoots[dim+shiftDim];

and finally the matrix matBaseHerm:

matBaseHerm = hermiteP[dim,xgrid];

This is the alternative form which can be used;

matBaseHerm= N/@Table[((((-1))^(k-1))/(π^(1/4))) (1/(Sqrt[(2^(k-1))(k-1)!]))HermiteH[k-1,xgrid],{k,1,dim}];

An example of computation :

energy[1.,g,V,{xgrid,weights},matBaseHerm]

We plot the energy eigenvalues as a function of α to detect a stability zone in α.

Plot[energy[α,g,V,{xgrid,weights},matBaseHerm],{α,0.1,3}]

It is seen that for 1 and 2.5 the eigenvalue  is very stable, if we want a better determination, with this value of dim, we can use FindMinimum

min=FindMinimum[energy[αvar,g,V,{xgrid,weights},matBaseHerm],{αvar,1,3.0}]

The result is in perfect agreement with the one found above, and this must not be a surprise as Gaussian integrals are exact for this interaction:

NumberForm[min[[1]] - enex4,20]