Discretization method

Quantum Mechanics
K. Konishi, G. Paffuti

Study the connection between the discretization methods and the variational principle.

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

Finite differences

Short review

In the first chapters of this book we used sometimes a very simple approximation to solve numerically the Schrödinger equation. We imagine our (one dimensional) system enclosed in a box of length L,  [xL,xR].  The Schrödinger (in natural units) equation

-(1/2)((d^2)/(dx^2))ψ[x] + V[x]ψ[x]= E ψ[x]

was supplied with the boundary conditions ψ[xL]=ψ[xR] = 0, which mimic the request of asymptotically vanishing solutions. The interval was discretized in n+1 intervals, with n internal points: due to boundary conditions only ψ at these points are unknown quantities.  The grid has the form

(Subscript[x, k]) = xL + k h; k =1,...,n;   h = (L/(n+1))

With a two points discretization of second derivative on the grid:

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

equation () was transformed in a linear system, with unknown quantities ψ[(Subscript[x, s])]≡ (Subscript[ψ, s]):

-(1/(2(h^2)))(Subscript[D2, ij])(Subscript[ψ, j]) + (Subscript[V, i]) (Subscript[ψ, i]) = E (Subscript[ψ, i]);   or   -(1/(2(h^2)))D2.ψ +V. ψ = E ψ;

where (Subscript[V, i])= V[(Subscript[x, i])], V denote the corresponding diagonal matrix and

D2 =( ({{"-2", "1", "0", ".", ".", ".", "."}, {"1", "-2", "1", ".", ".", ".", "."}, {"0", "1", ".", ".", ".", ".", "."}, {".", ".", ".", ".", ".", ".", "."}, {".", ".", ".", ".", ".", "1", "0"}, {".", ".", ".", ".", "1", "-2", "."}, {".", ".", ".", ".", "0", ".", "-2"}}))

Details and examples can be found in the notebook [NB-3.1-discretization.nb].

Variational version

We now show that an approximation very similar to () can be found from within a variational approach. In general we have to find stationary points of the quadratic functional (the energy is introduced as a Lagrange multiplier)

Q[ψ] = (1/2)∫(((∇ψ))^2) + ∫V[x] (ψ^2) - E ∫(ψ^2)

The idea is to choose a set of functions such that the derivatives pick up the values of the function at a predefined grid of points, which we assume uniform for simplicity.

Consider functions of the following type:

(Subscript[T, i])[x] =(({{"     0", "x ≤ \!\(\*SubscriptBox[\(x\), \(i - 1\)]\)  "}, {"\!\(\*FractionBox[\(x - \*SubscriptBox[\(x\), \(i - 1\)]\), \(\*SubscriptBox[\(x\), \(i\)] - \*SubscriptBox[\(x\), \(i - 1\)]\)]\) = \!\(\*FractionBox[\(x - \*SubscriptBox[\(x\), \(i - 1\)]\), \(h\)]\)", "\!\(\*SubscriptBox[\(x\), \(i - 1\)]\) ≤ x ≤ \!\(\*SubscriptBox[\(x\), \(i\)]\)  "}, {"\!\(\*FractionBox[\(\*SubscriptBox[\(x\), \(i + 1\)] - x\), \(\*SubscriptBox[\(x\), \(i + 1\)] - \*SubscriptBox[\(x\), \(i\)]\)]\) = \!\(\*FractionBox[\(\*SubscriptBox[\(x\), \(i + 1\)] - x\), \(h\)]\)", "\!\(\*SubscriptBox[\(x\), \(i\)]\) ≤ x ≤ \!\(\*SubscriptBox[\(x\), \(i + 1\)]\)  "}, {"     0", "\!\(\*SubscriptBox[\(x\), \(i\)]\) ≤ x  "}, {"    ", "   "}})

These functions are centered at grid points and vanishes after neighbors points. In Mathematica it can be constructed with the function Cases or with Piecewise (see Help).

This is a grid (we add extrema points for simplicity)

nP=10;
Table[xgrid[i]=i  h,{i,0,nP+1}];

In the function the index i appear shifted because vector's index start from 1

fT[i_,x_]:= Piecewise[{{0,x≤ xgrid[i-1]},
{((x-xgrid[i-1])/(xgrid[i]-xgrid[i-1])),xgrid[i-1]< x ≤ xgrid[i]},
{((xgrid[i+1]-x)/(xgrid[i+1]-xgrid[i])),xgrid[i]< x ≤ xgrid[i+1]},
{0,xgrid[i+1]< x}}];

Plot[Table[fT[i,x]/.h-> 1,{i,4,6}],{x,0,10},PlotRange-> {All,{0,1.5}},Filling->Axis]

A function f[x] is approximated by a broken line, fitting f[x] at grid' s points.

f[x_]:= (3/10)(x^2) -(3/95) (x^3)

Block[{h=1},
Plot[Sum[f[xgrid[i]](fT[i,x]),{i,1,nP}],{x,0,nP},Epilog-> Point/@ Table[{(xgrid[i]),(f[xgrid[i]])},{i,1,nP}]]]

In the general framework of variational principle we are taking a trial function of the form

ψ[x] = (∑_(s=1)^n)(Subscript[ψ, s]) (Subscript[T, i])[x]

Her we call (Subscript[ψ, i])the coefficient as they are just the values of ψ at the grid point.

In the interval (Subscript[x, s]), (Subscript[x, s+1])a function f and its derivative take the form (remember that (Subscript[f, 0])= (Subscript[f, N+1])= 0)

f[x]∼(Subscript[f, s]) + ((x-(Subscript[x, s]))/h)((Subscript[f, s+1])-(Subscript[f, s]));  Df ∼ (1/h)((Subscript[f, s+1])-(Subscript[f, s]));

On this set of functions the derivative part of the functional Q, eq.() take the form (the derivatives are constant and the integral in each interval is their values times the extension of the interval, i.e. h):

(1/2)∫(((∇ψ))^2) = (1/2)(∑_(s=0)^n) (1/(h^2))h ((((Subscript[ψ, s+1])-(Subscript[ψ, s])))^2) ≡ - (h/(2 (h^2))) (Subscript[ψ, i])(Subscript[D2, ij])(Subscript[ψ, j])

D2 is the matrix defined above.

Let us now consider the potential and norm term in Q. The potential term gives the contribution:

(∑_ij) (Subscript[ψ, i])(Subscript[ψ, j]) ∫ V[x] (Subscript[T, i])[x] (Subscript[T, j])[x] 𝕕x

As basis functions with indexes |i-j| ≥ 2 do not overlap, the matrix V is tridiagonal. A further simplification can be obtained if one uses trapezoidal approximation for integral. On the interval (Subscript[x, s])(Subscript[x, s+1]), of length h, the integral is approximated by

(∫_(Subscript[x, s])^(Subscript[x, s+1])) F[x]𝕕x ∼ (h/2)( F[(Subscript[x, s+1])] + F[(Subscript[x, s])])

and on the whole interval (our functional vanish at end points due to (Subscript[T, i])factors):

∫F[x] = h (∑_(s=1)^n)F[(Subscript[x, s])]

Due to the supports of the functions (Subscript[T, i])this implies

(Subscript[V, ii])= h V[(Subscript[x, i])]  and    (Subscript[V, ij])= 0     for    i≠j,

as the reader can easily check.

The same procedure can be applied in principle to (ψ^2)term, but the integrals can also be computed analytically and one easily found that

h (Subscript[N, ij])=∫(Subscript[T, i])[x](Subscript[T, j])[x]𝕕x =h (Piecewise[{{2/3, i=j}, {1/6, |i-j]=1}}])

i.e. we have again a tridiagonal matrix with 2/3 on the main diagonal and 1/6 on adjacent diagonals. As a check:

Table[Integrate[Evaluate[fT[i,x]fT[j,x]],{x,0,(nP+1)h},Assumptions-> h>0],{i,1,3},{j,1,3}]//MatrixForm

In conclusion the quadratic form becomes

Q = - (h/(2 (h^2))) (Subscript[ψ, i])(Subscript[D2, ij])(Subscript[ψ, j]) + (Subscript[V, ij]) (Subscript[ψ, i])(Subscript[ψ, j]) - E h (Subscript[N, ij]) (Subscript[ψ, i])(Subscript[ψ, j])

The variational principle implies than:

- (1/(2 (h^2))) (Subscript[D2, ij])(Subscript[ψ, j])+ (1/h)(Subscript[V, ij]) (Subscript[ψ, j]) = E (Subscript[N, ij])(Subscript[ψ, j])

which coincides with eq.() if integral are performed with a trapezoidal rule.

Potential

If the exact form of the potential matrix is required, it is sufficient to compute (Subscript[V, ss ])and  (Subscript[V, s,s+1]) due to the Hermiticity and tridiagonal form of the matrix. For a grid x = xL + s h one can perform the change of variables

x = xL + s h + z h;   (Subscript[x, s-1])≤ x≤ (Subscript[x, s+1]) ⇒ -1 ≤ z ≤ 1

and the raising and lowering part of (Subscript[T, s]) , (Subscript[T, s+1])functions.  For this interval and this variable, they are given  by

(Subscript[T, s]):   1+z;  -1≤ z ≤ 0;      1-z ; 0≤ z ≤ 1;    (Subscript[T, s+1]) :  z  for   0≤ z ≤ 1

It follows then that

(Subscript[V, ss])= h ((∫_(-1)^0)V[xL+s h + z](((1+z))^2)dz  + (∫_0^1)V[xL+s h + z](((1-z))^2)dz)
(Subscript[V, s,s+1])= h (∫_0^1)V[xL+s h + z](1-z)z dz

as a check let us compute the matrix N, which correspond to V= 1.

V[x_]:= 1;

diag =h Table[Integrate[V[xL+s h + z h](((1+z))^2),{z,-1,0}]+
Integrate[V[xL+s h + h z](((1-z))^2),{z,0,1}],{s,1,4}]

upperdiag = h Table[
Integrate[V[xL+s h + h z](1-z)z,{z,0,1}],{s,1,4}]

As a second example which will be useful later we can consider a harmonic potential and compute the formal expression for matrix elements: (we do not write the overall h term)

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

diagonal term, s - s

Integrate[V[xL+s h + z h](((1+z))^2),{z,-1,0}]+
Integrate[V[xL+s h +h z](((1-z))^2),{z,0,1}]

off diagonal term, s - s + 1

Integrate[V[xL+s h + h z](1-z)z,{z,0,1}]

Variational parameters

Let us ask a legitimate question : which are the variational parameters in the above scheme? The choice of the basis is a choice for the form of trial functions, and as usual we expect that for smooth functions the set (Subscript[T, i]) give convergence. But even with an infinite numbers of T - functions we do not  recover the exact answer to our problem: we are always in a finite volume, the length L of the box. This is our real parameter and obviously L→∞ is the stationary point.

How is the approach to the this point? This is an interesting problem by itself, we ask how the infinite volume limit is approached in a quantum system. We have not a simple and general answer to this question, so let us make some comments on this problem.

In a finite volume all eigenvalues of our matrices describe obviously bound states, and the spectrum is always discrete. This is even true for an infinite set of basis functions, i.e. for the exact solution in a finite volume. Take the free particle as an example: the eigenvalues are

(Subscript[E, n])= (((ℏ^2)(π^2))/(2m))((n^2)/(L^2))

We see that if one take the limit L→∞ at fixed n, eigenvalues go to zero. The correct, continuum, spectrum is obtained in the double limit L→∞, n→∞, n fixed, which means to keep fixed values for continuum momentum.

The question is different for bound states. If there are, for example, bound states for negative energy, then the gap between these states and the continuum must survive the limit. An analogous mechanism must work for cases like harmonic oscillator where we have only a discrete spectrum. In other words  the limit L→∞ must exists and one has

(Underscript[Lim, L→∞]) (Subscript[E, n])[L] = (Subscript[E, n])

A support to this claim is the following. Let us suppose that we consider the Hamiltonian at a scale λ (our box size) H[λ], and call E[λ] its eigenvalues. An infinitesimal change of scale is described by the infinitesimal generator D, where

λ = (e^α);   D = (ℏ/(2i))(x (∂/(∂x))+(∂/(∂x)) x) = (1/2)(x p + p x).

From Feynman Hellman theorem we have

((δE[λ])/(δα)) = ⟨ψ|((δH)/(δα))|ψ⟩ ∼ ⟨ψ|[D,H]|ψ⟩

We will see in chapter [] that commutators of D in finite volume can give some problem, but in the infinite volume limit everything must be well defined. The right hand side of () vanishes on stationary states and this confirm us that L→∞ is the stationary point. From the numerical point of view, especially in the local version of the discretization,  things can be quite involved: if one go to large L at fixed number of points in the grid, the spacing h = L/n become large and the error in the discretization procedure relevant. A better way to approach the limit is to take h fixed, depending somewhat on the required precision, while L gets large. An example will be shown below.

An example of variational calculation

Let us take as an example the simplest system, i.e. a harmonic oscillator.

First we present  a step by step procedure, than we will write a function to do systematic analysis on L dependence and on nP dependence(the number of grid points).

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

The potential is

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

nP, L are respectively the number of internal grid points and the length of our box, taken between -L/2 and L/2.

nP = 50; L = 10;
h = L/(nP+1);
xL = -L/2;
xgrid =xL + Range[nP]h;

The matrix D2  is always the same:

D2=Normal@SparseArray[{Band[{1,1}]-> -2,Band[{1,2}]-> 1,Band[{2,1}]-> 1},{nP,nP}];

For matrices D and N we have two possibilities, diagonal or tridiagonal, we call them matVD, matV, matND, matN

matVD = N@DiagonalMatrix[V[xgrid]];
matND = IdentityMatrix[nP];
matN = N@Normal@SparseArray[{Band[{1,1}]-> 2/3,Band[{1,2}]-> 1/6,Band[{2,1}]-> 1/6},{nP,nP}];

For the potential we use the formulas ()

diag = Table[((h^2)/30)+(((h^2) (s^2))/3)+((2 h s xL)/3)+((xL^2)/3),{s,1,nP}];
upperdiag = Table[((h^2)/40)+(((h^2) s)/12)+(((h^2) (s^2))/12)+((h xL)/12)+((h s xL)/6)+((xL^2)/12),{s,1,nP-1}];
matV = N@Normal@SparseArray[{Band[{1,1}]-> diag,Band[{1,2}]-> upperdiag,Band[{2,1}]-> upperdiag},{nP,nP}];

We will form 3 different Hamiltonians, starting from the simple discretization to the full fledged variational version and find corresponding eigenvalues (the first 5)

numStates = 5;

Sort[Eigenvalues[-(1/(2(h^2)))D2 + matVD]][[Range[numStates]]]

Sort[Eigenvalues[{-(1/(2(h^2)))D2 + matVD,matN}]][[Range[numStates]]]

Sort[Eigenvalues[{-(1/(2(h^2)))D2 + matV,matN}]][[Range[numStates]]]

We see that correctly the third result is an upper limit for the correct values, while the first one, for this model, approach the result from below.

There is some, marginal, difference between second and third results, indicating some small effect due to trapezoidal approximation of integrals. It is worth to say that within our approximation the difference is not crucial. The reader can judge on the accuracy, with this relatively small basis of 50 elements. In variational calculation based on other methods the accuracy with this number of basis element will be much better, but we have to stress the simplicity in the code and the universality of the approach.

Now let us compare systematically the point-wise approximation and the variational (the third hamiltonian).

The function getEigenvalues produces to list of eigenvalues, for the discretized and for the variational version of the problem. Each list contain a number numStates of eigenvalues, the value of this parameter is passed to the function, we will choose to compute 5 states.

getEigenvalues[nP_,L_,numStates_]:= Module[{xgrid,matV,D2,h,matN,diag,upperdiag,xL,matVD,enePoint,eneVar},
h = L/(nP+1);
xL = -L/2;
xgrid =xL + Range[nP]h;
D2=N@SparseArray[{Band[{1,1}]-> -2,Band[{1,2}]-> 1,Band[{2,1}]-> 1},{nP,nP}];
matVD = N@SparseArray[{Band[{1,1}]-> V[xgrid]}];
matN = N@SparseArray[{Band[{1,1}]-> 2/3,Band[{1,2}]-> 1/6,Band[{2,1}]-> 1/6},{nP,nP}];
diag = N@Table[((h^2)/30)+(((h^2) (s^2))/3)+((2 h s xL)/3)+((xL^2)/3),{s,1,nP}];
upperdiag =N@ Table[((h^2)/40)+(((h^2) s)/12)+(((h^2) (s^2))/12)+((h xL)/12)+((h s xL)/6)+((xL^2)/12),{s,1,nP-1}];
matV = N@SparseArray[{Band[{1,1}]-> diag,Band[{1,2}]-> upperdiag,Band[{2,1}]-> upperdiag},{nP,nP}];
enePoint=Sort[Eigenvalues[-(1/(2(h^2)))D2 + matVD,-numStates,Method-> Arnoldi]];
Off[Eigenvalues::"argpd"];
eneVar=Sort[Eigenvalues[{-(1/(2(h^2)))D2 + matV,matN},-numStates,
Method-> {"Arnoldi",BasisSize-> 30}]];
On[Eigenvalues::"argpd"];
{enePoint,eneVar}]

L fixed

We first get result for fixed L. As the typical length of a harmonic oscillator grows like (Sqrt[n]), n is the level, we expect that L=10 is large enough for the first 5 states.

numStates=5;
L=10;

eneLfixed = Table[getEigenvalues[nP,L,numStates],{nP,50,350,50}];

The matrix obtained ha 3 indexes, the first run on nP, the second is 1 or 2 and select the two set of data, the last run from 1 to numStates and select the eigenvalue. To compare the bound state energies:

discretized Hamiltonian

eneLfixed[[All,1,1]]

Variational Hamiltonian

eneLfixed[[All,2,1]]

We confirm that variational results are upper bound while in this case point - wise approximation gives lower bounds.

It is interesting that their mean has a much better precision than individual results, if this would be true for general models it would give an even better shape to this kind of approximation:

NumberForm[(1/2)(eneLfixed[[All,1,1]]+eneLfixed[[All,2,1]]),10]

Theoretical energies are well known :

eneTh = 0.5+Range[0,numStates];

If the errors are plotted as a function of grid points we see that to reach the same level of precision the m-th state one needs a number of steps growing with m, and this is one of the bad features of the method: higher states are difficult to catch.

ListPlot[Table[Log[Abs[eneLfixed[[All,1,k]]-eneTh[[k]]]],{k,1,5}],
PlotRange-> All]

The situation is similar for variational calculation, in this basis.

Changing L

We  select some L, starting from obviously too small ones, to show how the asymptotic limit is approached

varL = Table[k,{k,3,15,0.5}];

Here we compute the eigenvalues taking h fixed, this means to have a number of grid points proportional to L:

eneLvariable = Table[getEigenvalues[Round[100*varL[[k]]],varL[[k]],numStates],{k,1,Length[varL]}];

As we are interested to errors we collect differences between computed energies and theoretical result. We take as an example the variational computation:

data=Table[Transpose[{varL,Abs[eneLvariable[[All,2,k]]-eneTh[[k]]][[1;;Length[varL]]]}],{k,1,5}];

ListLogPlot[data[[1]],PlotRange-> {All,{(10^(-6)),0.2}}]

It is seen that we reach an asymptotic plateau, the error is not zero, this depends on our h. In this computation h∼ L/nP ∼ 1/100. The  error expected for discretization of the derivative is 𝒪((h^2)) and this explains the graph. Let us make the same computation with a smaller h:

eneLvariable = Table[getEigenvalues[Round[300*varL[[k]]],varL[[k]],numStates],{k,1,Length[varL]}];

data=Table[Transpose[{varL,Abs[eneLvariable[[All,2,k]]-eneTh[[k]]][[1;;Length[varL]]]}],{k,1,5}];

ListLogPlot[data[[1]],PlotRange-> {All,{(10^(-8)),0.2}}]

With a factor 1/3 on h we obtain an order of magnitude better results. We stress that for L=10 this computation has used a basis of 3000 elements, and this is not very economic from computational point of view.

We leave to the reader the task of studying excited states.