Uniform WKB approximation

Quantum Mechanics
                      K. Konishi, G. Paffuti

Apply the method of uniform approximation to a quartic oscillator.

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

Introduction

In file [WKBresults.nb] a brief account on uniform WKB approximations has been presented. The main result is that an arbitrary Schrödinger equation

ψ''[x] + (k^2)[x] ψ[x] = 0;

can be uniformly approximated by a model equation

u''[z] + (q^2)[z] u[z] = 0;

as far as q[z] has the same number of classical turning points, i.e. roots of q[z] = 0.

The correspondence is in a mapping x → z tailored in such a way that turning points correspond and the result is

ψ[x] = u[z]/(((z'[x]))^(1/2))

Here we present some simple examples of uniform approximation of bound state wave functions, using linear and quadratic comparison equations.

In the first case the solutions of () are given by Airy functions and the method is a uniformization of usual WKB approximation near the turning points.

In the second example the relevant solutions are the parabolic cylinder functions which, for bound states, simplify to Hermite polynomials and exponentials.

The results will be compared to numerical results obtained with usual variational techniques. The relative functions, collected in the section Functions, are those used in the chapter on variational calculus.

Some standard WKB functions, to compute turning points, action integral etc., are also inserted in that section, the reader is referred to notebooks [NB-11.1], [NB-11.2] for their discussion. The only directed needed is getEneWKB[n,V] where n=0,1,.. is the quantum number of WKB quantization and V is the potential.

Functions

Functions for variational calculations

We report here some functions useful to perform variational computations. For detailed informations on their use see notebooks in the chapter on variational techniques.

The main function is hermiteMatrices which compute in matrix form operators (p^2), (x^2), etc. 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}]];

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:

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 the recursion equation for Hermite polynomials.  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 standard WKB calculations

getTurningPoints[ene_,V_]:= Module[{aL,aR},
aR=x/.FindRoot[ene - V[x],{x,10}];
aL=x/.FindRoot[ene - V[x],{x,-10}];
{aL,aR}];

calcJ[ene_,V_]:= Module[{aR,aL,Jval},
{aL,aR}= getTurningPoints[ene,V];
Jval =(1/π)Re[NIntegrate[(Sqrt[2(ene - V[x])]),{x,aL,aR},AccuracyGoal-> 8]]];

quantize[ene_?NumberQ,V_,n_]:=calcJ[ene,V] - (1/2) -n;

getEneWKB[n_,V_]:= Module[{minV},
minV=First[FindMinimum[V[x],x]];
 y/. FindRoot[quantize[y,V,n] ==0,{y,minV+0.1Abs[minV],minV,200}]]

Uniform approximation with Airy functions

This approximation is not the best suited for studying bound state problems but we will present it for its simplicity and to point out some particular aspects of WKB approximation.

The comparison equation used is

u''[z] - (z - α) u[z] = 0

which corresponds to a linear potential, with forbidden region in z > 0. If a is the corresponding turning point in the original problem, the mapping is given by

(∫_α^z) 𝕕z q[z]= (∫_a^x)k[x] 𝕕x;  ⇒ (2/3) (((z-α))^(3/2))=  (∫_a^x)k[x] 𝕕x = S[a,x]

The general solution of simplified problem is a sum of Ai and Bi functions. If we require an exponential decrease at infinity Ai[z] is the only choice.

In the following we will consider for simplicity symmetric potentials with symetric turning points, ± a.

For the left turning point the comparison equation will be

u''[z] + (z - β) u[z] = 0

We take the freedom in the parameters α and β to simplify the analysis for bound states in a symmetric potential.

To have an overall approximation we need a mapping x → z both for positive and negative x, i.e. we must match the solutions of equations () and (). By simmetry we will use the point x = 0 to impose the matching, z will be positive for x > 0 and negative for x < 0. Taking accont the various signs and phase we can write explicitly for the mapping:

z- α = ((((3/2)(∫_a^x)(Sqrt[2(V-E)]) dx ))^(2/3)); for   x > a;     
z - α = - ((((3/2)(∫_x^a)(Sqrt[2(E-V)]) dx ))^(2/3)); for 0< x < a

By choosing z = 0 for x = 0 we have

α = ((((3/2)(∫_0^a)(Sqrt[2(E-V)]) dx ))^(2/3)) = (((((3π)/4)(1/π)(∫_(-a)^a)(Sqrt[2(E-V)]) dx ))^(2/3))= (((((3π)/4)J ))^(2/3))

where J is the classical integral action for a period: in usual WKB quantization  J = n + 1/2.

The solution in the region x> 0 will be written as

ψ[x] = Ai[z[x]-α]/(((z'[x]))^(1/2))

For x < 0 similarly we will have a solution (Subscript[z, 2]):

(Subscript[z, 2])- β = -((((3/2)(∫_x^(-a))(Sqrt[2(V-E)]) dx ))^(2/3)); for   x < - a;     
(Subscript[z, 2]) - β =  ((((3/2)(∫_(-a)^x)(Sqrt[2(E-V)]) dx ))^(2/3)); for -a< x < 0

By the simmetry of V[x] it follows that β = -α and (Subscript[z, z])[x] = - z[-x]. Then the solution in this region is

ψ[x] = Ai[-(z[x]-β)]/((((Subscript[z, 2])'[x]))^(1/2)) = Ai[(z[-x]-α)]/(((z'[-x]))^(1/2))

This means, as expected, that a prolongation by parity of right hand side solution will provide a solution also for the left hand side, i.e.

(Subscript[ψ, Ai])[x] = Ai[z[x]-α]/(((z'[x]))^(1/2)), x >0; 
 (Subscript[ψ, Ai])[x]= (((-1))^n)(Subscript[ψ, Ai])[-x] for x < 0.

The solution is by construction continuous for n even, and with equal left and right derivatives for n odd. The problem arises: where quantization conditions come in?

Consider odd states: the function must vanish at x = 0, instead we have:

(Subscript[ψ, Ai])[0] = Ai[-α]

Then the "quantization condition" would require -α to be a root of an Airy function. Similarly for even states (Subscript[ψ, Ai]^′)[0] = 0. z[x] is an odd function of x, then z'is even. For x ∼ 0 it behaves like a + b (x^2), i.e. its derivatives vanishes at x=0. The condition ψ'= 0 then implies

Ai'[-α] = 0

i.e. for even states -α must be a root of the derivative of an Airy function. These zeros are known, asymptotically (see for example [1]).

With n even and

ϵ = 3(π/4)(n+(1/2))

the zeros of Ai' are given by the negative of:

(ϵ^(2/3))(1 - (7/48)(1/(ϵ^2))+(35/288)(1/(ϵ^4))-(181223/207360)(1/(ϵ^6))+(18683371/1244160)(1/(ϵ^8))+ ...)

Similarly for n odd the successive (negative of ) the zeros of Ai are given by

(ϵ^(2/3))(1 + (5/48)(1/(ϵ^2))-(5/36)(1/(ϵ^4))+(77125/82944)(1/(ϵ^6))-(108056875/6967296)(1/(ϵ^8))+ ...)

From these results we see that for large n the matching of Airy functions would imply α = ϵ and form () J = n+1/2, the usual WKB quntization. For finite n the two conditions do not match, this means that we have two possibilities: or to adopt Airy condition, and change semiclassical estimation of energy levels, or keep usual WKB condition, and admit a small discontinuity at x = 0. This is another example from which it is apparent how Bohr-Sommerfeld type of quantization conditions holds better and better for high excited states, while for low n they can give wrong results.

Clearly the Airy type of match is exact for linear potentials but, for example, gives the wrong result for harmonic oscillator. A nice example for the linear potential well is in the help of Mathematica (select the word AiryAiPrime and go to the help via menu, ⇑⌘F on a Mac).

Energy levels

The function getEneWKB[n,V], introduced in notebook [], simply solves the quantization condition, in the form J = (n+1/2), where J is the action integral. Let us see what happens if we change the quantization conditions following eq.() and (). From () we have that

(((((3π)/4)J ))^(2/3)) = zero of Ai or Ai'  for odd/even states

We compute the first zeros of these functions (the zeros of Ai are obtained by a built in Mathematica function):

ζAiPrime[n_]:=x/.FindRoot[AiryAiPrime[x],{x,-(3Pi/8(Max[0,4n-1]))^(2/3),-(3Pi/8(4n+3))^(2/3)}];

maxnumQ = 20;
zerosOfAiPrime = -Table[ζAiPrime[k],{k,0,maxnumQ}];
zerosOfAi = -N@AiryAiZero[Range[1,maxnumQ]];

The quantization condition will be

quantizeAi[ene_?NumberQ,V_,n_]:=Module[{s,ζ},
If[Mod[n,2]== 0,s = n/2+1; ζ = zerosOfAiPrime[[s]]; prin=0,
s = (n+1)/2; ζ =zerosOfAi[[s]];prin=1];
(* Print[{prin,s,ζ}]; *)
calcJ[ene,V] -(4/(3π)) (ζ^(3/2))]

and the corresponding energy function :

getEneWKBAi[n_,V_]:= Module[{minV},
minV=First[FindMinimum[V[x],x]];
 y/. FindRoot[quantizeAi[y,V,n] ==0,{y,minV+0.1Abs[minV],minV,200}]]

Here is the example of harmonic oscillator.

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

The usual quantization conditions give exact results

getEneWKB[#1,V]&/@Range[0,10]

The analogous results for this airy - type quantization are

getEneWKBAi[#1,V]&/@Range[0,10]

The same computation can be done for an arbitrary potential. Let us consider an anharmonic oscillator and compute energy levels with the known variational procedure.

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

dim=200;
{p2,x2,x4,x1}=N/@ hermiteMatrices[dim];
Sort[Eigenvalues[(1/2)p2 + (1/2)x2 + (1/2)x4]][[1;;11]]

Off[FindRoot::"lstol"];
getEneWKB[#1,V]&/@Range[0,10]
On[FindRoot::"lstol"];

Off[FindRoot::"lstol"];
getEneWKBAi[#1,V]&/@Range[0,10]
On[FindRoot::"lstol"];

We see that for low levels there is disagreement. In general both conditions give only approximate results, but the Airy - type approximation introduces additional errors, therefore we discard this possibility. This means that we will take in the following energies computed with usual WKB formula, and this implies that our wave functions will have a small discontinuity at the origin (or their derivatives): we will avoid singularities at the turning points but we have to pay this price. This is the main reason which advises to avoid this approximation for bound states

Wave functions

Let us now construct a uniform approximation to wave functions. First we have to compute classical action. It is convenient to compute this quantity by solving a differential equation 𝒮'[x] = |p[x]| . We compute the action following the conventions in eqs.() and (). The base point will be the left turning point and we can limit ourselves to the computations in the region x < 0:

calcAction[ene_,V_]:=  Module[{aL,aR,T,C,σ,u},
{aL,aR} = getTurningPoints[ene,V];
σ =u/.First[NDSolve[{u'[z] ==  
(Sqrt[2Abs[(-V[z]+ene )]]),u[aL]== 0},u,{z,10 aL, 0}]];
{σ,aL,ene}]

From the action it is immediate to compute the mapping x → z, i.e. the function z[x]. We will use an interpolating function:

calcZXAi[numQ_,V_]:= Module[
{aL,eneWKB,xgrid,𝒮,temp,zval,zofx,α,xgrid1,xgrid2},
eneWKB=getEneWKB[numQ,V];
{𝒮,aL,temp} = calcAction[eneWKB,V];
α = (((((3.π)/4)(numQ+(1/2))))^(2/3));
xgrid1 = Range[0.,aL,-0.012];
xgrid2= Range[aL,10aL,-0.012];
xgrid =Reverse[Join[xgrid1,xgrid2]];
zval =Reverse[ -α + Join[(((((((3/2) 𝒮[#1]&/@xgrid1))^2)))^(1/3)),-(((((((3/2) 𝒮[#1]&/@xgrid2))^2)))^(1/3))]];
zofx =Interpolation[Transpose[{xgrid,zval}]];
zofx]

Let us make an example using again an harmonic oscillator

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

We compute a whole set of mappings :

maxNumQ =10;
zofx= calcZXAi[#1,V]&/@Range[0,maxNumQ];

This is an example

Plot[zofx[[1]][x],{x,-5,0}]

The wave function is now immediate :

ψWKBAi[n_,x_]:= Module[{α},α=(((((3.π)/4)(n+(1/2))))^(2/3));

If[x≤ 0,AiryAi[-(zofx[[n+1]][x] + α)]/zofx[[n+1]]'[x],(((-1))^n)ψWKBAi[n,-x]]]

Here is an "artistic" plot to show that there are no singularities at the turning points:

eneWKB=getEneWKB[#1,V]&/@Range[0,maxNumQ];

Plot[Evaluate[Prepend[Table[eneWKB[[n+1]]+ψWKBAi[n, x], {n, 0, maxNumQ}],V[x]]],{x, -5,5},Frame->True,Axes->False,PlotRange-> {Automatic,{0,eneWKB[[maxNumQ+1]]+1}}]

Let us now normalize the wave functions and compare them with exact results. Inside the structure of an interpolating function is contained, as first component, the range of the coordinates, this will allow to use in the normalization only computed values, avoiding annoyng warnings due to extrapolations if the integral limits are outside the interpolation range. Here is an example:

Flatten[First[zofx[[5]]]]

Here is the computation of norms :

calcNorm[n_,ψ_]:= Module[{a,temp,norm2},
a= Flatten[First[zofx[[n+1]]]][[1]];
norm2 = NIntegrate[((ψ[n,x])^2),{x,a,-a}];
1/Sqrt[norm2]]

norms = calcNorm[#1,ψWKBAi]&/@Range[0,maxNumQ]

The known wave functions are

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

We see that errors are of a few percent and, as expected, they are greater around x = 0.

Plot[Evaluate[Table[ψHO[k,x]-(((-1))^k)norms[[k+1]]ψWKBAi[k,x],{k,0,6}]],{x,-5,5},PlotRange-> All]

Quadratic approximation

This approximation is better tailored for two points inversion problems. As explained in notebook [WKBresults.nb] the comparison equation is

u''[z]+ (q^2)[z]=0;      (q^2)[z] = t - (z^2);

The parameter t is a function of energy and is fixed by requiring coincidence of turning points between original problem and ():

(∫_(Subscript[x, 1])^x)k[x]𝕕x = (∫_(-(Sqrt[t]))^z)q[z] 𝕕z

This fixes the value of the parameter t as a function of energy as to match second turning point we must have

(∫_(Subscript[x, 1])^(Subscript[x, 2]))k[x]𝕕x  =  (∫_(-(Sqrt[t]))^(Sqrt[t]))q[z] 𝕕z = t (π/2)

In particular for bound states in the semiclassical approximation the quantization condition gives

(∫_(Subscript[x, 1])^(Subscript[x, 2]))k[x]𝕕x  = π (n+(1/2))

and in this case

t = 2n+1

The general solution of the comparison equation is

u[z] = (Subscript[c, 1]) (Subscript[D, ((t-1)/2)])(-(Sqrt[2]) z)+ (Subscript[c, 2]) (Subscript[D, (1/2) (-t-1)])(i (Sqrt[2]) z)≡ (Subscript[c, 1])(Subscript[u, 1])+ (Subscript[c, 2])(Subscript[u, 2])

Where (Subscript[D, ν])[z]are parabolic cylinder functions. If the wave function must be bounded at z → ∞ , we have (Subscript[c, 2])= 0. For t = 2n+1 D[-(Sqrt[2]) z] is bounded also at z→ - ∞. In this limit the parabolic cylinder functions reduce to Hermite polynomials times an exponential, i.e. as expected, to wave functions of harmonic oscillator:

(Subscript[D, n])[-(Sqrt[2]) z] = (((-1))^n)(2^(-n/2)) (e^(-((z^2)/2))) HermiteH[n,z]

This means that, apart an overall normalization factor, the semiclassical approximation for a bound state is

(Subscript[ψ, n])[x] = (1/(Sqrt[z']))(Subscript[H, n])[z] Exp[-(z^2)/2]

To write explicitly th eapproximation we first compute the action in z valiables, i.e. the integral of q[z]. By rescaling z → (Sqrt[t])y we have, in the allowed and forbiden region respectively:

(Sqrt[t])Integrate[(Sqrt[1-(y^2)]),{y,-1,zt},Assumptions-> -1< zt< 1]

(Sqrt[t])Integrate[(Sqrt[-1+(y^2)]),{y,-1,zt},Assumptions-> zt<-1]

This allow us to write the function  σz for the action:

σz[z_?NumberQ,n_]:= Module[{t,zt},
t = 2n+1;
zt= z/Sqrt[t];
If[zt ≤  -1,(t/2) (zt (Sqrt[-1+(zt^2)])+Log[-zt+(Sqrt[-1+(zt^2)])]),
(t/4) (π+2 zt (Sqrt[1-(zt^2)])+2 ArcSin[zt])]];

We can always use our previous function calcAction, which we repeat for the reader:

calcAction[ene_,V_]:=  Module[{aL,aR,T,C,σ,u},
{aL,aR} = getTurningPoints[ene,V];
σ =u/.First[NDSolve[{u'[z] ==  
(Sqrt[2Abs[(-V[z]+ene )]]),u[aL]== 0},u,{z,10 aL, 0}]];
{σ,aL,ene}]

and write with only few changes (due to different choices of signs in σz ) the function calcZXQuad

calcZXQuad[numQ_,V_]:= Module[
{aL,eneWKB,xgrid,𝒮,temp,zval,zofx},
eneWKB=getEneWKB[numQ,V];
{𝒮,aL,temp} = calcAction[eneWKB,V];
xgrid = Range[0.,10aL,-0.01];
zval= (zvar/.FindRoot[σz[zvar,numQ]== 𝒮[#1],{zvar,0.1aL,0}])&/@xgrid;
zofx =Interpolation[Transpose[{xgrid,zval}]];
zofx]

The wave function is simply :

ψWKBQuad[n1_,x_]:= If[x≤ 0,(1/(Sqrt[zofx[[n1+1]]'[x]]))HermiteH[n1,zofx[[n1+1]][x]]Exp[-0.5((zofx[[n1+1]][x])^2)],
(((-1))^n1)ψWKBQuad[n1,-x]]

For the harmonic oscillator the change of valiables x → z is just a change of scale, so to make a nontrivial check let us consider an anharmonic oscillator

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

We compute a whole set of mappings :

maxNumQ =10;
zofx= calcZXQuad[#1,V]&/@Range[0,maxNumQ];

This is an example

Plot[zofx[[7]][x],{x,-4,0}]

Now let us plot the eigenfunctions

eneWKB=getEneWKB[#1,V]&/@Range[0,maxNumQ];

norms = calcNorm[#1,ψWKBQuad]&/@Range[0,maxNumQ];

Plot[Evaluate[Prepend[Table[eneWKB[[n+1]]+norms[[n+1]]ψWKBQuad[n, x], {n, 0, maxNumQ}],V[x]]],{x, -5,5},Frame->True,Axes->False,PlotRange-> {Automatic,{0,eneWKB[[maxNumQ+1]]+1}}]

There are no discontinuities at turning points, as expected, and no more internal discontinuity, as in the Airy case.

Let us now compute numerical wave functions with the function getPsi and compare them with our WKB result:

dim=200;
{p2,x2,x4,x1} =N/@ hermiteMatrices[dim];
{eigH,eigV}= Eigensystem[(1/2)p2 + (1/2)x2+(1/2)g x4];
order=Ordering[eigH];
eigH=eigH[[order]];
eigV=eigV[[order]];

We take a grid on which we evaluate both numerical ψ and WKB approximations

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

ψNum=getPsi[eigV[[#1]],rgrid]&/@Range[1,maxNumQ+1];

listψWKBQuad=Table[(norms[[num+1]]ψWKBQuad[num,#1]&/@rgrid),{num,0,maxNumQ}];

After fixing the relative signs

signs= Table[Sign[ψNum[[k,30]]*listψWKBQuad[[k,30]]],{k,1,maxNumQ+1}];

we plot the differences :

ListPlot[
Table[Transpose[{rgrid,ψNum[[num]]-signs[[num]]listψWKBQuad[[num]]}],
{num,1,maxNumQ+1}],Joined-> True,PlotRange-> All]

We see that everywhere and for every state the accuracy is around 0.002

References

1) M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions, Dover Publications Inc. (1965).