Chebyshev/Lagrange polynomials in spectral methods

Question

I am currently trying to familiarise myself with (Pseudo-)Spectral Methods for solving differential equations. Now, I am struggling to understand some obviously crucial concept of this approach. The notes that I have read so far proceed their explanation of spectral methods roughly as follows:

The solution $f(x)$ to a differential equation $A f = b$ is approximated by an expansion in terms of basis functions $\phi_m(x)$ such that

$$f(x) \approx p_N(x) = \sum_{m=0}^N a_m \phi_m(x)$$

The notes then always highlight Chebyshev polynomials $T_m(x)$ as a particularly suitable choice for $\phi_m(x)$ in nonperiodic cases. Also the collocation points are chose as the roots of Chebyshev polynomials (I assume the roots of the highest order Chebychev polynomial?).

Now, all these notes then suddenly switch the representation of $p_N(x)$ to cardinal functions (e.g. Lagrange interpolations $C_m(x)$), such that

$$p_N(x) = \sum_{m=0}^N f(x_m) C_m(x)$$

I assume the major advantage here is that the Lagrange interpolations satisfy

$$C_m(x_i) = \delta_{im}$$,

and that therefore the requirement that the differential equation is satisfied at the collocation points uncouples the equations, such that it reduces to

$$A f(x_i) = b \quad , \qquad i = 0, \ldots, N$$

---

Is this true so far?

-> If so, then why do all my notes at first introduce Chebychev polynomials as a suitable choice if then Lagrange polynomials turn out to be a much more suitable choice. Why the detour?

-> If not, what am I missing here? Where does my understanding fail?

Answer 1

Your understanding is perfectly fine, except for the last statement that Lagrange polynomials turn out to be a more suitable choice. In fact, both methods, the modal and the nodal expansion, have their strengths and weaknesses.

Moreover, under certain conditions, both methods are unitarily equivalent. So you could start the educational journey from either one, and arrive at the other. Those conditions are exactly that your collocation nodes are the zeros of an orthogonal polynomial, and most lectures seem to commence from this fact and thus introduce orthogonal polynomials first.

But, as said, one could start also from the Lagrange perspective. There, in order to approximate the solution to $Af = b$, one would choose a suitable grid $\{x_k\}$, define Lagrange functions $L_k(x)$ on the grid, and project the equation onto the thus obtained space,

\begin{align} \sum_k c_k \langle L_j , A L_k \rangle = \sum_k b_k \langle L_j , L_k \rangle \end{align}

If you can then evaluate the scalar products (i.e. integrals), you have a matrix equation that can be solved by a computer, and everything's fine.

However, evaluating those integrals can be tedious (particularly when you go to higher dimensions), and the resulting matrices are not quite comfortable. Already the overlap or mass matrix $\langle L_j , L_k \rangle$ is not diagonal, which means that you get generalized eigenvalue problems and the like.

In order to simplify, one usually introduces approximation methods by choosing \begin{align} \int w(x) g(x) dx = \sum w_i g(x_i) \end{align} and with this the Lagrange polynomials become orthogonal $\langle L_j , L_k \rangle = \delta_{jk}$. Now, however, one made an approximation to the original equation the accuracy of which is not clear. (The collocation method, on the other hand, turns this into an assumption, and exclusively focuses on the gridpoints.)

This is where orthogonal polynomials kick in, with which one gets nice equations and high accuracy at the same time. For example, by choosing $\{x_k\}$ the abscissas of a Gaussian integration, the mass matrix is diagonal by construction (without approximation). Moreover, again due to the high integration accuracy, one can expand the orthogonal polynomials exactly into the Lagrange polynomials and thus attain the nodal representation, i.e. the two representations are unitary equivalent in this case.

So far this holds for any orthogonal polynomial. The reason why particularly the Chebyshev polynomial is chosen in practice is that it has a variety of nice properties. Among those is that one can obtain the abscissas and weights analytically and that one can use the Fourier transform for expansions, Clenshaw-Curtis integration, and so on.

Answer 2

To complete david's answer:

References: Canuto et al., Spectral Methods Fundamentals in Single Domains

---

We consider the Burgers equation

\begin{align} \text{advective form}: \qquad \frac{\partial u}{\partial t} &+ u\frac{\partial u}{\partial x} - \nu \frac{\partial^2 u}{\partial x^2} = 0, \hspace{1cm} (1)\\ \text{divergence form}: \qquad\frac{\partial u}{\partial t} &+ \quad \frac{\partial \mathcal{F}(u)}{\partial x} \hspace{0.7cm} = 0, \hspace{1.05cm} (2) \end{align} where the flux is given by $$ \mathcal{F}(u) = \frac{1}{2} u^2 - \nu \frac{\partial u}{\partial x}. $$

---

A polynomial expansion is given by $$ u(x)=\sum_0^N \hat{u}_k \phi_k(x), $$ where $\phi_k$ may be the Legendre$^1$, Chebyshev$^1$, Fourier$^1$ or even Lagrange$^2$ basis.

Furthermore we define

\begin{align} ^1 \mathbf{\text{modal}}~\text{values}:& \qquad \hat{u}_k \hspace{2.4cm} \text{(spectral)}, \\[1em] ^2 \mathbf{\text{nodal}}~\text{values}:& \qquad u_k \equiv u(x_k) \quad \text{(collocated)}. \end{align}

---

Historically we call an approximation using $(1)$ + Fourier basis:

(a) Galerkin approximation (aliasing free): $$ \frac{\text{d}\hat{u}_k}{\text{d}t}+\sum_{m+n=k} \hat{u}_m\hat{v}_n +\nu k^2\hat{u}_k =0, \qquad \text{with} \qquad \hat{v}_k = ik \hat{u}_k.$$

(b) Pseudo-Spectral approximation (aliased): $$ \frac{\text{d}\hat{u}_k}{\text{d}t}+\sum_{m+n=k} \hat{u}_m\hat{v}_n + \sum_{m+n=k\pm N}\hat{u}_m\hat{v}_n + \nu k^2\hat{u}_k =0, \qquad \text{with} \qquad \hat{v}_k = ik \hat{u}_k.$$

(c) Collocation approximation (approximation takes place in physical space only): $$ \frac{\partial u^N}{\partial t}+ u^N v^N - \nu \frac{\partial^2 u^N}{\partial x^2}\bigg|_{x=x_j} =0, \quad j=0,..., N, \qquad \text{with} \qquad v^N = \frac{\partial u^N}{\partial x}.$$

---

The pseudo-spectral approximation uses a fully aliased transform method to evaluate the convolution sum. Note that $ \sum_{m+n=k\pm N}\hat{u}_m\hat{v}_n $ is the aliasing error. The collocation approximation is identical to the pseudo-spectral approximation in the sense that they yield the same solution (except for round-off errors). The same equivalence holds if we would approximate the divergence form of Burgers $(2)$ or even more complicated systems such as incompressible Navier-Stokes equations.