Elementary variational problems

Quantum Mechanics
 K. Konishi, G. Paffuti

Use the Mathematica function NVariationalBound to study the variational method in systems like harmonic oscillator, anharmonic oscillator, the hydrogen atom and the helium atom.

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

The method

Schrödinger equation can be obtained by the variational principle

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

i.e. by optimizing the energy functional on normalizable states. The minimum of the functional gives energy for the ground state of the system.

Mathematica provides a general function, NVariationalBound which solve this problem. The version VariationalBound can also give analytical solution when possible. These functions have been loaded with the package VariationalMethods at the beginning of this notebook.

The general command is

NVariational[{f, g}, {x, xmin, xmax}, funct, parameters]

This function minimize the ratio

((∫f dx)/(∫g dx))

xmin and xmax are extrema of the integration, funct is the trial function, which depends on some parameters, listed in the commands with or without bounds for searching minimum (the default are -∞, +∞). The reader is referred to Mathematica help for more precise and exhaustive information, the examples which we give below could probably suffice to get the reader acquainted with the command.

Harmonic oscillator

Let us start with a problem with a known solution, the harmonic oscillator. We will use natural units, ℏ = m = 1.

NVariationalBound[{(1/2)((u'[x])^2)+(1/2) (x^2) ((u[x])^2),((u[x])^2)},u[x],{x,-∞,∞},(E^(-a (x^2))) (1+b (x^2)),{a,1.0},{b,0.1}]

We see that the value of energy is correct while there is a small error in coefficients. As stressed in the text, in variational methods eigenvalues are more reliable than eigenfunctions. We see also a warning message issued by FindMinimum, which is internally used by the function to minimize the functional. To avoid the message we limit the precision with  AccuracyGoal -> 4 (or 6) (in any case we do not pretend to have a better precision than this one).

In the first list of the input there is the energy and the norm, respectively:

{(1/2)((u'[x])^2)+(1/2) (x^2) ((u[x])^2),((u[x])^2)}

The second argument is the unknown function, u[x], then there is the range of independent variable, x. The crucial point is the next input, which states the kind of functions used to minimize the functional. Here we have used a function with two parameters, a,b:

(e^(-a (x^2))) (1+b (x^2))

The last arguments are the name of the parameters and the starting value used in the search of minimum.

Let us note that energy can be given both as a quadratic form in first derivatives or as quadratic form with second derivatives

NVariationalBound[{-(1/2)u[x]u''[x]+(1/2) (x^2) ((u[x])^2),((u[x])^2)},u[x],{x,-∞,∞},(E^(-a (x^2))) (1+b (x^2)),{a,1.0},{b,0.1},AccuracyGoal->4]

Hydrogen atom

As a second example let us consider radial equation for hydrogen atom. Here the reader must not forget to put in the functional the right weight for radial coordinates. Energy and norm are written as

∫f dr=(∫_0^∞)𝕕r (r^2)  u[r](-(1/2)(1/(r^2))(∂_r) (r^2)(∂_r)u[r] - (1/r)u[r]);     ∫g dr=(∫_0^∞)𝕕r (r^2)  ((u[r])^2)

then the correct formulation of the problem is

NVariationalBound[{u[r](-(1/2)(1/(r^2))(∂_r) ((r^2)(∂_r)u[r])- (1/r)u[r])(r^2),((u[r])^2) (r^2)},u[r],{r,0,∞},(E^(-a r)) ,{a,1.0},AccuracyGoal->4]

We can try to find excited states also, with a judicious choice of trial function and initial values

VariationalBound[{u[r](-(1/2)(1/(r^2))(∂_r) ((r^2)(∂_r)u[r])- (1/r)u[r])(r^2),((u[r])^2) (r^2)},u[r],{r,0,∞},(E^(-a r)) (b- r),{a},{b}]

This example has been taken from the help of Mathematica and uses the version "analytical" for the search of minimum, i.e. VariationalBound instead of NVariationalBound.

Some simple applications

Anharmonic oscillator

Consider the Hamiltonian

-(1/2)((d^2)/(dx^2)) + (1/2)(x^2)+ (1/2)g (x^4).

This model is studied in some detail in other notebooks of this chapter. For g=10 an approximate value of the energy of the ground state is

ene10 =1.224587036;

Let us try NDVariational on this example

NVariationalBound[{(1/2)((u'[x])^2)+(1/2) ((x^2)+10(x^4)) ((u[x])^2),((u[x])^2)},u[x],{x,-∞,∞},(E^(-a (x^2))) (1+b (x^2)),{a,1.0},{b,0.1}]

We can  improve the approximation using a larger set of parameters

NVariationalBound[{(1/2)((u'[x])^2)+(1/2) ((x^2)+10*(x^4)) ((u[x])^2),((u[x])^2)},u[x],{x,-∞,∞},(E^(-a (x^2))) (1+Sum[b[k](x^(2k)),{k,1,3}]),{a,1.0},Sequence@@Table[{b[j],0.},{j,1,3}],AccuracyGoal-> 6]

The result is correct up to the fourth digit, but the execution time has grown.  The reader is referred to Mathematica help for the instruction Sequence, otherwise he can write explicitly {b[1],0} etc. in the previous instruction.

Helium

Mathematica function  VariationalBound allows a one-instruction solution for the helium approximate calculation exposed in the text:

risHe1=Rationalize[Simplify[VariationalBound[{u[r1,r2](-(1/2)(1/(r1^2))(∂_r1) ((r1^2)(∂_r1)u[r1,r2])-(1/2)(1/(r2^2))(∂_r2) ((r2^2)(∂_r2)u[r1,r2])- (Z/r1)u[r1,r2]-(Z/r2)u[r1,r2]+((u[r1,r2])/(Max[r1,r2])))(r1^2) (r2^2),((u[r1,r2])^2) (r1^2) (r2^2)},u[r1,r2],{{r1,0,∞},{r2,0,∞}},(E^(-a (r1+r2))) ,{a}]]]

We have used multipole expansion to write Coulomb interaction between electrons:

(1/(|(Subscript[r, 1])-(Subscript[r, 2])|))= (∑_(ℓ=0)^∞)(1/(Subscript[r, >]))(((((Subscript[r, <])/(Subscript[r, >]))))^ℓ)(Subscript[P, ℓ])[Cos(γ)];

with  (Subscript[P, ℓ]) Legendre polynomials and

(Subscript[r, >])=Max[(Subscript[r, 1]),(Subscript[r, 2])];  (Subscript[r, <])=Min[(Subscript[r, 1]),(Subscript[r, 2])];   (Subscript[r, 1])·(Subscript[r, 2]) = (Subscript[r, 1])(Subscript[r, 2])Cos[γ].

For s - states and for reduced wave functions the formula () implies

J[i,j,k,l]= ∫(Subscript[dr, 1])(Subscript[dr, 2]) (Subscript[f, i])[(Subscript[r, 1])](Subscript[f, j])[(Subscript[r, 2])] (1/(Subscript[r, >]))(Subscript[f, k])[(Subscript[r, 1])](Subscript[f, l])[(Subscript[r, 2])]

as the reader can easily check, and this is what is written in the formula for risHe1.

Double well potential

The following example is not numerically accurate but can give an idea on how a numerical solution for a variational problem can be found.

Let us consider a "double well" potential

V[x] = - (1/2)(x^2) + (g/2) (x^4) = (g/2)((((x^2) - (v^2)))^2)- g ((v^4)/2);      (v^2)= (1/(2g))

The shifted potential has two minima at x = ± v. We try to approximate the system taking as basis functions two gaussian centered on the two classical minima.

Clear[ α,g,x,v,x,f];

f[1,x_]:= Exp[- α (((x-v))^2)/2];
f[2,x_]:= Exp[- α (((x+v))^2)/2];

The two functions are not orthogonal then the mean value of the energy on the subspace generated by them is

ℰ=(((Subscript[c, i])(Subscript[H, ij])(Subscript[c, j]))/((Subscript[c, i])(Subscript[N, ij])(Subscript[c, j])));       (Subscript[N, ij]) = ∫(Subscript[f, i])(Subscript[f, j]) 𝕕x

The c's are the expansion coefficients. The variation with respect c for ℰ gives the eigenvalues equations

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

The eigenvalues will depend on the parameter α, and a minimization with respect to this parameter must be performed.

First we write the generic form of matrix elements

matH = ConstantArray[0,{2,2}];
matN=matH;

matH[[1,2]]= Integrate[(1/2)((∂_x)f[1,x])((∂_x)f[2,x])+ (g/2)((((x^2)-(v^2)))^2)f[1,x]f[2,x],
{x,-∞,+∞},GenerateConditions-> False];
matN[[1,2]]= Integrate[f[1,x]f[2,x],
{x,-∞,+∞},GenerateConditions-> False];
matH = matH+Transpose[matH];
matN = matN+Transpose[matN];

Table[matH[[i,i]]=Integrate[(1/2)((∂_x)f[i,x])((∂_x)f[i,x])+ (g/2)((((x^2)-(v^2)))^2)f[i,x]f[i,x],
{x,-∞,+∞},GenerateConditions-> False];
matN[[i,i]]=Integrate[f[i,x]f[i,x],
{x,-∞,+∞},GenerateConditions-> False];
,{i,1,2}];

matH= Simplify[(matH/.v-> 1/(Sqrt[2g]))];
matN= Simplify[(matN/.v-> 1/(Sqrt[2g]))];

matH//MatrixForm

Here is the Hamiltonian matrix

The eigenvalue problem can be solved by the generalized form of Eigenvalues in Mathematica :  Eigenvalues[{H,N}]. This form accept only machine numbers, so let us define

eig[β_?NumberQ]:=Eigenvalues[{matH/.{g-> 0.1,α->β},matN/.{g-> 0.1,α->β}}]

fEig[β_?NumberQ,k_Integer]:= eig[β][[k]]

The result of fEig will be one of the two eigenvalues.

min1=FindMinimum[fEig[β,1],{β,0.1,2}]

ene1=First[min1];

min2=FindMinimum[fEig[β,2],{β,0.1,2}]

ene2=First[min2];

We have an almost degenerate level

δE = ene1-ene2

The result is not very accurate as the correct one, for this value of g, is δE = 0.056216, as will be shown in the notebooks devoted to anharmonic oscillator, but this example show a general procedure:

This procedure will be amply explored in other notebooks of this chapter.