Chebyshev basis and pseudospectral method

Quantum Mechanics
 K. Konishi, G. Paffuti

Use the  Chebyshev polynomials to implement  a variational calculation. Study the pseudospectral method with this basis.

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

Introduction

Review of the theory

General remarks

In notebook [4] a general overview of collocation and pseudospectral methods has been presented.  Here we specialize the subject to the case of Chebyshev polynomials. Let us repeat some of the ideas introduced in notebook [4] to make the notebook self-contained.

In general a collocation method is defined by the choice of a gridpoints (Subscript[x, a]) and by the discretization of the equation on this grid:

-(1/2)(Subscript[𝒟, ab]^((2)))(Subscript[ψ, b]) + V[(Subscript[x, a])](Subscript[ψ, a]) = E (Subscript[ψ, a])

The derivative matrix 𝒟 has in any finite difference scheme elements different from zero for |a-b| ≤ k, for example in the nearest neighbors approximation of second derivative on a uniform grid

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

we have a tridiagonal matrix, i.e. k=1. The number of points involved in the computations of derivatives can be related to the degree of a polynomial which interpolate the function (in general a local interpolation): the approximation are made so as to be exact for these polynomials. For instance polynomials of second degree generate formula (). A collocation method is defined with a maximum degree polynomial, i.e. for a grid of (Subscript[n, G]) points by a polynomial of degree (Subscript[n, G])-1. In this case the matrix 𝒟 is in general a full matrix and all points concur to define the derivative at each point of the grid.

The best grid from the numerical point of view are related to the grid points used by Gaussian integration formulas, i.e. are zeros of an n-degree polynomial, orthogonal with respect some measure on the given interval. In this case collocation method is also called pseudospectral method. The reason being that equations () are seen to be equivalent to usual spectral equations, as those used in notebooks [6, 9], if matrix elements are computed by Gaussian integration formulas, as we have done in notebooks [7], [9].

In contrast to the spectral methods equation () is defined directly on "real" space, not in the space of coefficients of a spectral expansion. Moreover it is often much faster to solve: for an n-grid points we have to compute only n values of the potential, whereas one must compute of the order of n^2/2 values of the corresponding matrix elements in the spectral method.

In this notebook the interval [-1,1] will be used, and we will show how more usual infinite or semi-infinite domains can be reduced to this case.

The set of polynomials we will use is made of Chebyshev polynomials of the first kind, (Subscript[T, n])[x]. They are simply defined by

(Subscript[T, k])[Cos[θ]] = Cos[k θ].

They are orthogonal with respect the measure (((1 - (x^2)))^(1/2)) :

(∫_(-1)^1) 𝕕x (Sqrt[1-(x^2)]) (Subscript[T, a])[x] (Subscript[T, b])[x] = 0    if a≠ b

The set of polynomials of degree k ≤ n describe a linear space (Subscript[𝒱, n]) of dimension n+1.

The grid : internal points and boundary

Let us take a grid consisting of  n internal points in the interval [-1 , +1].  The values of the wave function on these points will be our unknown quantities. We have to embed the grid and the values in some interpolation scheme which allows the use of derivatives, acting on the interpolated function. We have to choose between two possibilities:

The second procedure is more flexible, as can account for different boundary conditions, and is also more easily programmable, so we follow this second alternative. Our extended grid will contain

(Subscript[n, E])= n + 2

grid points (the two boundary points are ± 1). The number of intervals is

(Subscript[n, I])= n + 1 ≡ N

This kind of grids are usually named the Gauss-Lobatto grids. We will use a Gauss-Lobatto grid based on the Chebyshev polynomials of the first kind.

The N+1 = (Subscript[n, E])values of a function at points of the extended grid can be associated to the set of polynomials of degree ≤ N. In the case of Chebyshev polynomials an interpolating polynomial  of the function can be written as

f[x] = (∑_(i=0)^N)(Subscript[c, i])(Subscript[T, i])[x]

where the N + 1 coefficients are fixed by the N+1 values of the function on the grid.

The Gauss-Lobatto grid is defined as the set of roots of the N+1 degree polynomial

(1-(x^2))(Subscript[T, N]^('))[x]

and it consists of the N+1 points

(Subscript[x, i]) = i (π/N);       i=0,...N

The collocation method

The pseudospectral collocation method consists in discretizing the derivatives on the grid of points with a maximal degree of interpolation, i.e. defining the discretization as to be exact for polynomials of degree N, in our notation. The value of a function on the grid can be most easily expressed through a set of cardinal functions, in our case polynomials, with the property

(Subscript[C, i])[(Subscript[x, a])]= (Subscript[δ, ia])

The unique polynomial which interpolates a function ψ is then trivially given by

ψ = (∑_(b=0)^N)(Subscript[ψ, b]) (Subscript[C, b])[x];    (Subscript[ψ, b])≡ ψ[(Subscript[x, b])]

The cardinal functions can be in general expressed in form of Lagrange interpolating polynomials

(Subscript[ℒ, i])[x]=(∏_(j≠i)^N)(((x-(Subscript[x, 0]))(x-(Subscript[x, 1]))..(x-(Subscript[x, i-1]))(x-(Subscript[x, i+1])).. (x-(Subscript[x, N])))/(((Subscript[x, i])-(Subscript[x, 0]))((Subscript[x, i])-(Subscript[x, 1]))..((Subscript[x, i])-(Subscript[x, i-1]))((Subscript[x, i])-(Subscript[x, i+1])).. ((Subscript[x, i])-(Subscript[x, N]))))

but in the Chebyshev case it is also possible to write a simpler closed form

(Subscript[C, j])[x]= (((-1))^(j+1))(((1-(x^2)))/((Subscript[c, j])(N^2)(x-(Subscript[x, j]))))(((Subscript[dT, N])[x])/dx)= (2/(N (Subscript[p, j])))(∑_(m=0)^N)(1/(Subscript[p, m]))(Subscript[T, m])[(Subscript[x, j])](Subscript[T, m])[x].

with

(Subscript[c, N])=2,     otherwise        (Subscript[c, j])=1;               (Subscript[p, 0])=(Subscript[p, N])=2;   otherwise    (Subscript[p, j])=1 .

The two forms of the function in () are related by trigonometric identities through ().

From () it follows that the derivatives on grid points are easily written in terms of derivatives of cardinal functions

Dψ [x] = (∑_(b=0)^n)(Subscript[C, b]^('))[x](Subscript[ψ, b]);   ⇒ Dψ [(Subscript[x, a])] = (∑_(b=0)^n) (Subscript[DC, b])[(Subscript[x, a])](Subscript[ψ, b])≡ (Subscript[𝒟, ab])(Subscript[ψ, b])

Similarly for any derivative :

(Subscript[𝒟, ab]^((L)))= ( (D^((L)))(Subscript[C, b]))[(Subscript[x, a])]

Derivative matrics 𝒟 can be easily computed, but we will use a Mathematica function which does this. The sequence of instructions is:

grid= x0 +(1/2) L(1.0 -  Cos[π Range[0,N]/(N)]);

fdf = NDSolve`FiniteDifferenceDerivative[ord, grid, DifferenceOrder -> "Pseudospectral"];

smat = fdf[DifferentiationMatrix[]]

We will embed these functions in some more user-friendly interface, adapted to our case.

Once we have 𝒟 matrices the Schrödinger equation on the grid is easily written:

-(( 1)/2)(∑_b)(Subscript[𝒟, ab])(Subscript[φ, b]) + (Subscript[V, a]) (Subscript[φ, a]) = E (Subscript[φ, a]);   (Subscript[V, a])= V[(Subscript[x, a])].

V is the potential term. We will use only the zero boundary conditions for this equation. This means that values at the first and last points of the grid do not appear in (). The output of above Mathematica functions gives the full matrix, taking into account the boundary points, i.e. are (Subscript[n, E])×(Subscript[n, E])=(N+1)×(N+1) matrices.  Only the internal part of these matrices, from row 2 up to row N, and similarly for columns, enter in ().

Integration formulas

Once we solved () we need a tool to compute the norms, scalar products etc. This can be done by using the Gauss integration formulas. For the Gauss-Lobatto grid given above these formulas reads

(∫_(-1)^1)f[ξ] (Sqrt[1-(ξ^2)]) dξ = (∑_(i=1)^n)(Subscript[w, i])(Subscript[f, i]);

n is the number of internal points of the grid. Abscissas and weights are accordingly:

(Subscript[ξ, i])= Cos[i (π/(n+1))];  (Subscript[w, i])= (π/(n+1))((Sin[i (π/(n+1))])^2)= (π/(n+1))(1-(Subscript[ξ, i]^2))

The reader has to remember that N = n + 1, the abscissas are exactly our grid points.

We seldom have an explicit root in front of our functions, in its absence we can write

(∫_(-1)^1)g[ξ]  dξ= (∑_(i=1)^n)(Subscript[w, i])(Subscript[g, i]) (1/(Sqrt[1-(Subscript[ξ, i]^2)])) = (π/(n+1))(∑_(i=1)^n)(Subscript[g, i])(Sqrt[1-(Subscript[ξ, i]^2)])

Intervals

In Quantum Mechanics the usual interval to be considered is [-∞,+∞] for one dimensional problem, or [0,+∞] for radial problems. One strategy to reduce this study to a finite interval equation is to adopt a "truncation domain technique", i.e. to put a maximum for our coordinates and make possibly a rescaling to bring them to the interval [-1,+1]. This technique works, but introduces some spurious errors due to the truncation procedure. It is better, and easier, to change the variables in the Schrödinger equation. We will use an algebraic change of variables, different possibilities can be explored by the reader.

In the formulas below λ is a free (variational) parameter, ξ ∈ [-1,1];   x ∈ [-∞,∞] and r ∈ [0,∞]. A rescaling of the wave function simplifies the Schrödinger equation bringing it to a canonical form in the variable  ξ . We give the results which can be checked by the reader by doing explicitly the change of variable in the derivatives.

Infinite interval

x ∈ [-∞,+∞] ;          ξ ∈ [-1,+1]

x = ((λ ξ)/(Sqrt[1-(ξ^2)]));   ξ = (x/(Sqrt[(λ^2)+(x^2)]));      ψ= (1/(Sqrt[ξ']))φ[ξ] = ((Sqrt[λ])/(((1-(ξ^2)))^(3/4)))φ[ξ];

The Schrödinger equation for φ is:

((-1)/(2(λ^2)))(((1-(ξ^2)))^3)(∂_ξ^2) φ[ξ] - (1/(8(λ^2)))(6-9 (ξ^2)+3 (ξ^4)) φ[ξ] + V[x[ξ]]φ= E φ

The Jacobian from x to ξ variables is

(dx/(dξ)) = (λ/(((1-(ξ^2)))^(3/2)))

Using () the norms are written as (n is the number of internal points)

(∫_(-∞)^∞)𝕕x (ψ^2) = (∫_(-∞)^∞)𝕕ξ  (λ/(((1-(ξ^2)))^(3/2)))(ψ^2) ≃ (π/(n+1))(∑_(i=1)^n)(Subscript[ψ, i]^2) (λ/((1-(Subscript[ξ, i]^2))))

Semi- infinite interval

r ∈ [0,+∞] ;          ξ ∈ [-1,+1]

r = λ (((1+ξ))/((1-ξ)));   ξ = ((r-λ)/(r+λ));    ψ= (1/(Sqrt[ξ']))φ[ξ] = ((Sqrt[2λ])/(|1-ξ|))φ[ξ];

The Schrödinger equation for φ is:

((-1)/(2(λ^2)))(1/4)(((-1+ξ))^4)(∂_ξ^2) φ+ V[r[ξ]]φ= E φ

The Jacobian from r to ξ variables is

(dr/(dξ)) = ((2λ)/(((1-ξ))^2))

Using () the norms are written as (n is the number of internal points)

(∫_(-∞)^∞)𝕕x (ψ^2) = (∫_(-∞)^∞)𝕕ξ  ((2λ)/(((1-ξ))^2))(ψ^2) ≃ (π/(n+1))(∑_(i=1)^n)(Subscript[ψ, i]^2) ((2λ)/(((1-(Subscript[ξ, i])))^2))(Sqrt[1-(Subscript[ξ, i]^2)])

The functions getD2 and normalize will provide automatically derivatives for λ=1 and normalized wave functions, for the infinite and the semi-infinite intervals. These are the only functions that reader must use to compute the eigenvalues and eigenvectors.

The utility scalar is provided to compute scalar products between states, computing in this way also matrix elements for q-operators.

Functions

The function getDCheb compute differentiation matrix and grid for Chebyshev points. Only internal points are present in the output.

getDCheb[x0_,L_,dim_]:= Module[{D1,D2,grid,fdf},
grid= x0 +(1/2) L(1.0 -  Cos[π Range[0,dim+1]/(dim+1)]);
fdf = NDSolve`FiniteDifferenceDerivative[2, grid, DifferenceOrder->"Pseudospectral"];
     D2=Normal[fdf[DifferentiationMatrix[]]];
{grid[[2;;dim+1]],D2[[2;;dim+1,2;;dim+1]]}]

getD2 compute the effective matrix for second derivative, if a second argument different from zero is provided, the result is second derivative in radial coordinates. The output is {ξgrid, xgrid1, D2L1} containing respectively the grid points in ξ variable (in the interval [-1,1]), the grid points in x (or r) variable, with scale parameter  λ = 1, and the matrix of second derivatives, always for λ=1.

If a different λ is needed it is sufficient to put

xgrid = λ xgrid1;          D2 = 1/(λ^2)D2L1

getD2[dim_,rad_:0]:= Module[{matM2,matK,ξgrid,D2,D2X,xgridL1},
{ξgrid,D2}= getDCheb[-1,2,dim];
If[rad== 0,
matM2 = DiagonalMatrix[(((1-((ξgrid)^2)))^3)];
matK = DiagonalMatrix[(3/4)(2- ((ξgrid)^2))(1-((ξgrid)^2))];
D2X = matM2.D2 + matK;
xgridL1= (( ξgrid)/(Sqrt[1-((ξgrid)^2)]));
Return[{ξgrid,xgridL1,D2X}];];
matM2 = DiagonalMatrix[(1/4)(((1-ξgrid))^4)];
D2X = matM2.D2 ;
xgridL1= ((1+ ξgrid)/((1-ξgrid)));
{ξgrid,xgridL1,D2X}];

Normalize transforms a vector of values computed on the ξ grid into a normalized eigenfunction, by multiplication for relevant factors, as indicated in () and (). The input needed is  ξgrid, λ , φ, i.e. the grid in ξ variable, the value of scale parameter and the vector of wave function's values on the ξ grid. A fourth, optional parameter allows the computation for a radial variable.

normalize[ξgrid_,λ_,φ_,rad_:0]:= 
Module[{dim,norm2,ψ},
dim=Length[ξgrid];
If[rad== 0,
ψ = ((Sqrt[λ])/((((1+ξgrid)(1-ξgrid)))^(3/4)))φ;
norm2=(π/(dim+1))Tr[(λ/((1+ξgrid)(1-ξgrid)))(ψ^2)];,
ψ = (Sqrt[2λ/(((1-ξgrid))^2)]) φ;
norm2=(π/(dim+1))Tr[((2λ)/(((1-ξgrid))^2))(ψ^2) (Sqrt[(1-ξgrid)(1+ξgrid)])];];
ψ/(Sqrt[norm2])];

scalar is similar to normalize, it needs as input two wave functions, and as usual ξgrid and λ.

scalar[ξgrid_,λ_,φ_,ψ_,rad_:0]:= Module[{dim,ris},
dim=Length[ξgrid];
If[rad== 0,
ris = (π/(dim+1))Tr[(λ/((1+ξgrid)(1-ξgrid)))ψ φ];,
ris= (π/(dim+1))Tr[((2λ)/(((1-ξgrid))^2))ψ φ (Sqrt[(1-ξgrid)(1+ξgrid)])];];
ris];

Example 1

To show how the machinery works let us take the harmonic oscillator. We do not use any special knowledge of the system, and every part of the computation is valid for an arbitrary one dimensional potential defined on the interval [-∞, +∞].

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

dim=100;

Prototype computation

Here is the computation of grids and second derivative matrix

{ξgrid,xgridL1,D2} = getD2[dim];

Here is the choice of the scale parameter and the computation of eigenvalues.

λ = 2.0;

Sort[Eigenvalues[((- 1)/(2(λ^2)))D2 +DiagonalMatrix[ V[λ xgridL1]]]][[1;;10]]

λ dependence

In general systems it is worthwhile to have an idea on eigenvalues as  functions of the scale parameters λ, in order to look for a stability valley for λ.

rangeL = Range[1,20,1];
dim=100;
numEig=Floor[dim/3];

eneL= Table[Sort[Eigenvalues[((-1)/(2(λ^2)))D2 +DiagonalMatrix[ V[λ xgridL1]]]][[1;;numEig]],
{λ,rangeL}];

ListPlot[Table[Transpose[{rangeL,eneL[[All,j]]}],{j,1,numEig}],Joined-> True,
PlotRange-> {{0,25},{0,numEig+5}},AxesOrigin-> {0,0},AxesLabel-> {"λ","(Subscript[E, n])"}]

A robust stability is present for many eigenvalues. A more detailed computation can be done by choosing a definite scale:

λ=8;

{eigH,eigV}= Eigensystem[((-1)/(2(λ^2)))D2 +DiagonalMatrix[ V[λ xgridL1]]];
order = Ordering[eigH];
eigH = eigH[[order]][[1;;numEig]];
eigV=eigV[[order]][[1;;numEig]];

The eigenvectors consist of function values at the points of ξ grid. Using normalize we can transform them in normalized functions. For instance the ground-state function is:

ψ=normalize[ξgrid,λ,eigV[[1]]];

The theoretical one is :

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

We compare the two, note that even the absolute normalization of the numerical wave function has been fixed numerically.

xgrid = λ xgridL1;

sign= Sign[ψ[[Floor[dim/2]-3]]];

Show[{Plot[ψTh[0,x],{x,-15,15},PerformanceGoal-> "Speed",
PlotRange-> {{-5,5},{0,0.8}},PlotStyle-> Red],
ListPlot[Transpose[{xgrid,sign*ψ}]]}]

Tha maximum discrepancy of the functions on  the grid is :

Norm[ψ*sign- ψTh[0,xgrid],∞]

Mean values

With numerical wave functions we can compute the amplitudes, mean values etc. For instance let us consider the mean value of (x^2)in the ground state. In our units the virial theorem implies

⟨ (x^2)⟩/2 = 1/2 (Subscript[E, 0])= 1/4;  ⇒  ⟨(x^2)⟩ = 1/2

Using the function scalar and the eigenstate function just computed :

meanx2 = scalar[ξgrid,λ,ψ,(xgrid^2) ψ]

Example 2

As a second example we take a radial problem, the Coulomb potential. The only difference in the procedure is the word "rad" (or something else  you might prefer to use) attached to the call getD2 and normalize.

V[r_]:= ((- 1)/r);

dim=100;

Prototype computation

{ξgrid,xgridL1,D2} = getD2[dim,"rad"];

λ = 10.0;

ene1=Sort[Eigenvalues[((- 1)/(2(λ^2)))D2 +DiagonalMatrix[ V[λ xgridL1]]]][[1;;10]]

There is a  good agreement with known values - 1/2 (n^(2 ))as is easily checked:

ene1*((Range[10])^2)

λ dependence

We can look at the λ dependence of the eigenvalues

rangeL = Range[2,30,2];
dim=100;
numEig=Floor[dim/3];

{ξgrid,xgridL1,D2} = getD2[dim,"rad"];

eneL= Table[Sort[Eigenvalues[((-1)/(2(λ^2)))D2 +DiagonalMatrix[ V[λ xgridL1]]]][[1;;numEig]],
{λ,rangeL}];

We plot the quantity (n^2)(Subscript[E, n])+ (n-1) to separate different eigenvalues in the graph:

ListPlot[Table[Transpose[{rangeL,eneL[[All,j]]*j^2+(j-1)}],{j,1,numEig}],Joined-> True,
PlotRange-> {{0,25},{-1,numEig+1}},AxesOrigin-> {0,0},AxesLabel-> {"λ","(n^2) (Subscript[E, n]) + (n-1)"}]

The region of stability is quite large, especially for lower eigenvalues. To appreciate this point let us print with 12 decimal digits the first two eigenvalues:

NumberForm[eneL[[All,1]],{14,12}]

NumberForm[eneL[[All,2]],{14,12}]

It is clear that in this region make no sense to try an optimization on λ, all values are equivalent within errors.

Let us take a central value,  λ = 15, and study the accuracy of the ground state eigenfunction.

λ=15;

{eigH,eigV}= Eigensystem[((-1)/(2(λ^2)))D2 +DiagonalMatrix[ V[λ xgridL1]]];
order = Ordering[eigH];
eigH = eigH[[order]][[1;;numEig]];
eigV=eigV[[order]][[1;;numEig]];

To get the normalized ψ it is sufficient a call to the function normalize, with last argument "rad", or something not zero.

ψ=normalize[ξgrid,λ,eigV[[1]],"rad"];

This is the theoretical reduced radial wave function for hydrogen atom

ψTh[r_, n_, L_] := r*  (2/(n^2))(Sqrt[(((n-L-1)!)/((n+L)!))])*
(((((2 r)/n)))^L) * (E^(-(r/n)))* LaguerreL[n-L-1, 2 L + 1,((2 r)/n)] ;

We make a plot for both theoretical and numerical values :

xgrid = λ xgridL1;

sign= Sign[ψ[[Floor[dim/2]-3]]];

Show[{Plot[ψTh[x,1,0],{x,0,15},PerformanceGoal-> "Speed",
PlotRange-> {{0,12},{0,0.8}},PlotStyle-> Red],
ListPlot[Transpose[{xgrid,sign*ψ}]]}]

Quantitatively the maximum of their difference is :

Norm[ψ*sign- ψTh[xgrid,1,0],∞]

Even the absolute normalization of ψ comes out correctly.  It does not depend on the knowledge of ψTh.