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

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?

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.

• You mentioned that sckemes for solving large dense linear systems are well documented in the literature, could you please give some literature recommendations? – user123 Nov 8 '16 at 5:55
• SIAM SC papers are a treasure trove of information on this. LU, QR and Cholesky decompositions have been the standard for a long time, so long as you have no information about the underlying matrix structure. Implementations of these are found in LAPACK/BLAS. Krylov subspace methods are also widely used: GMRES and CG, for example. – Spencer Bryngelson Nov 8 '16 at 18:44

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.

• Is the traditional iterative methods such as SOR, conjugate gradient or multigrid still work well for this kind of less sparse system? Could you please give a brief comparison w.r.t the problem I referred? – user123 Nov 8 '16 at 6:00
• CG does work well, but the issue is that you need to have a good preconditioner. This is not hard to do if you use discontinuous spectral methods (see the papers by Warburton and Hesthaven), but it becomes more cumbersome for continuous spectral methods. I don't have any first hand experience with things such as SOR applied to spectral method, but would imagine that these fixed point methods are really not good solvers for spectral methods; I would suspect that they are also not good preconditioners. – Wolfgang Bangerth Nov 8 '16 at 21:13
• As I mentioned, the papers by Tim Warburton and Jan Hesthaven will teach you a lot about these sorts of things. I would also recommend looking at the papers by Martin Kronbichler. – Wolfgang Bangerth Nov 8 '16 at 21:14