How is the dense system usually dealt with in spectral method?

Question

Unlike finite element (FEM) or finite difference methods (FDM), where the original PDE is transformed into a sparse linear system, spectral methods return a dense linear system. For a large system, it's not practical to use the a direct elimination method, so a natural choice is iterative methods. As far as I know, most are targeted to sparse matrices.

So, how this problem treated in the scientific computing world? In the case of Chebyshev polynomials, an FFT might give some help. What about other polynomials, such as Legendre or Lobatto? Is this difficulty the reason why this method is not so widely used as FEM or FDM?

Answer 1

There are a few important points to be made here.

Accuracy

Spectral methods exhibit, appropriately, spectral (super-algebraic) convergence. As a result, much smaller matrix sizes are required to produce the same accuracy as a FEM or FDM. Of course, the specific size and relative advantages depend on the specific scheme chosen and stability and resolution desired. An excellent example of this is shown below (S.W. Armeld et al. JCAM, 2009), where the eigenspectra of the Orr-Sommerfeld operator is computed and accuracy compared between a standard FDM and various spectral discretizations. Indeed this is a principal example of how spectral methods became ubiquitous in computation science, illustrated first by the hallmark paper Orszag JFM, 1971.

Matrix structure and solve

As you rightly point out, spectral methods yield dense and typically unstructured difference matrices. The expense of computing the eigensystem or a linear solve of such a matrix is no doubt more expensive than a symmetric, sparse system from a FDM scheme, were they both the same size. The exact expense, again, depends on the specific scheme and solution method. Methods for solving both such matrix systems are extensively studied and well documented in the literature.

Conditioning

This is a subtle but very important point. Spectral methods are typically horribly conditioned. Specifically, polynomial methods have condition number

$$\kappa \equiv \frac{\sigma_{max}}{\sigma_{min}} = \mathcal{O}(N^{2p})$$ where $p$ is the order of the highest derivative and $N$ is the number of collocation points. Accurately evaluating solutions to poorly conditioned matrices requires care and sometimes extended precision, making such operations more costly. This also places a limit on the accuracy you are able to obtain for a given precision (conveniently you can see this behavior in the above figure, also).

Interestingly, the above expression for $\kappa$ has yielded methods designed to split operators such that you have two coupled order $p/2$ equations, as opposed to a single order $p$ equation. There are other specialty methods to finagle your problem to be less computationally horrible, though I won't document them all here.

Answer 2

Even with spectral methods, the matrix is still sparse: if you fix the polynomial degree (even if to a large number) and refine the mesh, the number of nonzeros per row of the matrix stays constant. That's the definition of a sparse matrix. It's also independent of the kind of basis you use on each cell -- Legendre, Gauss-Lobatto, Chebyshev, all only couple locally on a cell and consequently lead to a sparse matrix.

What best to use is a separate question. First, the usual iterative solvers that work for low-order polynomial degrees still work. Matrix-vector products just cost more because the matrix is less sparse. In return, you can efficiently implement spectral methods matrix-free, i.e., where the matrix is not actually formed and only its action is implemented when multiplying with a vector.

A bigger question is how best to precondition the linear system. I like the experience to comment on that in sufficient detail.