Finite differences

Quantum Mechanics
 K. Konishi, G. Paffuti

Solve numerically  the eigenvalue problem

H ψ = E ψ ;      H = ((p^2)/(2m)) + V[x];

by discretizing the coordinate x. Check the method in the harmonic oscillator problem.

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

General remarks

Let us consider the usual Schrödinger's equation (in natural units m = ℏ = 1):

-1/2 ψ''[x] + V[x] ψ[x] = E ψ[x];                    b.c.:   ψ[xL] = ψ[xR] = 0;

The equation has been defined in the interval [xL.xR] and the usual Dirichlet boundary conditions have been assumed. The problem is to find eigenvalues and eigenvectors of eq.(). 
One simple procedure is to discretize the interval [xL,xR], i.e. to consider a grid of step h:

x[j] = xL + j*h;  j = 1,2,... ,nP;    Internal points in [xL,xR].

The unknowns in the problem become the values of ψ at the internal points of the grid, the values at the boundary are fixed. nP is the number of internal points and the grid has nP+1 intervals of width h. On the grid we can approximate the second derivative by finite differences:

With F[x[j]] = f[j],    F''[x[j]] ≃ ( f[j+1] - 2 f[j] + f[j-1])/(h^2);

and the Schrödinger equation becomes a system of linear homogeneous linear equations

-(1/(2(h^2)))(ψ[j+1]-2ψ[j] + ψ[j-1]) + V[j] ψ[j] = E ψ[j]

We have nP equations and using ψ[0] = ψ[nP+1] = 0:

-(1/(2(h^2)))(ψ[2] - 2 ψ[1]               ) + V[1]ψ[1] = E ψ[1];                                                  Eq. 1;
-(1/(2(h^2)))(ψ[3] - 2 ψ[2] + ψ[1]) + V[2]ψ[2] = E ψ[2];                                                   Eq. 2;
...  ...   ...   ... ;
-(1/(2(h^2)))(ψ[nP] - 2 ψ[nP-1] + ψ[nP-2]) + V[nP-1]ψ[nP-1] = E ψ[nP-1];      Eq. nP-1;-(1/(2(h^2)))(             - 2 ψ[nP]       + ψ[nP-1]) + V[nP]ψ[nP] = E ψ[nP];                            Eq.  nP;

We can write compactly the equations in the matrix form. Let us call D2 the matrix representing second derivative: from the system above it is a tri-diagonal matrix, with -2 on the main diagonal and +1 on diagonals shifted by one place. The potential term is represented by a diagonal matrix.

-(1/(2(h^2))) D2.ψ + V.ψ = E ψ;                         D2 = (({{-2, 1, 0, 0, 0}, {1, -2, 1, 0, 0}, {0, 1, -2, 1, 0}, {0, 0, 1, -2, 1}, {0, 0, 0, 1, -2}}));  (nP = 5);

Infinite domain

In usual quantum mechanics applications the interval [xL,xR] is unbounded. If the eigenfunctions decay sufficiently fast at large distances we can put the system in a box, for example:

xL = - L;         xR = L;

solve the system and then check for numerical stability as L become larger. This procedure is called domain truncation.

Introduction to domain truncation methods

In the example below we take as a check the harmonic oscillator potential V[x] = x^2/2, most of the methods are valid for a general potential. With this potential the domain of definition is [-∞,+∞]. In some cases (eg. radial equations in 3-dim problems) the interesting domain is [0,∞]. The methods can be trivially extended to this case.

V[x_]:= x^2/2;

Eigenvalues: simplest use

Define the grid, we put L = 10 as limit value for the interval.

L = 10; xL = -L/2; xR = L/2;
nP =105;
h = N[(xR - xL)/(nP+1)]; 
xP= xL + h* Range[nP];

N.B.:  nP internal points means nP+1 intervals. If nP is odd the point x = 0 belongs to the grid.

Define the kinetic matrix (D2) and the potential matrix (matV):

D2=SparseArray[{Band[{1,1}]-> -2,Band[{2,1}]-> 1,Band[{1,2}]-> 1},{nP,nP}];
matV = SparseArray[Band[{1,1}]-> V[xP],{nP,nP}];

The discretized  form of the hamiltonian is:

Ham = - 1/(2 (h^2)) D2 + matV;

Eigenvalues and eigenvectors can be computed separately or together using Mathematica functions Eigenvalues and Eigensystem.

We can compute all eigenvalues and choose those we want are interested to:

eigAll = Sort[Eigenvalues[Ham]];
eigAll[[Range[1,5]]]

Mathematica note

Function SparseArray 
Most elements of our matrices are 0. SparseArray gives an efficient way to store such objects, keeping only elements different from 0. Mathematica help gives detailed informations on the subject, see also Advanced Documentation. 
We used Mathematica 6 function Band[{i,j}] to fill the diagonal starting from position {i,j}.
Let us only remember that a standard way to fill the matrix is through patterns:  eg. {i_,i_} ->  -2 means that diagonals elements are -2,                                   the pattern {i_,j_}/;j-i== 1 means: for each i, j under the condition j-i = 1 the element is ... i.e. we define upper main diagonal etc.
For the matrix matV we can define the whole set of diagonal indices and then apply the rule indDiag-> V[x]. This is faster then {i_,i_} -> V[x[[i]]], as the reader can check. 
Using these instructions an alternative form for D2 and matV is

D2 = SparseArray[{{i_, i_} -> -2, {i_, j_} /; j - i == 1 -> 1,{i_, j_} /; i - j == 1 ->  1}, {nP, nP}];
indDiag = Table[{i, i}, {i, 1, nP}];
matV = SparseArray[indDiag -> V[xP]];

The function Normal transforms sparsearray structures in usual matrices.
The function SparseArray is available only with Mathematica 5. Users with previous version of the software can define matrices using Table, and / or other Mathematica functions. Eg.  
<< LinearAlgebra`MatrixManipulation`;
D2 = ZeroMatrix[nP];
Table[D2[[i, i]] = -2, {i, 1, nP}];
Table[D2[[i, i + 1]] = 1, {i, 1, nP - 1}];
Table[D2[[i + 1, i]] = 1, {i, 1, nP - 1}];

Function Eigenvalues 
Eigenvalues[M] compute the eigenvalues of M, it has several extensions to be used in other notebooks and the reader is referred to Mathematica Help for more details.
If M is a sparse matrix Eigenvalues use a recursive and efficient method to find some of the eigenvalues, if the matrix is full or more than 20% of the eigenvalues are requested it is more efficient to switch to a general method, capable of dealing also with arbitrary precision numbers, based on LAPACK routines. In the last instruction we call the function with a sparse matrix and request all the eigenvalues, then automatically Mathematica switched on the full-matrix method and gives us a warning.
To avoid the warning we can either transform our matrix to a Normal form, i.e. non sparse, or request only a few eigenvalues, the last method is usually vastly more efficient. If a number k of eigenvalues is requested the syntax is 
Eigenvalues[M, k]
The eigenvalues are  given in  decreasing absolute value. and k is the k-th eigenvalue. If we want the first k smallest eigenvalues we write a negative entry, i.e.
Eigenvalues[M, -5]
compute the lowest 5 eigenvalues. Below we compare the two methods. The times can change significantly depending on the computer you use, the processor etc. as usually linear algebra routines are optimized at a very low level.

Full matrix time :

eig5 =Sort[ Eigenvalues[Normal[Ham],-5]]//Timing

Sparse matrices algorithm

eig5 =Sort[ Eigenvalues[Ham,-5]]//Timing

The method used by Mathematica is called Arnoldi and usually converges quite fast and to the correct result. In case of troubles you can use some sub- options, one is 
Eigenvalues[M, k, Method-> {"Arnoldi", "Shift" -> value}
where value is an approximate value for one of the eigenvalues. A second option is to enlarge the dimension of the basis, BasisSize used internally by the methods to compute eigenvalues and/or the number of iterations used, MaxIterations. A couple of examples follow:

eig50 =Sort[ Eigenvalues[Ham,-50,Method-> {"Arnoldi","Shift"-> 0}]]//Timing

eig50 =Sort[ Eigenvalues[Ham,-50,Method-> {"Arnoldi","Shift"-> 0,BasisSize-> 70}]]//Timing

Eigenfunctions

The Mathematica function Eigensystem can be used to compute eigenfunctions. The calling sequence is {eigs,eigV} = Eigensystem[matrix]  eigs is the list of eigenvectors while eigV is the list of eigenfunctions. With Eigensystem it is possible to use the same options used with Eigenvalues. To be a bit more realistic we use a larger number of points, nP = 1001, i.e. matrices formally will be 1000 × 1000, and try to compute the first 50 eigenvalues (we do not expect good results for high eigenvalues).

L = 20; xL = -L/2; xR = L/2;
nP =1001;
h = N[(xR - xL)/(nP+1)]; 
xP= xL + h* Range[nP];

D2=SparseArray[{Band[{1,1}]-> -2,Band[{2,1}]-> 1,Band[{1,2}]-> 1},{nP,nP}];
matV = SparseArray[Band[{1,1}]-> V[xP],{nP,nP}];

Ham = - 1/(2 (h^2)) D2 + matV;

The output of  Eigensystem is a list, whose first elements is the vector of eigenvalues, the second one the set of corresponding eigenvectors.

{ene,vecPsi} = Eigensystem[Ham,-50];

Usually  ordered eigenvalues are required. We will do this using the function Ordering, which gives the canonical order of a list.

order = Ordering[ene];
{ene, vecPsi} = {ene[[order]],vecPsi[[order]]};

Let us check as an exercise the correctness of the result. vecPsi[[i]] is the i-th eigenvector we check the procedure by computing the norm of
H ψ-E ψ for all our energies.

Max[Table[Norm[ Ham.vecPsi[[k]]-ene[[k]]vecPsi[[k]]],{k,1,Length[ene]}]]

Comments

Memory usage

It is important that the reader appreciate the difference between a structure like a sparse array and a usual matrix. Our Ham contains, in full form (10^6) entries while the numbers of effectvely stored values in the sparse version is about 3·(10^3) , as is confirmed by the instruction

Ham

The first number in the output is the number of elements different from zero and stored in Ham. The memory occupation is vastly different, as can be seen using ByteCount:

{ByteCount[Ham],ByteCount[Normal[Ham]]}

These numbers affect all linear operations. As an example for an operation like M.v  in dimension N, roughly speacking with sparse structures the number of elementary operations (sum, products etc.)  grows like N, while for full matrices like (N^2).

The form of the solution

vecPsi  is a list of eigenvectors ( a matrix whose rows are eigenvectors).  This form is very convenient for a number of algebraic manipulation on the whole set of eigenvectors. For example if we denote by matpsi e set of two eigenvectors

(matpsi = Array[psi,{2,3}])//MatrixForm

first eigenvector:

matpsi[[1]]

Suppose you want to multiply each vector by some coefficient c[[i]], the usual multiplication does the work:

(vecCoeff = Array[coeff,2])

The result is obtained simply by multiplication:

vecCoeff*matpsi

If you want to project a given vector on each element of the matrix (e.g. to find the expansion coefficient in the case of a Hilbert space) the usual dot multiplication works:

vec = {v[1],v[2],v[3]}

matpsi.vec

Clear[matpsi,vec,vecCoeff,psi];

Normalization and interpolation

In our case eigenfunctions are given on a finite set of points on the selected interval. Depending on your needs you can use the form produced by Eigensystem, or convert the numbers in interpolating functions etc. The second option is surely the most easy to handle in Mathematica so let us start with it. Our goal is to compute the mean value of (p^2) on the eigenstates.

Interpolation

To use interpolated functions it is convenient to have to whole interval, i.e. to add the two border points to x, x∈ [-L/2,+L/2]

xExtended = Join[{xP[[1]]-h},xP,{xP[[nP]]+h}];

We have also to add the limiting value ( 0 ) to eigenvalues. You can use the same function, Join, or PadRight.

The general form PadRight[vec, newdim, pad, left] extend the vector vec to have dimension newdim by adding the element pad (0 in our case).

These padded values is usually added to the right, but if the (optional) parameter left is specified, a number left of pad is put on the left.

We can now obtain in a single line the whole set of interpolating functions

ψ=Table[Interpolation[Transpose[{xExtended,PadRight[vecPsi[[k]],nP+2,0,1]}]],{k,1,Length[vecPsi]}];

We can use ψ for graphs, for derivatives etc. ψ[[k]] is the (unnormalized) wave function of the k-th level.

As an example

Plot[Re[ψ[[3]][x]],{x,-L/2,L/2},PlotRange-> All]

Mathematica note

Example for PadRight:

PadRight[{v1,v2,v3},5,0,1]

Interpolation:

f = Interpolation[{{x1, f1},{x2,f2},....{xN,fN}}] 
produces an interpolating function such that f[xN] = fN. The only thing to remember is that we must pass to the function couples of the form {xi , value}, while one usually dispose of two vectors, xValues, fValues. 
The couples are easily obtained by taking the Transpose of the list {xValues, fValues}.

Memory usage

With this relatively low numbers of eigenfunctions we have no wasted too much memory:

ByteCount/@{ψ,vecPsi}

Normalization of eigenvectors

We can easily compute norms for our ψ

Short[normψ = Table[NIntegrate[((Abs[ψ[[k]][x]])^2),{x,-L/2,L/2},AccuracyGoal-> 6],{k,1,Length[ψ]}]]

We see that all the norms are equal and are just given by  h

h

In effect Eigensystem produce normalized eigenvectors, let us call (Subscript[c, k]) the components of one of them. The norm (squared) on a grid of step h is (approximatively)

∑ dx ((ψ[(Subscript[x, k])])^2)= ∑ h (Subscript[c, k]^2) = h

This also suggest that with our present precision good results could also been obtained by applying finite approximations directly to vecPsi. This is the first option indicated at the beginning of the section and is left as an exercise to the reader to develop this strategy. We quote only one result which can be useful. An integral on an interval divided with steps of length h is approximated by the Trapezoidal formula

(∫_(Subscript[x, 0])^(Subscript[x, m]))f[x]𝕕x = h[((Subscript[f, 0])/2)+ (Subscript[f, 1])+(Subscript[f, 2])+ . . . (Subscript[f, m-1])+ ((Subscript[f, m])/2)] + 𝒪((h^3))

In our case each function (except derivatives) are zero on the boundary and in the previuous formula the first and last terms vanish. So with good approximation the norms can be in effect computed as finite sum and for eigenstates their values is (Sqrt[h]).

Mean values

Let us consider for example the computation of the mean value of (p^2). We have, for the first 10 levels

meanP2 = Table[(1/(normψ[[k]]  ))*NIntegrate[((Abs[(∂_x) ψ[[k]][x]])^2),{x,-L/2,L/2},AccuracyGoal-> 6],{k,1,10}]

For an harmonic oscillator of unit frequency the expected values are 2 (<(p^2)/2>) = 2 (Subscript[E, n])/2 = (Subscript[E, n])

Table [( n+1./2),{n,0,9}]

Check for waves functions

The known exact result for eigenfunctions is  (remember here n = 0, 1, ... )

psiHarmonic[n_,x_]:= (1/(π^(1/4)))(1./(Sqrt[(2^n)(n)!]))HermiteH[n,x]Exp[-x^2/2];

We also know that on the level with quantum number n  we have <(x^2)>  = (n+1/2), then for large n the wave function must extend at least up to (Sqrt[n]), in a bounded inteval we expect higher states are more affected by these finite size effects.

To compare the analytic and numeric result we have to fix a phase, i.e. the sign of ψ. let us choose for example at a point near to 0  of the grid, the relative signe is

ptcheck = Floor[nP/2]+10;
sign= Table[Sign[psiHarmonic[k-1,xP[[ptcheck]]]/vecPsi[[k]][[ptcheck]]],{k,1,Length[vecPsi]}];

The error in the computation can be estimated by the approximation in (𝕃^2)sense in the interval, i.e.

Table[NIntegrate[ (((psiHarmonic[k-1,x] - sign[[k]] (1/(Sqrt[normψ[[k]]]))ψ[[k]][x]))^2),{x,-L,L},AccuracyGoal-> 7],{k,1,Length[vecPsi]}]

For the first 20 eigenstates the approximation is good, for the rest very bad, as expected. Locally, i.e. for each point of the grid, the first 20 states differ from exact ones for at most 2 / 1000, usually much less:

Table[Max[Abs[vecPsi[[k]]*sign[[k]]/(Sqrt[normψ[[k]]])
 - psiHarmonic[k-1,xP]]],{k,1,20}]

As expected the error grows with principal quantum number.

Generic case

If one is interested in the study of the approximation it is better to write a function which embodies all the steps of the computation. Input parameters will be the length of the box (assumed symmetric) , the number of grid points, the number of eigenvalues and the potential. The output will be {x, ene, ψ}, i.e. the grid points, the energy eigenvalues and the normalized eigenfunctions (we approximate the norm with Sqrt[h] ). The shift value is choosen near the minimum of the potential, in order to get the lowest eigenvalues.

solveSch[V_,numEig_,{nP_,L_}]:=
Module[{xL,xR,h,indDiag,D2,Ham,matV,ene,ψ,norm,xP,minV},
xL = -L/2; xR = L/2;
h = N[(xR - xL)/(nP+1)]; 
xP= xL + h* Range[nP];
D2 = SparseArray[{{i_,i_} ->  -2,{i_,j_}/;j-i==1-> 1,
{i_,j_}/; i-j== 1->  1},{nP,nP}];
indDiag = Table[{i,i},{i,1,nP}];
     matV = SparseArray[indDiag-> V[xP]];
minV = Min[V[xP]];
Ham =- 1/(2 (h^2)) D2 + matV;
{ene,ψ}=Eigensystem[Ham,numEig, Method-> {"Arnoldi","Shift"->( minV-0.01)}];
order = Ordering[ene];
ene = ene[[order]];
ψ = ψ[[order]]/Sqrt[h];
{xP,ene,ψ}]

For example for an anharmonic oscillator with coupling constant g = 10, the solution for the first 5 states on a grid with L = 10 and nP = 1001 is:

V[x_]:= x^2/2 + 10*x^4/2;

{x,ene,ψ} = solveSch[V,5,{1001,10}];

ene

Exercises

Compute the lowest eigenvalue f the Hamiltonian in the quartic potential V[x] = x^2/2 + g x^4/2 (anharmonic oscillator) with g = 10.. The approximate result is:

ℰ = 1.224587036059193459134396953093865213110138999308673549845104365885027406655096801688754290649704430

Compute the lowest eigenvalue for the Schrödinger equation in the potential V[x] = - U for |x| < a. Choose values for a and U.

Compute the lowest eigenvalue for the Schrödinger equation in the potential V[x] = - U/Cosh[x]^2. The exact value is [1]

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

The Schrödinger equation for the radial part of hydrogen atom eigenstates has the form

-(1/2) u''[r] + (- (1/r)+ ((L(L+1))/(2 (r^2))))u[r] = E u[r]; u[0] = 0;

L (the angular momentum) is an integer. Apply the domain truncation method to this equation.

References

1) L.Landau, E. Lifchitz, Quantum Mechanics: Non Relativistic Theory, Butterworth-Heinemann; 3 edition (1981).