Gauss

Quantum Mechanics
 K. Konishi, G. Paffuti

Study the Gauss integration formulas using Mathematica .

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

Gauss integration formulas

General remarks

Given an interval, possible infinite, (a,b) and a weight function μ[x], Gauss - type formulas give an approximate evaluation of integrals with a measure μ[x]𝕕x:

(∫_a^b)μ[x] F[x]𝕕x ≃ (∑_(i=1)^n)(Subscript[w, i])F[(Subscript[x, i])].

The points (Subscript[x, i]) are the abscissas and the numbers (Subscript[w, i]) the weights. This formula is based on the theory of orthogonal polynomials. The reader is referred to books on numerical analysis for the proof and for more details.

We shall be interested in the integration of simplest functions, the polynomials. An integration formula will be said of accuracy k if it gives exact results for polynomials of degree ≤ k. We have at our disposal n abscissas and n weights, i.e. 2n parameters. The integral of a polynomial is fixed once we know the integrals of simple x powers, by linearity. Using the 2n parameters at our disposal we can at most reproduce powers up to 2n-1 (i.e. (x^k)starting from k=0).  Thus  the  integration formula has accuracy 2n-1

at best.

Gauss' integration formulas are of maximum degree, i.e. an n-point formula is exact for polynomials of degree up to 2n-1.

Integration formulas are useful as we can approximate integrals of any function, once this function has been approximated by a polynomial. The simplest polynomial approximation is due to Lagrange and is briefly reviewed in next subsection.

Integration and interpolation

Lagrange formula

Consider a grid of n points (Subscript[x, 1])... (Subscript[x, n]) and let (Subscript[f, i])= f[(Subscript[x, i])] the values of a function f[x] on these points. It is possible to write a polynomial ϕ of degree n-1 which interpolates the function. This is most simply defined using Lagrange elementary polynomials:

(Subscript[ℒ, i])[x]=(∏_(j≠i)^n)(((x-(Subscript[x, 1]))..(x-(Subscript[x, i-1]))(x-(Subscript[x, i+1])).. (x-(Subscript[x, n])))/(((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])))).

by inspection it is clear that these polynomials act like "delta functions" on the grid:

(Subscript[ℒ, i])[(Subscript[x, j])] = (Subscript[δ, ij])

These are polynomials of degree n-1. The Lagrange interpolating polynomial follows :

ϕ[x] = (∑_(i=1)^n) (Subscript[f, i]) (Subscript[ℒ, i])[x]

It is possible to write this polynomial in a more compact form, which will be useful in some formal manipulations. Let consider the symmetric n-degree monic (i.e. with coefficient 1 for the higher order term) polynomial with roots (Subscript[x, i])

π(x) = (∏_(i=1)^n)(x-(Subscript[x, i]))

We have

ϕ[x] = (∑_(i=1)^n)(( π(x))/((x-(Subscript[x, i]))π'((Subscript[x, i])))) (Subscript[f, i])

In fact  π (x)/(x-(Subscript[x, i]))π'((Subscript[x, i])) is evidently a polynomial of degree n-1 which vanishes for all points of the grid except (Subscript[x, i])and on these points take the value 1, as follows from de l'Hopital formula, i.e. it coincides with (Subscript[ℒ, i])(x).

Integration

From any grid of points it is possible to extract an integration formula. If μ(x) is the (possible) weight, i.e. how our measure is defined, we can use eq.() to write

∫μ(x)𝕕x f(x) ≃ ∫μ(x)𝕕x ϕ(x)= (∑_(i=1)^n)∫μ(x)𝕕x(( π(x))/((x-(Subscript[x, i]))π'((Subscript[x, i])))) (Subscript[f, i]) ≡  (∑_(i=1)^n) (Subscript[W, i])(Subscript[f, i]);

i.e. weights are given by

(Subscript[W, i])= ∫μ(x)𝕕x(( π(x))/((x-(Subscript[x, i]))π'((Subscript[x, i]))))≡  ∫μ(x)𝕕x (Subscript[ℒ, i])[x]

This formula in general is accurate for polynomials of degree n - 1, or n.  All usual formulas of Newton-Cotes type on uniform grids are formulas of this nature.

An important improvement, discovered by Gauss, is obtained by using non uniform grids. In particular if the grid is built from  the zeros of an orthogonal polynomial of degree n with respect the measure μ than the degree of accuracy is maximal, i.e. 2n-1.

We review below some simple properties of orthogonal polynomials before proceeding to the integration formulas themselves.

Orthogonal polynomials

Given a measure μ[x] there is an associated set of orthogonal polynomials which form a complete set of functions with respect to that measure. Let us call them p[j,x], j being the degree. They are orthogonal, and it is assumed that they are normalized:

(∫_a^b)μ[x](Subscript[p, j])[x] (Subscript[p, k])[x] = (Subscript[δ, jk]) ;

We assume also that (Subscript[p, 0])[x]=1.

As stated above it can be shown that if one chooses as abscissas the zeros of the polynomial of degree n (the n which appears in () ) then there exist weights (Subscript[w, i])such that formula () is exact for polynomials up to order 2n-1. A theorem assures that in [a,b] the polynomial of degree n has n real zeros in this interval.

With a change of variable one can always reduce a finite interval to [-1,1], or [0,1].  We have essentially three different cases and for each case an infinite range of weight functions. The most common and useful cases are listed below, with the name of associated polynomials

The roots of the n-th polynomial can be found by functions like NSolve in Mathematica but it is interesting to note that they can be found by formulating the problem in a quantum mechanical language, in this case the roots are the eigenvalues of x operator in the subspace of polynomials of degree n-1, this is briefly sketched in the following subsection and in problem [1].

Quantum mechanics

It is interesting to note, from Quantum Mechanics point of view,  that Gauss formulas can be understood as a finite dimensional approximation of the Hilbert space constructed on the linear space spanned by the polynomials, (𝕃^2)[μ;a,b], or with the usual measure 𝕕x, by basis functions (Subscript[f, n])= (Sqrt[μ[x]])p[n,x].

Let us consider the operator x, and assume that its matrix elements in this basis form an (infinite) matrix (Subscript[X, ij]), i,j=0,1,...  Let|(Subscript[x, i])⟩ be the (formal) eigenstates. Then a generic integral can be written as ( |1⟩ denote the identity polynomial )

∫μ[x] F[x]𝕕x = (∑_(i,j)) ⟨1|(Subscript[x, i])⟩⟨(Subscript[x, i])|F[x]|(Subscript[x, j])⟩⟨(Subscript[x, j])|1⟩ = (∑_i)((|⟨1|(Subscript[x, i])⟩|)^2) F[(Subscript[x, i])].

If we take a finite dimensional approximation of Hilbert space, i.e. we consider only polynomials of degree ≤ n-1, the (Subscript[x, i])'s are the n eigenvalues of the matrix (Subscript[X, ij]) and ⟨1|(Subscript[x, i])⟩ are the first component of the corresponding eigenvectors.  Formula () is the Gauss formula.  It can be shown, see problem [1] that (Subscript[x, i])are both the eigenvalues of (Subscript[X, ij]) and the roots of p[n,x].  The weights (Subscript[w, i]) are the square of first component of corresponding eigenvectors.

A crucial point in the proof, and in practical implementation of the algorithm, is that a polynomial basis necessarily satisfies a recursion relation of the form (for normalized polynomials):

x (Subscript[f, n]) = (Subscript[γ, n+1])(Subscript[f, n+1]) + (Subscript[δ, n]) (Subscript[f, n]) + (Subscript[γ, n-1]) (Subscript[f, n-1])

The associated matrix X is hermitian and tridiagonal :

(Subscript[X, ij]) = ({{, "\!\(\*SubscriptBox[\(f\), \(0\)]\)", "\!\(\*SubscriptBox[\(f\), \(1\)]\)", ".", ".", ".", "\!\(\*SubscriptBox[\(f\), \(n\)]\)", ""}, {"\!\(\*SubscriptBox[\(f\), \(0\)]\)", "\!\(\*SubscriptBox[\(δ\), \(0\)]\)", "\!\(\*SubscriptBox[\(γ\), \(1\)]\)", "0", ".", ".", ".", "."}, {"\!\(\*SubscriptBox[\(f\), \(1\)]\)", "\!\(\*SubscriptBox[\(γ\), \(1\)]\)", "\!\(\*SubscriptBox[\(δ\), \(1\)]\)", "\!\(\*SubscriptBox[\(γ\), \(2\)]\)", ".", ".", ".", "."}, {".", "0", ".", ".", ".", ".", ".", "."}, {".", ".", ".", ".", ".", ".", ".", "."}, {".", ".", ".", ".", ".", ".", "\!\(\*SubscriptBox[\(γ\), \(n\)]\)", "0"}, {"\!\(\*SubscriptBox[\(f\), \(n\)]\)", ".", ".", ".", ".", "\!\(\*SubscriptBox[\(γ\), \(n\)]\)", "\!\(\*SubscriptBox[\(δ\), \(n\)]\)", "."}, {"", ".", ".", ".", ".", "0", ".", "."}})

This allows an extremely simple implementation of the Gauss formulas : it is sufficient to use the Eigenvalues or Eigensystem commands. It is also possible to give explicit analytical formulas for the weights, and this can be useful in some cases.

In the following sections we will show how the formalism works for certain class of polynomials.

N.B. We are using as a lowest order polynomial p[x] = 1 and this means that our measure is normalized, i.e. its integral is 1. Sometimes it can be convenient to adopt another normalization, for example for Laguerre polynomials one use Exp[-(x^2)] dx or for Legendre polynomials in [-1,1] one use dx. In these cases, if (Subscript[μ, 1]) is the normalized measure we have instead:

∫μ 𝕕x = c ;  Laguerre:  (∫_0^∞)Exp[-(x^2)]dx = (Sqrt[π]);   Legendre:  (∫_(-1)^1)1dx = 2

and it suffices to multiply normalized weights by the constant c.

Recursive relations

In literature recursive relations are usually presented for non normalized polynomials. In the form

(Subscript[p, i+1])[x] = ((Subscript[A, i])x + (Subscript[B, i]))(Subscript[p, i])[x] - (Subscript[C, i])(Subscript[p, i-1])[x];    i ≥ 1.

Defining the normalization as

(Subscript[h, i]) = ((Subscript[p, i]),(Subscript[p, i])) = ∫μ 𝕕x (((Subscript[p, i])[x])^2).

we have explicitly

(Subscript[A, i])= ((Subscript[a, i+1])/(Subscript[a, i]));     (Subscript[B, i])= ((Subscript[a, i+1])/(Subscript[a, i]))(((Subscript[b, i+1])/(Subscript[a, i+1]))- ((Subscript[b, i])/(Subscript[a, i]))) ;     (Subscript[C, i])= (((Subscript[a, i+1]) (Subscript[a, i-1])(Subscript[h, i]))/((Subscript[a, i]^2) (Subscript[h, i-1]))).

where

(Subscript[a, i])= coeff. max. degr. (i.e. i) in   (Subscript[p, i]);     (Subscript[b, i]) coeff. of term  (x^(i-1))  in  (Subscript[p, i]).

In this notation the matrix J, whose eigenvalues give the roots, has the form

J = (({{(Subscript[α, 1]), (Subscript[β, 1]), ., ., ., .}, {(Subscript[β, 1]), (Subscript[α, 2]), ., ., ., .}, {., ., ., ., ., .}, {., ., ., ., ., .}, {., ., ., ., ., (Subscript[β, n-1])}, {., ., ., ., (Subscript[β, n-1]), (Subscript[α, n])}}))  ;              where:   |({{(Subscript[α, i])= ((Subscript[b, i-1])/(Subscript[a, i-1]))- ((Subscript[b, i])/(Subscript[a, i]))}, {(Subscript[β, i])= ((Subscript[a, i-1])/(Subscript[a, i]))(Sqrt[((Subscript[h, i])/(Subscript[h, i-1]))])}})  .

Christoffel - Darboux relation

A very useful identity, which can be proven by induction and using the recursion relation () is the following, usually known as the Christoffel - Darboux relation:

(x-ξ) (∑_(i=0)^(n-1))(1/(Subscript[h, i])) (Subscript[p, i])[x](Subscript[p, i])[ξ] = ((Subscript[a, n-1])/((Subscript[a, n])(Subscript[h, n-1])))((Subscript[p, n])[x](Subscript[p, n-1])[ξ] - (Subscript[p, n])[ξ](Subscript[p, n-1])[x]).

The weights

For orthogonal polynomials it is possible to write in closed form the weights of the Gaussian integration.

Let (Subscript[x, i])the roots of the n-th polynomial (Subscript[p, n]). The π(x) function is nothing but  (Subscript[p, n]) normalized in such a way that its highest power has coefficient unity, i.e. the monic form of (Subscript[p, n]):

π(x) = (1/(Subscript[a, n]))(Subscript[p, n])(x)

By substitution in  () we have

(Subscript[W, i]) = (1/((Subscript[p, n]^('))[(Subscript[x, i])]))∫μ(x)𝕕x (((Subscript[p, n])[x])/((x-(Subscript[x, i]))))

Consider the Christoffel - Darboux relation for ξ = (Subscript[x, i]). As this is a root of (Subscript[p, n]) one has  (Subscript[p, n])[ξ]=0 and

(((Subscript[p, n])[x])/(x-(Subscript[x, i]))) = (((Subscript[a, n])(Subscript[h, n-1]))/(Subscript[a, n-1]))(1/((Subscript[p, n-1])[(Subscript[x, i])]))(∑_(i=0)^(n-1))(1/(Subscript[h, i])) (Subscript[p, i])[x](Subscript[p, i])[(Subscript[x, i])]

If one integrates this relation with respect to the μ measure one obtains a term proportional to the weight on the left hand side.  On the right hand side only the constant term (Subscript[p, 0])survives, due to the orthogonality relations. This integral, by definitions of the h's values is ((Subscript[p, 0]) is a constant, equal to (Subscript[a, 0]))

∫μ(x)𝕕x (Subscript[p, 0]) = ∫μ(x)𝕕x (Subscript[a, 0]) = ((Subscript[h, 0])/(Subscript[a, 0]))

(Subscript[W, i]) = (((Subscript[a, n])(Subscript[h, n-1]))/(Subscript[a, n-1]))(1/((Subscript[p, n]^('))[(Subscript[x, i])](Subscript[p, n-1])[(Subscript[x, i])]))(1/(Subscript[h, 0]))((Subscript[h, 0])/(Subscript[a, 0]))= (((Subscript[a, n])(Subscript[h, n-1]))/((Subscript[a, n-1])(Subscript[a, 0])))(1/((Subscript[p, n]^('))[(Subscript[x, i])](Subscript[p, n-1])[(Subscript[x, i])]))

This formula can be simplified in particular cases using relations for the derivative of orthogonal polynomials.

Gauss - Lobatto's  extension

Gauss' formulas give integrals in term of the grid points internal to the interval of integration [a, b]. For some applications it can be useful to extend the grid on some fixed points (Subscript[y, α]). It can be shown that in this case one has to construct a sequence of orthogonal polynomials with respect to  μ(x) r(x), where r(x) is the polynomial (x-(Subscript[y, 1])) ... (x- (Subscript[y, m])), to obtain maximal degree integrations. Formulas of this mixed nature are very useful in numerical integrations and some of the routines used internally by Mathematica are based on such an extension (the Konrod formulas). Here we shall be interested in the simplest case, where the extension concerns the boundary of the interval. The associated formulas are called the Gauss - Lobatto formulas.

As abscissas are fixed a priori the addition of a point in the grid will add only 1 in accuracy (we can fix only the weight). For a finite interval [a,b], the Gauss Lobatto formula (with n internal points) will have an accuracy (2n-1)+2 = 2n+1. For a seminfinite interval, [0,∞], the addition of 0 to the grid will produce an accuracy of (2n-1)+1 = 2n.

Weights can be computed by the general formula (), but in most cases this can explicitly be  transformed to formulas like ().  Below we shall give Gauss Lobatto roots and weight for most important class of polynomials, making an explicit computation only for the Laguerre polynomials, a case where the Gauss Lobatto extension is not easily found in books.

Laguerre polynomials

The Laguerre polynomials are orthogonal with respect the measure μ[x] = Exp[-x] in the interval [0,∞], and they will be particularly useful in studying radial problems. In the normalization and conventions of Mathematica (which are not universally adopted):

(∫_0^∞)𝕕x Exp[-x] (Subscript[L, i])[x] (Subscript[L, j])[x] = (Subscript[δ, ij]).

Roots and weights

The Laguerre polynomials satisfy the following recursion formula

x (Subscript[L, n])[x] = -(n+1)(Subscript[L, n+1])[x] + (2n+1)(Subscript[L, n])[x] - n (Subscript[L, n-1])[x]

Then vectors δ and γ of formula () are for n points formula (i.e. polynomials up to degree n-1): (we use n=5 in the example below and call J the matrix).

The leading coefficients of the polynomials can be read directly from ():

(Subscript[a, n+1])=- (( 1)/(n+1))(Subscript[a, n]); (Subscript[a, 0])=1  ⇒ (Subscript[a, n])= ((((-1))^n)/(n!))

as a check :

Table[Coefficient[LaguerreL[k,x],(x^k)],{k,1,5}]

All norms (Subscript[h, i]) as defined in () are units and for weights we have, from eq.():

(Subscript[W, i]) =  (((Subscript[a, n])(Subscript[h, n-1]))/((Subscript[a, n-1])(Subscript[a, 0])))(1/((Subscript[L, n]^('))[(Subscript[x, i])](Subscript[L, n-1])[(Subscript[x, i])]))= -(1/n)(1/((Subscript[L, n]^('))[(Subscript[x, i])](Subscript[L, n-1])[(Subscript[x, i])]))=  -(1/n)((Subscript[x, i])/((Subscript[x, i])(Subscript[L, n]^('))[(Subscript[x, i])](Subscript[L, n-1])[(Subscript[x, i])]))

Using the relation

x (d/dx)(Subscript[L, n])[x] = n (Subscript[L, n])[x] - n (Subscript[L, n-1])[x]

for x = (Subscript[x, i]) (root of (Subscript[L, n])) we can write the weights as

(Subscript[W, i]) =  (1/(n^2))((Subscript[x, i])/(((Subscript[L, n-1])[(Subscript[x, i])])^2))

An alternative form is obtained using the recursion relation for x = (Subscript[x, i])

(n + 1) (Subscript[L, n+1])[(Subscript[x, i])] = n (Subscript[L, n-1])[(Subscript[x, i])]

giving rise to

(Subscript[W, i]) =(1/(((n+1))^2))((Subscript[x, i])/(((Subscript[L, n+1])[(Subscript[x, i])])^2))

This is the form which will be used explicitly below.

Numerical routines for roots and weights

n = 5;
δ = 2Range[0,n-1]+1; γ = -Range[1,n-1];

(J = Normal[SparseArray[{Band[{1,1}]-> δ,Band[{1,2}]-> γ,Band[{2,1}]-> γ}]])//MatrixForm

Let us check that its eigenvalues are indeed the roots of L[n, x];

Sort[Eigenvalues[N@J]]==( x/.NSolve[LaguerreL[n,x]== 0,{x}])

Now we write two different functions:  one computes weights with eigenvectors, the other with analytical expressions

laguerreRootsEig[n_]:=Module[{J,δ,γ,radEig,w,eigs,eigV},
δ = 2Range[0,n-1]+1; γ = -Range[1,n-1];
J = Normal[SparseArray[{Band[{1,1}]-> δ,Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
{eigs,eigV}=Eigensystem[N@J];
order = Ordering[eigs];
radEig = eigs[[order]];
w=eigV[[order,1]]^2;
{radEig,w}]

Using explicit formulas for the weights :

laguerreRoots[n_]:=Module[{J,δ,γ,radEig,w},
δ = 2Range[0,n-1]+1; γ = -Range[1,n-1];
J = Normal[SparseArray[{Band[{1,1}]-> δ,Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
radEig=Sort[Eigenvalues[N@J]];
w=(radEig/((((n+1))^2)((LaguerreL[n+1,radEig])^2)));
{radEig,w}]

The two results are the same, within our precision (the reader consult Mathematica help for MapThread ). The function compare written below gives the ∞ norm, i.e. the maximum of absolute value of the difference between two list, Max[Abs[ a- b]]. We use Flatten to transform matrices in vectors, if needed.

compare[a_,b_]:= Norm[Flatten[a-b],∞]

MapThread[compare,{laguerreRootsEig[5],laguerreRoots[5]}]

Lagrange polynomials

We could in principle compute the weights using the general formula (). Severe rounding errors appear using Lagrange polynomials for high degrees in that form, but it is a good exercise to check these formulas for low degrees.

Lagrange elementary polynomials (Subscript[ℒ, i]) are easily constructed by defining an interpolation polynomial with value 1 on (Subscript[x, i]) and zero on other grid points.

This is the result obtained with our routine:

n = 3;
{xgrid, weights} = laguerreRoots[n];

These are Lagrange interpolating polynomials

ℒ[x_]:=Table[InterpolatingPolynomial[Transpose[{xgrid,SparseArray[{k}-> 1,{n}]}],x],{k,1,n}];
SetAttributes[ℒ,Listable];

From equation ():

weightsFromLagrange = Table[Integrate[ℒ[x][[k]] Exp[-x],{x,0,∞}],{k,1,n}];

to be compared to previous result

compare[weights,weightsFromLagrange]

They correctly agree.

Stability

Let us show at first why it is much better to use this matrix formalism instead of simply calling NSolve for computing roots. For small n the two procedures are equivalent:

n = 10;
radNS10=x/.NSolve[LaguerreL[n,x]== 0,x];
rad10 = laguerreRootsEig[n][[1]];

compare[radNS10,rad10]

For high enough n, NSolve with usual precision settings, has some trouble, for example you can get imaginary results:

n = 100;
radNS100=x/.NSolve[LaguerreL[n,x]== 0,x];
radNS100[[15;;20]]

The problem  can be solved by working with higher precision, but the time rapidly grows:

(radNS100 = x/.NSolve[LaguerreL[n,x]== 0,x,WorkingPrecision-> 55];)//Timing

The results are real :

radNS100[[15;;20]]

and they coincide within a reasonable precision with those obtained by our functions :

rad100 = laguerreRootsEig[n][[1]];//Timing

Max[Abs[rad100-radNS100]]

In the following we present a case in which it is better to use LaguerreRoots than LaguerreRootsEig.

Some integrals

Maximum degree

First of all let us check the maximum degree property of the Gauss Laguerre formulas. Integrals of monomials are easy to work out

(∫_0^∞)Exp[-x] 𝕕x (x^k) = k!

Let us consider n = 4. The formulas must be correct up to 2 n - 1 = 7, i.e. up to integrals of (x^7) :

Table[k!,{k,0,8}]

n = 4;
{roots,weights} = laguerreRoots[n];
Table[(roots^k).weights,{k,0,8}]

This is the expected result, the last integral is only approximate, as it is supposed to be.

Integrals of non polynomial functions

Let us check the formula with respect to the precision for some integrals. We will do always a double check, using both versions of weights:

n=10;

{rootsEig,weightsEig} = laguerreRootsEig[n];
{roots,weights} = laguerreRoots[n];

f0[x_]:= x^5;

{weightsEig.f0[rootsEig],weights.f0[roots], Integrate[Exp[-x]f0[x],{x,0,∞}]}

i.e. in this case we have an exact result. Let us take another function

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

{weightsEig.f0[rootsEig],weights.f0[roots], NIntegrate[Exp[-x]f0[x],{x,0,∞}]}

We see an error of about (10^(-4))which can be lowered augmenting n

n = 50;
{rootsEig,weightsEig} = laguerreRootsEig[n];
{roots,weights} = laguerreRoots[n];
{weightsEig.f0[rootsEig],weights.f0[roots], NIntegrate[Exp[-x]f0[x],{x,0,∞}]}

It is obvious that one can also compute integrals without μ[x], just considering f Exp[x] in the weight formula:

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

{(weightsEig*Exp[rootsEig]).f0[rootsEig],(weights*Exp[roots]).f0[roots], NIntegrate[f0[x],{x,0,∞}]}

While laguerreRoots gives the correct result, laguerreRootsEig gives a completely wrong result. To understand why let us have a look at weights, and at the same quantities multiplied by the exponentials:

weightsEig

weightsEig*Exp[rootsEig]

We see that a relatively small error in weights, computed by Eigensystem, is enhanced by 80 orders of magnitude in the last entry.

This is why it is safer, in these cases, to consider the analytical version of the computation.

Extended Laguerre grids

For various applications it can be useful to define a set of grid points which include origin, i.e. to perform a Gauss-Lobatto extension of the grid. Following the ideas exposed in the first section we have to construct a series of polynomials orthogonal with respect to the measure Exp[-x] x, x being the polynomial whose root is the new point to be added. There is no particular work to be done, these polynomials are, as it is well known, generalized Laguerre polynomials of index 1. In general

(∫_0^∞)Exp[-x] 𝕕x  (x^λ) (Subscript[L, n]^λ)[x](Subscript[L, m]^λ)[x] = 0 if n≠ m;

With the normalization used by Mathematica

(∫_0^∞)Exp[-x] 𝕕x  (x^λ) (((Subscript[L, n]^λ)[x])^2)= (((n+λ)!)/(n!))

These formulas, with λ = 1, give us polynomials and normalization factors (Subscript[h, i]), needed for computing weights.

For n interior points and one point at the border, i.e. x=0, the grid to be chosen must be composed by 0 and the roots of (Subscript[L, n]^1)[x]. To compute the roots of Laguerre polynomial we act as before.

The general recursion relation for Laguerre polynomials reads

(n+1) (Subscript[L, n+1]^α)[x] - (2n + α + 1 - x)(Subscript[L, n]^α)[x] + (n+α)(Subscript[L, n-1]^α)[x] = 0

We can use the formulation () or define a normalized basis. Let us follow this last way, as we have done for the usual Laguerre basis. The normalized polynomials are

(Subscript[f, k]) = (1/(Sqrt[k]))(Subscript[L, k-1]^1)[x];  k = 1,2,...

With α = 1 and with the shift n → n-1 relation () can be transformed in

x (Subscript[f, n]) = 2n (Subscript[f, n]) - (Sqrt[n(n+1)])(Subscript[f, n+1]) - (Sqrt[n(n-1)])(Subscript[f, n-1]) ;

which defines vectors δ and γ to be used to define matrix J giving roots as eigenvalues, see eq.().

Let us pass to weights. First we note that from () we can compute coefficient of highest order term (Subscript[a, n]) in (Subscript[L, n]^1):

(n+1)(Subscript[a, n+1]) =- (Subscript[a, n])  ⇒ (Subscript[a, n]) =(((-1))^n) (1/(n!));   (Subscript[a, 0]) = 1.

As a check:

Table[Coefficient[LaguerreL[k,1,x],(x^k)],{k,1,5}]

We have used (Subscript[L, 0]^1)[ x] = 1. The normalizations are as in ():

(Subscript[h, n]) = n+1.

To compute the weights we start from the general definition (). If (Subscript[x, i])are the roots of (Subscript[L, n]^1) the characteristic polynomial π(x), with 0 added as a root, is

π(x) = x (∏_(i=1)^n)(x-(Subscript[x, i])) = (x/(Subscript[a, n]))(Subscript[L, n]^1)[x] ≡ (x/(Subscript[a, n]))(Subscript[p, n])[x]

We write  (Subscript[p, n])for the polynomial for typographical reasons. Let us consider separately internal points weights and first weight (Subscript[W, 0])( relative to x=0).  For i ≥ 1 we have

(Subscript[a, n])π'((Subscript[x, i])) = (Subscript[p, n])[(Subscript[x, i])] + (Subscript[x, i]) (Subscript[p, n]^('))[(Subscript[x, i])] = (Subscript[x, i])(Subscript[p, n]^('))[(Subscript[x, i])]

being (Subscript[x, i]) a root of (Subscript[p, n]). For (Subscript[x, 0])=0 instead

(Subscript[a, n])π'((Subscript[x, 0])) = (Subscript[p, n])[0] = n+1

The value in 0 can be checked from the recursion relation, or from the explicit expression for the polynomials:

(Subscript[L, n]^λ)(x)=((Γ(n+λ+1) )/(Γ(n+1)))(∑_(k=0)^n)(((Subscript[(-n), k]) (x^k))/(Γ(k+λ+1) k!)) ;         (Subscript[(a), k])=a (a+1) (a+2)...(a+k-1);    ⇒     (Subscript[L, n]^λ)(0)=((Γ(n+λ+1))/(Γ(n+1) Γ(λ+1)))=(((n+λ)!)/(n! λ!))=(({{n}, {λ}})).

For the weights we have then to compute (for i ≥ 1)

(Subscript[W, i])=  ∫μ(x)𝕕x((x (Subscript[p, n])[x])/((x-(Subscript[x, i]))(Subscript[x, i])(Subscript[p, n]^('))((Subscript[x, i]))))

Now note that ω (x) = x μ[x] = x Exp[-x] is exactly the measure with respect to (Subscript[p, n]) is an orthogonal polynomial (they have been constructed with such a task), and this is the measure which is used to define norms etc. So let us write:

(Subscript[W, i])=  ∫ω(x)𝕕x(((Subscript[p, n])[x])/((x-(Subscript[x, i]))(Subscript[x, i])(Subscript[p, n]^('))((Subscript[x, i]))))

This integral can be computed using the Christoffel - Darboux relation. For ξ = (Subscript[x, i]) we have see (), ()

(((Subscript[p, n])[x])/(x-(Subscript[x, i]))) = (((Subscript[a, n])(Subscript[h, n-1]))/(Subscript[a, n-1]))(1/((Subscript[p, n-1])[(Subscript[x, i])]))(∑_(i=0)^(n-1))(1/(Subscript[h, i])) (Subscript[p, i])[x](Subscript[p, i])[(Subscript[x, i])] ⇒ (Subscript[W, i]) = (1/(Subscript[x, i])) (((Subscript[a, n])(Subscript[h, n-1]))/((Subscript[a, n-1])(Subscript[a, 0])))(1/((Subscript[p, n]^('))[(Subscript[x, i])](Subscript[p, n-1])[(Subscript[x, i])])).

Inserting the computed values for the coefficients

(Subscript[W, i])= -(1/n)n (1/((Subscript[x, i])(Subscript[p, n]^('))[(Subscript[x, i])](Subscript[p, n-1])[(Subscript[x, i])]))= -(1/((Subscript[x, i])(Subscript[p, n]^('))[(Subscript[x, i])](Subscript[p, n-1])[(Subscript[x, i])]))

For the derivatives of general Laguerre polynomials

x (d/dx)(Subscript[L, n]^α)[x] = n (Subscript[L, n]^α)[x]- (n+α) (Subscript[L, n-1]^α)[x] ⇒ (Subscript[x, i])(Subscript[p, n]^('))[(Subscript[x, i])] = - (n+1)(Subscript[p, n-1])[(Subscript[x, i])]

and

(Subscript[W, i])= (1/(n+1))(1/(((Subscript[p, n-1])[(Subscript[x, i])])^2))

From () on grid points

(n+1)(Subscript[L, n+1]^1)[(Subscript[x, i])] = -(n+1)(Subscript[L, n-1]^1)[(Subscript[x, i])]

and we can write

(Subscript[W, i])= (1/(n+1))(1/(((Subscript[L, n+1]^1)[(Subscript[x, i])])^2))

For the root x = 0 we have, from ()

(Subscript[W, 0])=  ∫μ(x)𝕕x((x (Subscript[p, n])[x])/(x(n+1)))=  ∫μ(x)𝕕x(( (Subscript[p, n])[x])/((n+1)))= (1/(n+1))(∫_0^∞)𝕕x Exp[-x] (Subscript[L, n]^1)[x]

The integral can be done with integration by parts. Using

(Subscript[L, n]^1)[x] =- (d/dx)(Subscript[L, n+1])[x]

and orthogonality

(Subscript[(Subscript[W, 0]) =  - (1/(n+1))(∫_0^∞)𝕕x Exp[-x](d/dx)(Subscript[L, n+1])[x]=  - (1/(n+1))(Exp[-x](Subscript[L, n+1])[x])|, 0]^∞) + (1/(n+1))(∫_0^∞)𝕕x Exp[-x](Subscript[L, n+1])[x] = (1/(n+1))(Subscript[L, n+1])[0] = (1/(n+1))

Numerical routines for roots and weights

From equations (), () and () we can write the following function for roots and weights on the extended grid

laguerreRootsL1[n_]:=Module[{J,δ,γ,radEig,w,eigs},
δ = 2Range[1,n]; γ = -(Sqrt[Range[n-1]Range[2,n]]);J = Normal[SparseArray[{Band[{1,1}]-> δ,Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
radEig=Sort[Eigenvalues[N@J]];
w= (1/(n+1))(1/((LaguerreL[n+1,1,radEig])^2));
w=Prepend[w,N[1/(n+1)]];
radEig = Prepend[radEig,0];
{radEig,w}]

We can also write a function which computes directly the weights from their definition through Lagrange interpolating polynomials, as we have done in Laguerre case:

rootCardLag[dim_,m_]:=Module[{ℒ,xgrid,dim1,w1,xgrid1},
xgrid= x/.NSolve[LaguerreL[dim,m,x]== 0,{x}];
dim1 = dim+1;
xgrid1 = PadLeft[xgrid,dim1];
ℒ[x_]:=Table[InterpolatingPolynomial[Transpose[{xgrid1,SparseArray[{k}-> 1,{dim1}]}],x],{k,1,dim1}];
SetAttributes[ℒ,Listable];
w1=Table[NIntegrate[ℒ[x][[k]] Exp[-x],{x,0,∞}],{k,1,dim1}];
{xgrid1,w1}]

They agree:

n=4;

compare[laguerreRootsL1[n],rootCardLag[n,1]]

and the roots are correct :

x/.NSolve[LaguerreL[n,1,x]== 0,{x}]

laguerreRootsL1[n][[1]]

We note that the roots and weights are not the same as a simple extension of the internal Laguerre grid. In fact, rootCardLag can compute these quantities also for the usual grid, passing 0 as second index of Laguerre polynomials (we put off a warning message which is issued by NIntegrate as the first integral is zero (exactly) than to a better precision than the default one of NIntegrate)

Off[NIntegrate::"ncvb"]
rootCardLag[n,0]
On[NIntegrate::"ncvb"]

Integrals

The only check we perform is to show that the accuracy has effectively grown from 2 n - 1 to 2 n. Let us consider the test we made on the  internal points of Gauss Laguerre grid:

(∫_0^∞)Exp[-x] 𝕕x (x^k) = k!

If we take n = 4,   we expect correct results up to k=7 for the old internal grid and up to k=8 for the new extended grid:

n=4;

Table[k!,{k,0,9}]

Old grid

{roots,weights} = laguerreRoots[n];
Table[(roots^k).weights,{k,0,9}]

New grid

{roots,weights} = laguerreRootsL1[n];
Prepend[Table[(roots^k).weights,{k,1,9}],Total[weights]]

We see that in fact the k = 8 value is exact, while for k=9 both grids give only approximate results.

Hermite polynomials

We can repeat without modification the previous analysis for the Hermite polynomials, the ones appearing in the solution of harmonic oscillators. The polynomials are orthogonal with respect the measure Exp[-(x^2)] but with the usual definition (used also in Mathematica ) they are normalized as:

(∫_(-∞)^∞)𝕕x Exp[-(x^2)] (Subscript[H, n])[x] (Subscript[H, m])[x] = (Subscript[δ, nm]) (2^n) n! (Sqrt[π]).

Hermite polynomials satisfy the recursion relation

x H[n,x] = (1/2)H[n+1,x] + n H[n-1,x]

To keep things simple let us define normalized polynomials

(Subscript[p, n])[x] = (1/(π^(1/4)))(1/(2^(n/2)))(1/(Sqrt[n!]))(Subscript[H, n])[x]

For these polynomials

x (Subscript[p, n])[x] = ((Sqrt[n+1])/(Sqrt[2]))(Subscript[p, n+1])[x] + ((Sqrt[n-1])/(Sqrt[2]))(Subscript[p, n-1])[x]

We recognize here the matrix elements of x operator for a harmonic oscillator.  We can write at once (we include the (Sqrt[π])only at the final stage of the normalization for weights)

hermiteRootsEig[n_]:= Module[{J,δ,γ,radEig,w,eigs,eigV},
 γ = Sqrt[Range[1,n-1]/2];
J = Normal[SparseArray[{Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
{eigs,eigV}=Eigensystem[N@J];
order = Ordering[eigs];
radEig = eigs[[order]];
w=(Sqrt[π])eigV[[order,1]]^2;
{radEig,w}]

hermiteRoots[n_]:= Module[{J,δ,γ,radEig,w},
  γ = Sqrt[Range[1,n-1]/2];J = Normal[SparseArray[{Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
radEig =Sort[ Eigenvalues[N@J]];
w=(((2^(n-1))(Sqrt[π]))/(n^2))((n!)/((HermiteH[n-1,radEig])^2));
{radEig,w}]

Some checks

Consistency and agreement with zeros of Hermite polynomials :

compare[hermiteRootsEig[5],hermiteRoots[5]]

compare[(x/.NSolve[HermiteH[5,x]== 0,x]),hermiteRoots[5][[1]]]

Integrals

This one must be exact :

n = 20;
f[x_]:= (x^6)

{roots,w}= hermiteRoots[n];

NumberForm[{NIntegrate[Exp[-(x^2)]f[x],{x,-∞,+∞}],w.f[roots]},10]

An approximate integral

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

NumberForm[{risE=NIntegrate[Exp[-(x^2)]f[x],{x,-∞,+∞}],w.f[roots]},10]

This is the convergence rate :

TableForm[Table[{roots,w}= hermiteRoots[n];NumberForm[#1,10]&/@{n,w.f[roots],risE},{n,20,120,20}],TableHeadings-> {None,{"n","Gauss","exact"}}]

Legendre Polynomials

The Legendre polynomials are orthogonal in the interval [-1,1] and, as is well known from angular momentum theory, satisfy:

(∫_(-1)^(+1))𝕕x (Subscript[P, n])[x](Subscript[P, m])[x] = (2/(2n+1)).

These polynomials satisfy the recursion relation

x  (Subscript[P, n])[x] =(( n+1)/(2n+1))(Subscript[P, n+1])[x] + (n/(2n+1))(Subscript[P, n-1])

For normalized polynomials p[k,x] = LegendreP[k,x] (Sqrt[(2k+1)/2])

x  (Subscript[P, n])[x] =(( n+1)/(Sqrt[(2n+1)(2n+3)]))(Subscript[P, n+1])[x] + (n/(Sqrt[(2n+1)(2n-1)]))(Subscript[P, n-1])

The function legendreRoots is similar to previous ones, the factor 2 which is needed if one use the usual measure dx, is added for the weights:

legendreRootsEig[n_]:= Module[{J,δ,γ,radEig,w,eigs,eigV},
 γ = Table[(( j+1)/(Sqrt[(2j+1)(2j+3)])),{j,0,n-2}];
J = Normal[SparseArray[{Band[{2,1}]-> γ,Band[{1,2}]-> γ}]];
{eigs,eigV}=Eigensystem[N@J];
order = Ordering[eigs];
radEig = eigs[[order]];
w=2eigV[[order,1]]^2;
{radEig,w}]

{rad,w} = legendreRootsEig[6]

This is the original Gauss rule, usually reported for the interval [0,1]. To change the interval form [0,1] to [a,b] it is sufficient the change of variables

x → a + (b-a) x;    roots -> a+(b-a)roots;    w -> (b-a)w

This rule is used as a particular submethod by NIntegrate and Gauss points and weights can be extracted as shown below (the second argument of GaussRuleData is the precision):

myGR[n_,a_,b_]:= {a+(b-a)#[[1]],(b-a)#[[2]]}&[Take[NIntegrate`GaussRuleData[n,6],2]]

This can be compared with our calculation :

myGR[6,-1,1]

Chebyshev polynomials

Chebyshev is one of the possible transcrptions of a cyrillic name, another one is Tsebysheff.

These polynomials are of particular importance, especially for finite domain problems, and are, in a sense, a Fourier series in disguise. There are two types of Chebyshev polynomials, T and U, in Mathematica ChebyshevT[n,x] and ChebyshevU[n,x]. Their definition is

(Subscript[T, n])[cos[θ]] =Cos[n θ]  ;    (Subscript[U, n])[cos[θ]] = ((Sin[(n+1)θ])/(Sin[θ])).

The respective measures are

μ[x] = (1/(Sqrt[1-(x^2)]));      μ[x] = (Sqrt[1-(x^2)]);

Their normalizations are easily computed by the change of variables x = Cos[θ]:

(∫_(-1)^1)((𝕕x)/(Sqrt[1-(x^2)]))  (Subscript[T, n])[x] (Subscript[T, m])[x] = (π/2) (Subscript[δ, nm]) ;  n≠ 0;        (∫_(-1)^1) 𝕕x (Sqrt[1-(x^2)])  (Subscript[U, n])[x] (Subscript[U, m])[x] = (π/2) (Subscript[δ, nm]) ;

For these polynomials there is no need of all our machinery, as the roots and weights are found by simple inspection of () and from the elementary properties of Fourier series. In any case let us carry out our usual computation:  it gives a check of the method in any case.

ChebyshevT

Let us have a look at the form of first polynomials

TableForm[Table[ChebyshevT[n,x],{n,0,4}]]

The roots of n - th polynomial are clearly given by n θ = (2k+1)π/2 i.e.

θ = (π/(2 n))(2k+1), k=0... n-1;    x = Cos[(π/(2 n))(2k+1)]

The recursion relation follows from

Cos[θ] Cos[n θ] = (1/2) Cos[(n+1)θ] + (1/2)Cos[(n-1)θ];

The weights, taking account the normalization condition () are given by

(Subscript[w, k]) = (π/n)

The Gauss formula is

(∫_(-1)^1)((𝕕x)/(Sqrt[1-(x^2)]))F[x] = (∑_(k=0)^(n-1))F[(Subscript[x, k])] (Subscript[w, k]) = (π/n)(∑_(k=0)^(n-1))F[(Subscript[x, k])]

In θ variables this formula reads

(∫_0^π)f[θ]dθ = (π/n)  (∑_(k=0)^(n-1))f[(π/(2 n))(2k+1)]

which is obviously true. The integral of identity gives π on both sides, while the integral of Cos[j θ] gives zero on both sides, for j≤ 2n-1:

Clear[n]

Simplify[Sum[ (π/n)  (∑_(k=0)^(n-1))Cos[j(π/(2 n))(2k+1)],{k,0,n-1}],Assumptions-> j∈Integers]

Let us verify these relations within our framework.

The recursion relation follows from

Cos[θ] Cos[n θ] = (1/2) Cos[(n+1)θ] + (1/2)Cos[(n-1)θ];

But there is a tricky point.  Let us factor out in the measure  π/2, which is the common factor for all (Subscript[T, n])with n≥1. Then

(∫_(-1)^1)((𝕕x)/(Sqrt[1-(x^2)])) = (∫_0^π) dθ = π = (π/2) 2;

i.e. if we want normalized states we have to take as first element (Subscript[f, 0])= 1/(Sqrt[2]). The first two recursion relations then read

x (Subscript[f, 0]) = (1/(Sqrt[2]))x = (1/(Sqrt[2]))(Subscript[f, 1]);   x (Subscript[f, 1])= (x^2) = (1/2)(2(x^2)-1)+(1/(Sqrt[2]))(Subscript[f, 0]);

Taking account of this fact and with the normalization condition () (i.e. a factor π for the weights, the integral of 1):

rootsChebyshevT[n_]:= Module[{J,δ,γ,radEig,w,eigs,eigV},J = Normal[SparseArray[{Band[{2,1}]-> 1/2,Band[{1,2}]-> 1/2},{n,n}]];
J[[1,2]]=J[[2,1]]=1/(Sqrt[2]);
{eigs,eigV}=Eigensystem[N@J];
order = Ordering[eigs];
radEig = eigs[[order]];
w=π eigV[[order,1]]^2;
{radEig,w}]

chebyshevTRoots[n_]:= {Sort[N@Table[Cos[(π/(2 n))(2k+1)],{k,0,n-1}]],(π/n)Table[1,{n}]}

compare[N@chebyshevTRoots[7],N@rootsChebyshevT[7]]

Some integrals

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

{roots,w} = chebyshevTRoots[20];

NumberForm[{risE=NIntegrate[(1/(Sqrt[1-(x^2)]))f[x],{x,-1,+1}],w.f[roots]},10]

As will be clear in the study of pseudospectral techniques, where the Chebyshev polynomials and other form of interpolation will be used, it is important that sampling points are not uniformly distributed:

Plot[f[x],{x,-1,1},Epilog-> Point/@Transpose[{roots,f[roots]}]]

ChebyshevU

Let us have a look at the form of the first polynomials

TableForm[Table[ChebyshevU[n,x],{n,0,4}]]

The roots of n - th polynomial are clearly given by (n+1) θ = k π i.e.

θ = (π/((n+1)))k, k=1... n;    x = Sin[(π/((n+1)))k]/Sin[θ].

Note that θ = 0 is not a root.

The recursion relation follows from

Cos[θ] Sin[(n+1) θ] = (1/2) Sin[(n+2)θ] + (1/2)Sin[n θ];

The weights, taking account the normalization condition () are given by

(Subscript[w, k]) = (π/n)

The same attention as before must be paid to n=0, but in this case:

(∫_(-1)^1)𝕕x (Sqrt[1-(x^2)])= π/2;

and all normalizations agree. We write the result:

rootsChebyshevU[n_]:= Module[{J,δ,γ,radEig,w,eigs,eigV},J = Normal[SparseArray[{Band[{2,1}]-> 1/2,Band[{1,2}]-> 1/2},{n,n}]];
{eigs,eigV}=Eigensystem[N@J];
order = Ordering[eigs];
radEig = eigs[[order]];
w=π eigV[[order,1]]^2;
{radEig,w}]

We write - Cos[] which sorts automatically the results, Cos[θ] being  a decreasing function in [0,,π].

chebyshevURoots[n_]:= {N@Table[-Cos[(π/( n+1))k],{k,1,n}],(π/n)Table[1,{n}]}

rootsChebyshevU[5]

N[chebyshevURoots[5]]

Gauss Lobatto points.

If we add the extrema of the interval we obtain the Gauss-Lobatto points. These are obviously the zeros of  (1-(x^2))(Subscript[T, n]^('))[x].On a generic interval

CGLGrid[x0_, L_, n_Integer /; n > 1] := x0 +(1/2) L(1 -  Cos[π Range[0,n-1]/(n-1)])

N@CGLGrid[-1,2,7]

rootsChebyshevU[5][[1]]

We can generate them with :

internalPoints[x0_,L_,n_Integer]:= x0 +(1/2)L(1-Cos[(π/( n+1))Range[1,n]])

internalPoints[1.4,6.5,6]

or dropping the first and last elements of CGLgrid

CGLGrid[1.4,6.5,8]

CGLGrid[1.4,6.5,8][[2;;-2]]

Fourier series

Our last example is probably the more important. Polynomials can be considered not only on ℝ but also in other spaces, for instance on circles, and this brings us directly to the Fourier series, i.e. to the expansion of periodic functions in terms of the trigonometric functions. To simplify the notation and to bring out the similarities with previous cases we adopt the exponential version of the series. Take an interval [0, L] and a periodic function on it, its Fourier expansion is defined as

F[x] = (∑_(k=-∞)^(+∞))(Subscript[c, k]) (1/(Sqrt[L]))Exp[i k x ((2π)/L)]

In quantum mechanical jargon the exponential are eigenfunctions of a free particle on the circle of length L.  Or with L=1 and x interpreted as an angular variable, k is the eigenvalue of (Subscript[L, z ]) and the exponential function  its eigenstates. Without any approximation

(Subscript[c, k])= (∫_0^L)F[x](1/(Sqrt[L]))Exp[-i k x ((2π)/L)]dx

From () one can extract an approximant of order N, which we take in a symmetric form:

(Subscript[S, N])[x] = (∑_(k=-N/2)^(N/2))(Subscript[a, k])(1/(Sqrt[L]))Exp[i k x ((2π)/L)]

Functions of this type are called trigonometric polynomials.  Introducing the complex variable

ζ = Exp[i  x ((2π)/L)]

they have a polynomial form, negative powers arise because ζ lies on a unitary circle, and in fact 1/ζ = (SuperStar[ζ]). We take N even to fix the ideas. In the sum () there appear N+1 terms, and in this sense they are the analog of a polynomial of degree N (in previous sections we considered polynomials of degree N-1). The same conclusion can be drawn by factorizing (ζ^(-N/2)), the resulting sum is a polynomial of degree N in ζ.

We can ask which is a Gauss formula based on this kind of polynomials, if there is one. The answer might look unexpected: it is just the elementary trapezoidal approximation to the integrals.

Consider a uniformly spaced grid of N+1 points (i.e. the degree of the polynomial augmented by 1, as in the previous sections, or the number of independent coefficients):

(Subscript[x, s]) = -(L/2)+ (L/(N+1))s;   s=0,2,... N

The point x = L is identified with x = 0, so if one wishes one can take the index in the range 1,...N+1. Consider at first the values of the trigonometric polynomials on this grid:

S[(Subscript[x, s])]= (∑_(k=-N/2)^(N/2))(Subscript[a, k])(1/(Sqrt[L]))Exp[i k (s/(N+1)) 2π ]≡ (∑_(k=-N/2)^(N/2)) (Subscript[U, sk]) (Subscript[a, k])

The (N+1)×(N+1) matrix U is unitary (apart a factor)

(∑_(k=-N/2)^(N/2))(Subscript[U, sk])(Subscript[U, jk]^( * )) =((N+1)/L) (Subscript[δ, sj]);      (∑_(s=0)^N)  (Subscript[U, sk])(Subscript[U, sq]^( * )) =((N+1)/L) (Subscript[δ, kq])

This follows from the trivial but fundamental identities

(∑_(k=0)^N)(ζ^k)=(Piecewise[{{((1 - (ζ^(N+1)))/(1-ζ)), ζ≠1}, {N+1, ζ=1}}])  ;     and   (∑_(k=-N/2)^(N/2))(ζ^k)= (ζ^(-N/2))(∑_(k=0)^N)(ζ^k)=(Piecewise[{{(ζ^(-N/2))((1 - (ζ^(N+1)))/(1-ζ)), ζ≠1}, {N+1, ζ=1}}])

Using () one can invert () :

(Subscript[a, k])=(L/(N+1))(∑_(s=0)^N) S[(Subscript[x, s])](1/(Sqrt[L]))Exp[-i k (s/(N+1)) 2π ]

Proof of ()

We have, in this notation:

(Subscript[U, sk])(Subscript[U, jk]^( * )) = (∑_(k=-N/2)^(N/2)) (1/L)Exp[i (s-j)((2π)/(N+1))k] ;   ζ = Exp[i (s-j)((2π)/(N+1))];   (ζ^(N+1))= 1

then

(Subscript[U, sk])(Subscript[U, jk]^( * )) = ((N+1)/L)(Subscript[δ, sj])

Similarly

(∑_(s=0)^N)  (Subscript[U, sk])(Subscript[U, sq]^( * )) =  (∑_(s=0)^N)(1/L)Exp[i (k-q)((2π)/(N+1))s];  ζ = Exp[i (k-q)((2π)/(N+1))];   (ζ^(N+1))= 1;

and

(∑_(s=0)^N)  (Subscript[U, sk])(Subscript[U, sq]^( * )) =((N+1)/L)(Subscript[δ, kq])

Gauss formulas

Consider now the integral of a trigonometric polynomial (). We have clearly, using ()

(∫_0^L)(Subscript[S, N])[x] 𝕕x = (Subscript[a, 0]) = (L/(N+1))(∑_(s=0)^N) S[(Subscript[x, s])]

Thus we can define a Gauss formula on a uniform grid, of N+1 points, as

(∫_0^L)F[x] 𝕕x = (∑_(s=0)^N)(L/(N+1)) F[(Subscript[x, s])];    (Subscript[w, s]) = (L/(N+1)).

The reader will recognize the weights of ChebyshevT, which are the Cos[] part of the Fourier series. In that case L = 2, and we consider N instead of N+1.

The abscissas (Subscript[x, s]) are the roots of a polynomial equation of degree N+1 on the circle (this is the analogous of ordinary roots for polynomials in ℝ):

(ζ^(N+1))= 1.

As in other cases Gauss formulas are not only true for trigonometric polynomials of degree N, but up to degree 2N. Indeed for an arbitrary degree M we have exactly

(∫_0^L)(Subscript[F, M])[x] 𝕕x = (Subscript[a, 0])

Consider now the right hand side of () for each term in (Subscript[F, M]). For (Subscript[a, 0])it gives the correct result, and the whole integral, so we have to show that all other terms contribute zero to the sum. For the k-th term:

(∑_(s=0)^N) Exp[i k ((2π)/L)(Subscript[x, s])]= (∑_(s=0)^N) Exp[i k 2π (s/(N+1))]= (∑_(s=0)^N) (ζ^s) ;  ζ = Exp[i k  ((2π)/(N+1))];

The sum is zero as long as ζ ≠ 1, i.e. up to k≤ N, which correspond to trigonometric polynomials of degree 2N (with 2N+1 coefficients).

Relation () is nothing but the trapezoidal approximation of integrals. The trapezoidal rule for an interval of length L divided in N+1 intervals of length h = L/(N+1)  ((Subscript[x, 0])=0 and (Subscript[x, N+1])= L) approximate the integral as the sum of the areas of the trapezia formed by the graph of f[x] and the x-axis:

(∫_0^L)f[x]𝕕x ∼ (h/2)(f[(Subscript[x, 0])]+f[(Subscript[x, 1])])+ (h/2)(f[(Subscript[x, 1])]+f[(Subscript[x, 2])])+... + (h/2)(f[(Subscript[x, N])]+f[(Subscript[x, N+1])])

Taking into account periodicity of f[x]:

(∫_0^L)f[x]𝕕x = h (∑_(s=0)^N)f[h s] = (L/(N+1))(∑_(s=0)^N)f[(L/(N+1)) s].

The same kind of analysis shows that formulas  () for the coefficient of the approximant are the trapezoidal approximation of exact integrals for Fourier series, and, by the previous analysis, these formulas are exact for trigonometric polynomials.