1. Theory

1.1. The pseudospectral method

In the pseudospectral method a function \(f(x)\) is approximated as a linear combination

\[f(x) \simeq f_N(x) = \sum_{j=0}^N f(x_j) g_j(x)\]

of cardinal functions \(g_j(x)\) where the expansion coefficients are values of \(f(x)\) at the grid points \(x_j\). For Lobatto-type grids, the interior grid points \(\{ x_j \}_{j=1}^{N-1}\) are defined as the zeros of \(P_N^\prime(x)\), where \(P_N(x)\) is a \(N\)-th order polynomial [1],

\[P_N^\prime(x_j) = 0, \ \ j=1,\cdots,N-1,\]

while the end (boundary) points \(x_0\) and \(x_N\) are given.

The cardinal functions \(g_j(x)\) have the property that

\[g_j(x_i) = \delta_{ij}, \label{delta_property_g}\]

and can be written as

\[g_j(x) = -\frac{1}{N(N+1)P_N(x_j)}\frac{(1-x^2)P_N^\prime(x)}{x-x_j}.\]

The integral of a function is approximated by quadrature

\[\int_{x_0}^{x_N} f(x) dx \simeq \sum_{j=0}^N w_j f(x_j), \label{quadrature_rule}\]

where \(w_j\) are the quadrature weights.

Footnotes

[1]

Which type of polynomial to choose depends on the problem. Examples include Chebyshev, Legendre, Laguerre or Hermite polynomials.

1.2. Gauss-Legendre-Lobatto

In the Gauss-Legendre-Lobatto pseudospectral method, \(P_N(x)\) is taken as the \(N\)-th order Legendre polynomial. The end points are \(x_0 = -1\) and \(x_N=1\), while the interior grid points must be found as the roots

\[P_N^\prime(x_i) = 0,\]

and the quadrature weights are given by

\[w_i = \frac{2}{N(N+1)P_N(x_i)^2}.\]

The grid points and weights can be determined using numpy functions (REF). Furthermore,

\[\begin{split}\frac{d g_j}{dx} \Bigr|_{x=x_i} &= g_j^\prime(x_i) = \tilde{g}_j^\prime(x_i) \frac{P_N(x_i)}{P_N(x_j)}, \\ \tilde{g}_j(x_i) &= \begin{cases} \frac{1}{4}N(N+1), \ \ i=j=0, \\ -\frac{1}{4}N(N+1), \ \ i=j=N, \\ 0, \ \ i=j \text{ and } 1 \leq j \leq N-1, \\ \frac{1}{x_i-x_j}, \ \ i \neq j, \end{cases}\end{split}\]

and for \(i,j = 1,\cdots,N-1\)

\[\begin{split}\frac{d^2 g_j}{dx^2} \Bigr|_{x=x_i} &= g_j^{\prime \prime}(x_i) = \tilde{g}_j^{\prime \prime}(x_i) \frac{P_N(x_i)}{P_N(x_j)}, \\ \tilde{g}_j^{\prime \prime}(x_i) &= \begin{cases} -\frac{1}{3} \frac{N(N+1)}{(1-x_i^2)}, & i = j,\\ -\frac{2}{(x_i-x_j)^2}, & i \neq j. \end{cases}\end{split}\]