Keywords: pseudospectral differentiation; Chebyshev-Gauss-Lobatto grid; Legendre grid; equidistant grid; accuracy; floating-point arithmetic.
$\mathbf{2010~MSC:~65D25~(Primary),~65M70;~65N35~(Secondary).~}$

Fundamental results from approximation theory refers to the best uniform approximation of a smooth function by polynomials and the uniform approximation of a smooth $2\pi$ periodic function by trigonometric polynomials. These results are reviewed for instance in the monograph [1], Ch. 3. Thus, the Chebyshev equioscillation theorem states that a best approximation is unique for important classes of approximating functions but the set of nodes on which this approximation is realized is not necessarily unique. On the other hand, the spectral collocation methods essentially depend on the set of nodes on which the differential equation is collocated. In order to implement these methods one needs very accurate differentiation matrices of various orders.

Consequently, the accuracy at which the pseudospectral (collocation) differentiation matrices operate is of utmost importance in numerical analysis. In this note we will address the issue of the dependence of the accuracy of differentiation process on the node distribution in a grid covering a finite interval. We will consider an equispaced grid, an arbitrary one and then Legendre and Chebyshev grids. For non algebraic polynomials we will consider the Fourier interpolate.

It is well established that spectral collocation methods for solving differential equations is based on weighted interpolants of the form

where $u_{j}:=u\left(x_{j}\right),u:[a,b]\rightarrow\mathbb{R}$, the set of interpolating functions $\phi_{j}(x),j=1,\ldots,N$ satisfy $\phi_{j}\left(x_{k}\right):=\delta_{jk}$ (the Kronecker delta), the set of nodes $x_{j},j=1,\ldots,N$ in $[a,b]\subset\mathbb{R}$ is distinct but otherwise arbitrary and the weight $\alpha(x)$ is an arbitrary continuously differentiable positive function (see for instance our contribution [3]). Typically, $u(x)$ is the solution of an initial/boundary value problem.

The baricentric form of interpolation polynomial writes (see for instance the seminal paper [7])

where $w_{j}^{-1}:=\prod_{m=1,m\neq j}^{N}\left(x_{j}-x_{m}\right)$. This means that $p_{N-1}(x)$ defined above is an interpolant of the function $u(x)$ in the sense that

The collocation derivative operators are generated by taking various order derivatives of (1) and evaluating them at nodes $x_{k},k=1,\ldots,N$, i.e.,

Consequently, the $l-$ th order differentiation matrices associated to this operator are computed by

where $\phi_{j}(x)$ are given by Lagrange’s formula

The approximation theory dictates that the set of nodes $x_{j},j=1,\ldots,N$ cannot be just any set of nodes. The main aim of this note is to make this statement more clear. Thus we use the MATLAB code poldif.m from [7] in order to perform the differentiation in (2).

In order to quantify the performances of every set of nodes we compute two specific parameters, namely:

• •

  the norm of the error in approximating the zero vector, i.e., $\left\|D^{(1)}\cdot\mathbf{1}_{N}\right\|$ where $\mathbf{1}_{N}:=\operatorname{ones}(N,1);$

  Report issue for preceding element
• •

  the $\operatorname{rank}\left(D^{(1)}\right)$ for various values of approximation parameter $N$.

  Report issue for preceding element

Instead of the first parameter we could use another one reflecting the fact that the matrix $D^{(1)}$ has to satisfy

i.e., the derivative of a constant vanishes.

Let’s consider a set of equidistant nodes

where $N$ takes in turn the value from the first row of Table 3.1. Thus the interval [0,1] is successively divided in $N-1$ subintervals each of length $1/(N-1)$.

Along with the equidistant grid (3) we consider first an arbitrary one, namely

It is fairly clear that this new grid is a non uniform one with nodes clustering to the left margin 0 . The extent at which the differentiation matrices perform on both grids is depicted in Fig. 1. The error in approximating the zero vector takes huge values. The collocation derivatives of the hat function and of a highly oscillatory but fairly smooth function, i.e., $\exp(\sin(nx))$, with $n=10$ on 20 equispaced nodes are depicted in Fig. 2. They show large oscillations which cluster to the ends of the interval $[-1,1]$. They are the direct consequence of the numerical instability of the differentiation process.

It is well known that given $N$ nodes each of the differentiating matrices should be rank $N-1$, a differentiating has the constant vector as its null space. Thus, it is important to point out that this equispaced approach works accurately for a very small number of nodes. As differentiating matrix $D^{(1)}$ should be rank one deficient, this simple test shows that even for a very rough approximation $D^{(1)}$ has additional nullspaces. Moreover, in case of quadratically spaced grid $\widetilde{x}_{j}$, the differentiation matrices

become much more degenerated. Table 3.1 shows that the situation dramatically deteriorates as $N$ increases.

In this section we will analyze some well known grids such as Legendre, Chebyshev and Fourier applied to a finite interval. First, let’s reconsider two illustrative examples, the Chebyshev and Legendre derivatives of the hat function $h(x):=\max(0,1-abs(x))$ and of the smooth but fairly oscillatory function $\exp(\sin(nx)),n>$ 2.

perform satisfactorily well. Our numerical experiments have showed that all differentiation matrices keep a rank of order $N-1$. Most notably, for $N$ larger than 700 the Legendre differentiation process rapidly becomes unstable.

The best situation with polynomial differentiation is encountered when Chebyshev nodes of second kind

or equivalently Chebyshev-Gauss-Lobatto quadrature nodes (see for instance [5] or our contribution [4] p. 11) are used. The corresponding differentiation matrix has the entries (see for instance [2], p. 69)

The above differences $\left(x_{k}-x_{j}\right)$ may be subject to floating-point cancellation errors for large $N$. Various tricks based on simple trigonometric identities have been used in [7] in order to avoid such errors in floating-point arithmetic. Thus, the MATLAB code chebdif.m has been fairly stable algorithm in computing these matrices.

The upper curve in Fig 5 correspond to this situation.

Beyond this polynomial differentiation we will pay a particular attention to the Fourier differentiation matrices (see [7]). The Fourier interpolate reads

where

Using the baricentric form of the interpolate (see [6] Sect. 13.6) the MATLAB code fourdif.m from [7] provides the Fourier differentiation matrices on the nodes

Their performances are the best as it is apparent from Fig. 5 (see the lower curve). For $N=O\left(2^{6}\right)$ it is of utmost importance to underline that Fourier and even Chebyshev differentiation matrices work fairly close to the machine precision. Excellent approximations are also attained when the cut off parameter ranges up to $2^{11}$. It is also important to notice that all differentiation matrices conserve a correct rank. It is a practical illustration of the spectral accuracy.

In ([7]) p. 478 the authors state that no general error analysis applicable to differentiation process on arbitrary grids has been undertaken. We hope that the present note fills this gap at least partially. Fourier and Chebyshev differentiation matrices have proved to be fairly reliable. This facts explain to some extent the success of the collocation spectral methods based on them.

[1] Atkinson, K., Han, W., Theoretical Numerical Analysis. A Functional Analysis Framework, Third Ed., Springer 2009.
[2] Canuto, C., Hussaini, M.Y., Quarteroni, A., Zang, T.A., Spectral Methods in Fluid Dynamics, Springer Verlag, Berlin, 1988.
[3] Gheorghiu, C.I., Spectral Methods for Differential Problems, Casa Cartii de Stiinta Publishing House, Cluj-Napoca, 2007.
[4] Gheorghiu, C.I., Spectral Methods for Non-Standard Eigenvalue Problems. Fluid and Structural Mechanics and Beyond, Springer Briefs in Mathematics, 2014.
[5] Gottlieb, D., Orszag, S.A., Numerical Analysis of Spectral Methods: Theory and Applications, SIAM Philadelphia, 1977.
[6] Henrici, P., Applied and Computational Complex Analysis: Discrete Fourier Analysis-Cauchy Integrals-Construction of Conformal Maps-Univalent Functions, Vol. 3 John Wiley and Sons, Inc., New York, 1986.
[7] Weideman, J.A.C., Reddy, S.C. A MATLAB Differentiation Matrix Suite, ACM Trans. on Math. Software, 26, 465-519 (2000).