anharmonic oscillator - 1

Quantum Mechanics
 K. Konishi, G. Paffuti

Solve numerically the anharmonic oscillator problem using the Hermite polynomials. Study the double well. Try to compute the eigenvalues with 50-digit precision.

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

Introduction

In this notebook the units: ℏ = m = 1  are used.

In this notebook a more detailed study of a prototype model, the anharmonic oscillator, is given. We will use a variational method based on the expansion in oscillator-like wave functions, i.e. exponentials times Hermite polynomials, as explained in notebook [9]. The basis functions are equipped with a scale factor α which will be used as a variational parameter.  For α = 1:

(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,...

Set

(Subscript[F, k])[r] = (Subscript[f, k])[α r]

The factor α is reabsorbed in the scale transformation

r = x/α;        x = α r

The wave functions in the original non-rescaled coordinates are

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

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

The effective Hamiltonian

By a possible change of units the most general Hamiltonian for an anharmonic oscillator can be cast in the form

H = -(1/2)((((∂_r)))^2) + (μ/2)(r^2) + (g/2)(r^4)

In the rescaled x variable :

H = (1/2)(α^2)(p^2) + (μ/(2(α^2)))(x^2) + (g/(2(α^4)))(x^4)

The parameter μ can be reduced to ± 1, the minus sign describes the "double well" case and will be studied at the end of the notebook.  g is the coupling constant.

The scale parameter

The scale parameter α should be fixed, in principle, for each state separately, by minimizing each eigenvalue of Hamiltonian (). We will found that there is a range of α in which a relevant part of the eigenvalues remain stable. This region is larger for low-lying eigenvalues and smaller for highly excited states. The most natural thing is to fix α to the value for which the highest eigenvalue is stable.

This value can be estimated in the following way, which seems to be valid for any kind of potential which grows at infinity, or at least for potentials whose matrix representation is quasi -diagonal. One can estimate the best α value as the one for which the "harmonic oscillator" underlying the basis functions () best fits the (unperturbed) energy of the highest state. Functions () are simply the eigenstates of an oscillator of unit frequency and from known matrix elements for p and x it is easily shown that

ℰ(α)≡⟨n|H|n⟩ = ((2n+1)/4)((α^2)+(μ/(α^2))) + (3/8)g (1/(α^4))(2(n^2)+2n+1)

If we minimize with respect to α we find the equation:

((∂ℰ)/(∂α))=0   ⇒  (2n+1)(α - (μ/(α^3)))- (3/2)(g/(α^5))(2(n^2)+2n+1)

The positive solution of this equation will be taken as the overall scale parameter. In fact we will use also bases so large that any "reasonable" α must be a good choice.

We will use only half of the states, the even or odd states for even or odd eigenvalues. If the basis as dimension dim, the highest principal quantum number is (2 dim). We will see that only about 1/2 of total eigenvalues of the finite matrix are significant, so our effective highest state will be n = dim, and this is the value which will be used in (). An interesting semiclassical check of our choice is found in the problems.

The numbering of states

An annoying point of confusion can be the correspondence between the index of our vectors, which start always from 1, and the principal quantum number n, which conventionally start from zero. Another source of complications comes from the choice of the sequence of even or odd states, in this case as a vector-like index increase by one the corresponding principal quantum number increase by 2.

We will adopt the following convention:  we introduce a variable parity with value 0 or 1 to select even or odd states respectively. The reader can easily convince himself that as k runs in the sequence 1,2,3 ... we have

n = 2 k - 2 + parity;      k = (n+2-parity)/2 .

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 recurrence 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.

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[{p2,x2,x4}];];
x1=SparseArray[{{i_,k_}/;k-i==1->(Sqrt[(n+1)/2])},{dim,dim}];
x1 = (x1+Transpose[x1]);
{p2,x2,x4,x1}]];

Ths function give the solution of equation ()

getscale[n_,g_,μ_]:= x/.FindRoot[(2n+1)(x-(μ/(x^3)))-((3g)/(x^5))(2(n^2)+2n+1),{x,(((3g))^(1/3))}]

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

Anharmonic oscillator, μ > 0

The model

μ = 1;  g=1;

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

Plot[V[x],{x,-1.5(Sqrt[(μ/g)]),1.5(Sqrt[(μ/g)])},AxesLabel-> {"x","V[x]"}]

Example

Here we sketch the procedure for the computation of eigenvalues.

parity = 0;

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

α = 1.0;
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;

eigH = Sort[Eigenvalues[Ham]][[1]]

We will adopt this scheme at first to look at dependence of eigenvalues on dimension, at fixed scale, then at dependence on the scale.

Dependence on the dimensions

μ= 1; parity = 0;
g=1.0;

Let us look at states

state = {1,11,20,50};

i.e. levels :

levels = 2 state - 2 + parity

The dimension will be varied as :

dimRange =Join[Range[10,30,10], Range[40,300,20]];

We fix quite arbitrarily α = 1.

α =1.;

ene = Table[
state0 = Select[state,#≤ dim/2&];
{p2,x2,x4} =N/@ hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
eigH=Sort[Eigenvalues[Ham,-dim/2,Method-> Arnoldi]][[state0]];
PadRight[eigH,Length[state],"*"]
,{dim,dimRange}];

title = {" k=1, n=0 ", " k=11, n=20"," k=20, n=38"," k=50, n=98" };

ScientificForm[TableForm[ene,
TableSpacing-> {3,3},
TableHeadings-> {ToString/@dimRange,title}],{13,12}]

The stability is quite impressive, but we can do even better for higher states, choosing more correctly α. First let us give a look at α dependence at fixed dimension.

Dependence on the scale

dim = 200;
alphaRange = Range[0.5,5,0.5];

ene = Table[
{p2,x2,x4} =N/@ hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
eigH=Sort[Eigenvalues[Normal[Ham]]]
,{α,alphaRange}];

ListPlot[Table[Transpose[{alphaRange,ene[[All,j]]/(Sqrt[j])}],{j,1,dim}],Joined-> True,
PlotRange-> {{0,5},{0,300}}]

g=1.;

α=getscale[300,g,μ]

ene = Table[
state0 = Select[state,#≤ dim/2&];
{p2,x2,x4} = N/@hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
eigH=Sort[Eigenvalues[Ham,-dim/2,Method-> Arnoldi]][[state0]];
PadRight[eigH,Length[state],"*"]
,{dim,dimRange}];

ScientificForm[TableForm[ene,
TableSpacing-> {3,3},
TableHeadings-> {ToString/@dimRange,title}],{13,12}]

We record the last ground state eigenvalue for a later comparison

eneGroundState=Last[ene][[1]]

Error estimate

Energy eigenvalues are so stable that an estimate of the error with Temple method gives a too optimistic prediction. Let us explain this point. Here we compute eigenvalues and eigenvectors for a basis of 200 elements

dim = 200;
α=getscale[dim,g,μ];
{p2,x2,x4} = N/@hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
{eigH,eigV}=Eigensystem[Ham,-dim/2,Method-> Arnoldi];
order = Ordering[eigH];
eigH = eigH[[order]];
eigV= eigV[[order]];

The error estimate in Temple approximation on k-th level is given by

(Subscript[δETemple, k]) = ((Subscript[σ, k])/((Subscript[E, k+1])-(Subscript[E, k])));    (Subscript[σ, k]) = ⟨(Subscript[ψ, k])|(H^2)|(Subscript[ψ, k])⟩ - (((⟨(Subscript[ψ, k])|H|(Subscript[ψ, k])⟩))^2)

where (Subscript[ψ, k]) is our trial normalized eigenstate. To compute matrix elements of (H^2) it is sufficient to enlarge the basis, and this is an exact procedure for anharmonic potential as operators at most connect states differing by 4 units in principal quantum number. The original eigenstate (Subscript[ψ, k]) is padded with 0 to the right and becomes and approximate eigenvector in the enlarged Hilbert space.

Now the quadratic fluctuation σ on normalized states is identically equal to

(Subscript[σ, k])= ⟨(Subscript[ψ, k])|(((H - (Subscript[m, k])))^2)|(Subscript[ψ, k])⟩;   (Subscript[m, k])= ⟨(Subscript[ψ, k])|H|(Subscript[ψ, k])⟩

and this form is less sensible to rounding errors. Let us enlarge the basis and compute δE with the two different prescriptions:

dim2= dim+4;

Here we consider our eigenvectors as vectors in the larger vector space :

eigV=Table[PadRight[eigV[[j]],dim2],{j,1,Length[eigV]}];

We compute the new form of the Hamiltonian

{p2,x2,x4} = N/@hermiteMatrices[dim2,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;

We estimate  errors using quadratic fluctuations . Here is the first formulation of σ, eq.()

δET1=Table[vPsi = eigV[[i]]/Norm[eigV[[i]]];
mean = vPsi.(Ham.vPsi); vT=Ham.vPsi;
(vT.vT- mean^2)/Abs[eigH[[i]]-eigH[[i+1]]],{i,1,Length[eigV]-1}];

The vectors were indeed already normalized, we divided by the norm just to correct for rounding errors.

Here it is the second definition :

δET2=Table[vPsi = eigV[[i]]/Norm[eigV[[i]]];
mean = vPsi.(Ham.vPsi);
vT = Ham.vPsi - mean*vPsi;
vT.vT/Abs[eigH[[i]]-eigH[[i+1]]],{i,1,Length[eigV]-1}];

GraphicsRow[{ListPlot[δET1,PlotLabel-> "δET1"],
ListPlot[δET2,PlotLabel-> "δET2"]}]

We see that the second one, the correct one from the numerical point of view, gives a really small error, too small to be reliable within our numerical precision. In fact for a tridiagonal matrix the Temple definition, or similar ones, test the error induced by the truncation of Hilbert space, and this is in fact what we are asking for, but it is not sensible to errors induced by numerical approximation, as essentially it test the relevance of adding few vectors to our space, if the rest of the computation, the first one, has some numerical induced incertitude, this is not reflected in this estimate. The first graph, δET1, is more reasonable but it suffers from numerical incertitude coming from the incorrect (numerically) way of computing  σ, and this error has almost nothing to do with errors on the eigenvalues.

For these reasons we think it is better to rely on a simpler check: just test stability against the dimension and estimate errors from the variations in eigenvalues.

Let us take the first 100 eigenvalues, as we use Arnoldi method in Eigenvalues function we start the computation from dim = 200, and compare, at fixed α, eigenvalues computed from dim=200 and dim = 250:

dim = 200;
α=getscale[dim,g,μ];
{p2,x2,x4} =N/@ hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
ene1 =Sort[Eigenvalues[Ham,-100,Method-> Arnoldi]];

dim = 250;
{p2,x2,x4} =N/@ hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
ene2 =Sort[Eigenvalues[Ham,-100,Method-> Arnoldi]];

δE = ene2-ene1;

This is the plot  for the difference: we see an error of at most 6×(10^(-12)) which is very reasonable, and is not so different from δET1, showing that rounding errors produce the same "noise" in both computations. For smaller eigenvalues the error is an order of magnitude smaller.

GraphicsRow[{ListPlot[δE],ListPlot[δE[[1;;20]]]}]

If we want a smaller error we have to switch to a better arithmetic.

errorGroundState=δE[[1]]

Anharmonic oscillator, μ < 0

The model

μ = -1;  g=0.1;

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

Plot[V[x],{x,-1.5(Sqrt[((-μ)/g)]),1.5(Sqrt[((-μ)/g)])},AxesLabel-> {"x","V[x]"},PlotRange-> {{-5,5},{-1.5,2}}]

The system is similar to previous one for high excited states, but low lying eigenvalues can be almost degenerate, their splitting being due to tunneling effect: it is the well known "double well potential". We will see in the study of semiclassical approximation that this splitting is of the order of Exp[-1/g], this clearly points to a potential problem for numerical work: for small enough g the splitting can be less than our typical precision, (10^(-13)), this means that in this case some more sophisticated approach has to be taken. We will study this problem in the chapter devoted to WKB, her we present a more standard case, with moderate splitting and g not too small, i.e. we take g=0.1. The main aim of this computation is to show how to produce eigenfunctions.

Eigenvalues and eigenfunctions

g=0.1; μ = -1;

dim=100;

We will compute separately even and odd states

α = getscale[dim, g, μ];

N.B. If we want to use Arnoldi method we have to select correctly the eigenvalues. The instruction

Eigensystem[Ham, - n , Method -> Arnoldi]

select the smallest n eigenvalues, ordered with their absolute value. Our potential is not positive definite, as in the previous case, so we can:

We adopt the second alternative which is simpler and safer. The potential has a (double) minimum at

(x^2)= ((- μ)/(2g));   (Subscript[V, min])= (1/2)μ (((- μ)/(2g))) + (1/2)g (((((- μ)/(2g))))^2)=- (( (μ^2))/(8 g))

We shift the Hamiltonian with (the negative of) this constant value and subtract the same quantity from results.

Vshift = (( (μ^2))/(8g));

parity = 0;
{p2,x2,x4} = N/@hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4 + Vshift IdentityMatrix[dim];

{eigH,eigV} = Eigensystem[Ham,-dim/2,Method-> Arnoldi];
order = Ordering[eigH];
eigHeven=eigH[[order]]-Vshift;
eigVeven = eigV[[order]];

We do the same for odd states

parity = 1;
{p2,x2,x4} = N/@hermiteMatrices[dim,parity];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4 + Vshift IdentityMatrix[dim];

{eigH,eigV} = Eigensystem[Ham,-dim/2,Method-> Arnoldi];
order = Ordering[eigH];
eigHodd=eigH[[order]]-Vshift;
eigVodd = eigV[[order]];

If we plot the first few difference of eigenvalues we see that the first two are almost degenerate

ListPlot[(eigHodd-eigHeven)[[1;;10]]]

To plot eigenfunctions it is faster to consider a grid and use ListPlot.

rgrid = Range[-5,5,0.05];

On this grid we compute eigenfunctions using the base of fixed parity functions:

fB1 = getPsi[eigVeven[[1]],α rgrid,0];

Show[{ Plot[V[x],{x,-1.5(Sqrt[((-μ)/g)]),1.5(Sqrt[((-μ)/g)])},AxesLabel-> {"x","V[x]"},PlotRange-> {{-5,5},{-1.5,2}},PlotStyle-> {Red,Dashing[{0.02}]}],
ListPlot[Transpose[{rgrid,-4fB1}],Joined-> True]}]

We see the symmetric nature of ground state. A similar plot can be done for excited state.

fB1 = getPsi[eigVodd[[1]],α rgrid,1];

Show[{ Plot[V[x],{x,-1.5(Sqrt[((-μ)/g)]),1.5(Sqrt[((-μ)/g)])},AxesLabel-> {"x","V[x]"},PlotRange-> {{-5,5},{-1.5,2}},PlotStyle-> {Red,Dashing[{0.02}]}],
ListPlot[Transpose[{rgrid,2fB1}],Joined-> True]}]

High precision eigenvalues

There are several reasons for which is worth to compute at high precision eigenvalues for anharmonic oscillator. These can offer a test for resummation of perturbation theory, a test for weak coupling semiclassical expansion etc.

It is not difficult to do previous computation at arbitrary precision. We have to use input values of the parameter at high precision and use function Eigenvalue without the option Arnoldi, which works only with machine precision numbers.

Different methods exist, as to compute directly the secular determinant and finding its roots with FindRoot.

Let us start by defining an utility which compute eigenvalues with a prescribed precision, the notations must be obvious at this point of the notebook.

computeEneHighPrecision[state_,parity_,α_,dim_,prec_:MachinePrecision]:= 
Module[{p2,x2,x4,ene0},
{p2,x2,x4}=N[hermiteMatrices[dim,parity],prec];
Ham = ((α^2)/2)p2 + (μ/(2(α^2)))x2 + (g/(2 (α^4)))x4;
ene0=Sort[Eigenvalues[Normal[Ham]]][[state]]];

If our sought precision is precisionGoal, it is reasonable to work with higher and higher dimensions up the precision sought is obtained. The following function simply make a loop, with a while control which stop execution when the precision is obtained or when the dimension go beyond the limit dimMax. An optional boolean argument printResults print intermediate results if settled to True (the default)

energyEstimate[state_,parity_,α_,precGoal_,dimMax_,printResults_:True]:=
Module[{dim,prec,ene,deltaE,enenew,precEne,enesafe,risOld ,deltaEold,precEneold},
dim = 10;
shiftdim = 10;
prec = precGoal+10;
precEne = prec;
ene = computeEneHighPrecision[state,parity,α,dim,prec];
enesafe = ene;
deltaE = 1;
While[ 
dim ≤ dimMax -shiftdim&& precEne ≥ precGoal && deltaE ≥ 10^(-precGoal),
dim = dim+shiftdim;
enenew = computeEneHighPrecision[state,parity,α,dim,prec];
deltaEold =  Norm[ene - enesafe];
precEneold = Precision[ene];
deltaE = Norm[ene - enenew];
precEne = Precision[enenew];
enesafe = ene;
ene = enenew;
risOld = {enesafe,dim-shiftdim,precEneold,deltaEold};
If[printResults== True,
Print[TableForm[
{{"dim=  " <> ToString[dim],
NumberForm[enenew,Ceiling[precGoal]+1],
" prec.= "<>ToString[precEne],"δE:",ScientificForm[deltaE,4]
}}]];];];
If[dim ≤ dimMax && deltaE < deltaEold,{ene,dim,precEne,deltaE},risOld]
]

Let us compute for example the ground-state energy for the anharmonic oscillator with g=1 and μ = 1. The precision goal is 50.

state=1;
precGoal=50;
maxDim = 120;

Here we define our parameter and fix their precision. α is taken fixed as eigenvalues depend very little on it and in any case for high dimensions it is

expected that eigenvalues are α-independent.

dim=200; μ = 1;
g=1; parity = 0;

α = getscale[dim,g,μ];

precParameters=100;
{g,μ,α} = SetPrecision[#,precParameters]&/@{g,μ,α};

({ene,lastDim,precEne,deltaE}=energyEstimate[state,parity,α,precGoal,maxDim])//Timing

Some interesting lessons can be drawn from this computation

{ene-eneGroundState, errorGroundState}

These facts confirm the correctness of our procedure.

As a last test let us make something not completely trivial. For the double-well potential the ground state and first excited state are almost degenerate for small coupling. WKB perturbation theory will give for the energy splitting

ΔE = (1/(Sqrt[π]))(2^(9/4))((e^(-((Sqrt[2])/(3g))))/(Sqrt[g])) ( 1 + 𝒪(g) )

We will compute this expression in the chapter on WKB. The corrections are, up to order (g^2):

dEpert[g_]:= (1/(Sqrt[π]))(2^(9/4))((E^(-((Sqrt[2])/(3g))))/(Sqrt[g]))(1-((71 g)/(24 (Sqrt[2])))-((6299 (g^2))/2304))

For instance

dEpert[0.005]

This rules out any possibility of making a standard machine precision computation of the eigenvalues to obtain such a result. In high-precision computation

we can do the calculation without particular problems, except some extra CPU time.

The ground state

state=1;
precGoal=50;
maxDim = 120;

g=5/1000;

μ=-1;parity=0;
α = getscale[dim/2,g,μ];
{g,μ,α} = SetPrecision[#,precParameters]&/@{g,μ,α};

({ene0,lastDim,precEne,deltaE}=energyEstimate[state,parity,α,precGoal,maxDim])//Timing

The excited state

μ=-1;parity=1;
α = getscale[dim/2,g,μ];
{g,μ,α} = SetPrecision[#,precParameters]&/@{g,μ,α};

({ene1,lastDim,precEne,deltaE}=energyEstimate[state,parity,α,precGoal,maxDim])//Timing

The energy splitting and the comparison with WKB result :

ScientificForm[{ene1-ene0,N[dEpert[g]]},6]

The agreement is satisfactory.