A Mathematical Journal Vol. 6, No 4, (1-32). December 2004. Spectral Methods for Partial Differential Equations Bruno Costa 1 ABSTRACT In this article we present the essential aspects of spectral methods and their applications to the numerical solution of Partial Differential Equations. Starting from the fundamental ideas, we go further on building auxilliary techniques, as the treating of boundary conditions, and presenting accessory tools like mapping and filtering, finishing with a complete algorithm to solve a classical problem of fluid dynamics: The flow through a circular obstacle. We also present a short comparison with Finite Differences, showing the superior efficiency of spectral methods in problems with smooth solutions. Several equations like the wave equation in one spatial dimension, Burgers in a 2D domain and a simple multidomain setting for the Navier-Stokes 2D are solved numerically and the results are presented at the applications sections. We end with a brief presentation on the software PseudoPack2000 and a quick dis- cussion on the relevant literature. 1 Introduction Spectral methods have been used extensively during the last decades for the numerical solution of partial differential equations (PDE) due to their bigger accuracy when compared to Finite Differences (FD) and Finite Elements (FE) methods. The rate of convergence of spectral approximations depends only on the smoothness of the solution, yielding the ability to achieve high precision with a small number of data. This fact is known in literature as ”spectral accuracy”. In FD and FE, the order of 1This work was supported by CNPq grant # 300315/98–8. 2 Bruno Costa 6, 4(2004) accuracy depends on some fixed negative power of N , the number of gridpoints being used, even when only analytic functions are involved. The expression spectral methods has different meanings for several subareas of Mathematics, like Functional Analysis and Signal Processing. In this article, spec- tral methods has the meaning of a high accuracy numerical method to solve Partial Differential Equations. The numerical solution is expressed as a finite expansion of some set of basis functions. When the PDE is written in terms of the coefficients of this expansion, the method is known as a Galerkin spectral method. Spectral colloca- tion methods, also known as pseudospectral methods, is another subclass of spectral methods and are similar to Finite Differences methods due to direct use of a set of gridpoints, which are called ¨collocation points¨. A third class are the Tau spectral methods. These methods are similar to the Galerkin spectral methods, however the expanding basis is not obliged to satisfy boundary conditions, requiring extra equa- tions. In this article, we will mainly concentrate on collocation methods. The reader interested on Galerkin and Tau methods must check the works of Canuto et al [8] and Gotllieb and Orszag [20]. In the past several years, the activity on both theory and application of spectral methods have been concentrated on collocation spectral methods. One of the reasons is that collocation methods deals with nonlinear terms more easily than Galerkin methods. The nonlinear terms are treated on a FD manner, via gridpoints values multiplication. The underlying idea of a collocation spectral method is to approximate the unknown solution in the entire computational domain by an interpolating high order polynomial at the collocation points. The spatial derivatives of the solution are approximated by the derivatives of the polynomial and the time derivative, if it exists, is solved through classical finite differences schemes. Periodic and non-periodic problems are respectively treated with trigonometric and algebraic polynomials. Some of the methods commonly used in the literature are the Fourier collocation methods for periodic domains and the Jacobi polynomials for non-periodic domains, with the Chebyshev and Legendre polynomials as special cases. Pseudospectral methods have a wide range of applications going from 3-D seismic wave propagation to turbulence, combustion, non-linear optics and passing through aero-acoustics and electromagnetics. The material exposed in this article intends to enable the interested reader to solve the ubiquitous Navier-Stokes equations, whose solution is exemplified in the case of the two-dimensional flow through an obstacle. Thus, the necessary main ingredients as mapping and boundary conditions are treated in some extent along with several basic examples. This article is divided as follows: In Section 2, we state the fundamentals of spec- tral methods by describing the Galerkin spectral methods for linear and nonlinear problems. The collocation version of spectral methods comes in Section 3 where we introduce polinomial interpolating algorithms. In Section 4, we discuss mapping tech- niques for spectral methods in one and two dimensions in space. The subjects of fil- tering and boundary conditions are dealt with in Section 5, together with applications to basic equations of fluid dynamics. We finish the article with a quick description of the Library PseudoPack and a brief discussion on the relevant literature. 6, 4(2004) Spectral Methods for Partial Differential Equations 3 2 Preliminaries We start this section showing a direct comparison between spectral methods and Finite Diferences. The main idea to be passed is the following: Spectral Methods achieve a greater precision with a smaller number of points than Finite Difference methods. The table below, extracted from [8] shows the error decay of approximating the first derivative of esin x through a second and a fourth order FD method and a Fourier spectral collocation method: N FD 2nd order FD 4th order Fourier Spectral 16 3.4020e-001 1.3844e-001 4.3179e-003 32 9.3589e-002 1.5750e-002 1.7619e-007 64 2.5833e-002 1.1318e-003 2.3870e-014 128 6.5118e-003 7.5900e-005 7.2054e-014 The graph below in Figure (1) shows the decay of the error with varying number of points. Figure 1: Error comparison between Spectral and Finite Differences schemes. In [35] one can find quick Matlab codes to generate graphics like the one above. 4 Bruno Costa 6, 4(2004) 2.1 Linear periodic problems In spectral methods the solution of a PDE is approximated by a finite linear combi- nation of a chosen set of orthogonal functions { ϕj (x) }∞ j=0 as fN (x) = N∑ k=0 f̂kϕk(x), (1a) where the f̂k’s are called the spectral coefficients of f with respect to the basis {ϕk}; f̂k is a measure of how much of f is in the ”direction ” of ϕk and can be evaluated through an integral like f̂k = ∫ f (x) ϕk(x)w (x) dx, (2) where w (x) is a weight function associated to the basis {ϕk}. For instance, if we consider periodic functions in the interval [0, 2π], one such a basis is the Fourier basis formed by sines and cosines as {cos kx, sin kx}∞k=0 , which in its concise complex form becomes {eikx}∞k=−∞. If f satisfies some conditions, like continuity and bounded variation, one can prove that fN given by fN (x) = N∑ k=0 f̂ke ikx with f̂k = ∫ f (x) eikxdx, converges to f uniformly as N increases. The derivative of f is approximated by f′N which is given by f′N (x) = N∑ k=0 f̂kϕ ′ k(x). (3) f′N is named the Galerkin spectral derivative of f . Once we know the specifics about the derivatives ϕ′k (a recurrence formulae, in general), we can work out formula (3) and find out the coefficients f̂′k to write f ′ N (x) = ∑ f̂′kϕk(x). For instance, let f be an even periodic function in [0, 2π]. An appropriate basis to approximate f is the Fourier Basis {cos kx}∞k=0. If we write fN (x) = N∑ k=0 f̂k cos kx, we easily see that an approximation to the second derivative of f is given by: f′′N (x) = N∑ k=0 ( −k2f̂k ) cos kx, meaning that once we know the spectral coefficients f̂k of f we obtain an approxi- mation to f′′ by multiplying each fk by (−k2). In this way, double differentiation in spectral space is equal to multiplication by (−k2). The case of the first derivative is a little bit more ”complex”, for we need to use the basis {eikx}, since the first derivative of an even function is odd. 6, 4(2004) Spectral Methods for Partial Differential Equations 5 With the above technique to compute second derivatives, let us solve the heat equation below ut − νuxx = 0, x ∈ [−π, π], t > 0, (4) u(x, 0) = e−x 2 , x ∈ [−π, π], where u is the temperature and ν is the difusion coefficient. Using the approximation uN (x) = N∑ k=0 ûk cos kx (5) and substituting into the above equation one obtains N∑ k=0 ( ∂ûk ∂t + νk2ûk ) cos kx = 0. (6) Since the functions {cos kx}Nk=0 are orthogonal, which means that ∫ π −π cos kx cos lxdx = δkl, we can reduce equation (6) to the system of N + 1 ordinary diferential equations ∂ûk ∂t + νk2ûk = 0, for k = 0, .., N, (7) ûk (0) = (u0)k , for k = 0, .., N. The example above shows the essential steps to solve a PDE with spectral methods. First, is necessary to choose an appropriate basis for the problem to be solved, for instance, periodic problems are solved with the Fourier basis and nonperiodic ones with a polynomial basis. Next, we need the formula relating the spectral coefficients of f with those of its derivative. In this way, a system of ODEs analogous to the one in (7) can be formed and solved. Note that one does not work with the physical values of the solution, but with its set of spectral coefficients. In fact, the physical values will be necessary to approximate the integrals (2) through quadrature formulae, as we shall see later. We now show why the error |f − fN| decays so fast with the increase of number of points, as seen in Figure 1 above. If f is such that we can substitute it by its expansion series ∞∑ k=0 f̂kϕk(x), the error in approximating f by fN can be measured by the size of the tail of the above series, given by εN = |f − fN| = ∣∣∣∣ ∞∑ k>N f̂kϕk(x) ∣∣∣∣. Considering the Fourier basis { eikx } , we may write f̂k = 1 2π ∫ 2π 0 f (x) e−ikxdx = 1 2πik [f (2π) − f (0)] + 1 2πik ∫ 2π 0 f′ (x) e−ikxdx. 6 Bruno Costa 6, 4(2004) Therefore, we can see that ∣∣∣f̂k∣∣∣ = O (1k) and that εN = O ( 1N ) . However, if f̂ is periodic, the boundary term vanishes and we can integrate by parts once again to obtain ∣∣∣f̂k∣∣∣ = O ( 1k2 ) or εN = O ( 1N 2 ). Further iteration of this argument shows that the decay of the spectral coefficients of f depend on the diferentiability of f , together with the periodicity of its derivatives. In particular, if f is C∞, we can iterate the integration by parts above indefinitely to affirm that the decay of the spectral coefficients is faster than any negative power of N, or that εN = O ( e−cN ) , for some c > 0. This is what is known as exponential or spectral convergence and it is this property that has given its name to spectral methods. 2.2 Nonlinear Problems and Pseudospectral Methods Mounting the system of ODEs as in (7) gets an extra degree of complexity in the presence of nonlinear terms in the PDE. For instance, consider the same heat equation above with a diffusion coefficient that now depends on the temperature itself, say, ν (u) = cu, where c is a constant. Thus, the equation to be considered is ut −cuuxx = 0. This time, let us use the complex form of the Fourier expansion and substitute the truncated series (5) into the equation to obtain: N∑ k=0 ∂ûk ∂t eikx + ( c N∑ l=0 ûle ilx ) N∑ k=0 k2ûke ikx = N∑ n=0 ∂ûk ∂t einx + c N∑ n=0 N∑ k+l=n ( k2ûlûk ) einx = N∑ n=0 ∂ûk ∂t einx + c N∑ n=0 ŵne inx = 0. The computation of the coefficients ŵn takes O ( N 2 ) operations, which is too expen- sive if compared to the O (N ) cost of a finite difference algorithm to treat the same nonlinear term. In 3 dimensions, this cost rises up to O ( N 4 ) in the spectral method and O ( N 3 ) for finite diferences. It was this single fact that historically kept spectral methods from being used on applications involving a large number of data. This situation changed when trans- form methods between the spectral coefficients and the physical values were devel- oped independently by Orszag() and Eliasen, Machenhauer and Rasmussen (). These tranforms, known by the general name of Fast Fourier Transforms (FFT), allow the computation of the coefficients of the nonlinear term to be performed in O (N log N ) operations. The general idea of the transform technique to compute the spectral coefficients of a general nonlinear term as f (x) g (x) is as follows: Start by computing the physical 6, 4(2004) Spectral Methods for Partial Differential Equations 7 values of f and g through FFTs (O (N log N ) operations), perform the pointwise product in physical space (O (N ) operations) and go back to spectral space, applying the inverse transformation to generate the coeficients ˆ(f g)k (O (N log N )). Orszag (1971) termed the conjugation of the Galerkin spectral method with the transform technique to compute nonlinear terms as a pseudospectral method. This term is also used in the literature to denote collocation methods, which we will introduce in the next section. This happens, because depending on the way the collocation method is designed, it becomes equivalent to a pseudospectral method. For the interested reader, we advance that if the collocation points are chosen as nodes of the quadrature formulae used in the pseudospectral method, than this equivalence occurs (see Canuto and Hesthaven). 2.3 Non-periodic problems For non-periodic problems, the sets of functions choosen as bases are solutions of certain differential equations, the Sturm-Liouville problems. For instance, the set of Chebyshev polynomials defined as Tk (x) = cos (k arccos (x)) are solutions to the following differential equation in the interval [−1, 1]: − (√ 1 − x2T ′k (x) )′ = k2 √ 1 − x2 Tk (x) . (8) It is easy to prove that these functions are indeed polynomials, since T0 ≡ 1, T1 (x) = x and by using trigonometric identities one can arrive to the recurrence relation Tk+1 = 2xTk − Tk−1. The expansion assumes the form fN (x) = N∑ k=0 f̂kTk (x) with f̂k = ∫ 1 −1 f (x) Tk (x) w (x) dx and w (x) = ( 1 − x2 )− 12 . Note that, after the change of variable x = cos (θ), the Chebyshev polynomials become cosine functions, enabling many theoretical results to be transported from the Fourier system, including spectral convergence. If a function f is expanded on a Chebyshev series as f = ∞∑ k=0 f̂kTk, then its derivative can also be expressed on a Chebyshev series as f′ = ∞∑ k=0 f̂ (1) k Tk, whose coefficients f̂ (1)k are computed with the help of the recurrence relation 2Tk (x) = 1 k + 1 T ′k+1 (x) − 1 k − 1 T ′k−1 (x) , k ≥ 1. (9) For instance, when considering an approximation fN as above, we have from (9) that 2kf̂ (1)k = ck−1f̂ (1) k−1 − f̂ (1) k+1, for k ≥ 1, with ck = 2, for j = 0, N and 1 otherwise. We 8 Bruno Costa 6, 4(2004) also have f̂ (1)k = 0, for k ≥ N , and the remaining coefficients can be computed in decreasing order by 2kf̂ (1)k = ck−1f̂ (1) k−1 − f̂ (1) k+1, for 0 ≤ k ≤ N − 1. This relation can be easily generalized to obtain derivatives of higher order (see [8]). Another property inherited from the Fourier system is the possibility of applying fast transforms between spectral and physical spaces. This makes the Chebyshev polynomials the preferred system among the many possibilities of nonperiodic bases. Much more information on the Chebyshev polynomials can be found at the books of Fox and Parker (1968) and Rivlin (1974) . Another set of great use is the set of Legendre polynomials, defined as the solutions of the Sturm-Liouville problem − (( 1 − x2 ) L′k (x) )′ = k (k + 1) Lk (x) . (10) Although there is no fast transform for these polynomials, the spectral coefficients of the expansion f = ∞∑ k=0 f̂kLk are defined as f̂k = ∫ 1 −1 f (x) Lk (x) dx, where the weight function w (x) = 1 makes much easier to prove theorems for this set than for the Chebyshev polynomials. Here, the book of reference is the one by Szegö (1939) . Note that both functions in the left hand side of equations (8) and (10) vanish at ±1. It is this property that allows the vanishing of the boundary terms on the integration by parts process, yielding spectral convergence to the Chebyshev and Legendre bases. Sturm-Liouville problems with this property are named singular. Infinite domains can also be treated by polynomials. Hermite and Laguerre poly- nomials are the common bases used for infinite and doubly infinite domains respec- tively. It is also common to map infinite domains to finite ones and use the Chebyshev polynomials (see [8]). 3 Polinomial Interpolation In the collocation version of spectral methods, one requires the numerical approxima- tion to be exact on a set of pre-defined points, the collocation points. The general idea is to interpolate the solution on the collocation points and approximate its derivatives by the derivatives of the interpolating polinomial. As we shall see below, collocation spectral methods treat with equal cost problems with linear, variable coefficients or non-linear terms. All products are performed pointwise at the physical space and this makes collocation schemes very similar to a high order finite differences scheme (see [18]). 6, 4(2004) Spectral Methods for Partial Differential Equations 9 For the sake of clarity, let us start with the step of computing the numerical approximation of the derivative fx of a real function f , from which we only know its values {f (xj )} N j=0 at a set of collocation points {xj} N j=0 . We need to find a polynomial of degree N that interpolates f at these N + 1 points. Note that if the degree is less than N, such polynomial might not exist (not all sets of 3 points are colinear). It is also easy to show that we do not have uniqueness (monic uniqueness) if the degree is greater than N . A general formula for such a polynomial depending only on the values {f (xj )} N j=0 is given by (IN f ) (x) = N∑ j=0 f (xj )gj (x), (11) where gj (x) is a Nth-degree polynomial satisfying gj (xk) = δjk = { 1 if j = k, 0, otherwise. The polynomials {gj (x)} N j=0 are called cardinal functions. It is easily seen that IN f (x) is a Nth-degree polynomial interpolating f , since it is a linear combination of the Nth-degree polynomials gj (x), and it satisfies (IN f ) (xk) = f (xk), k = 0, ..., N. (12) In this way, given a set of collocation points, we only need to find the N + 1 Nth- degree polynomials gj (x), which can be easily computed by means of Lagrangian interpolation through the formula gj (x) = P (x) P ′(xj )(x − xj ) , where P (x) = N∏ j=0 (x − xj ). The derivative of the interpolating polynomial is simply given by the (N − 1)th- degree polynomial (IN f ) ′ (x) = N∑ j=0 f (xj )g ′ j (x), (13) providing an easy way to obtain the approximated values (IN f ) ′ (xi): we only need to multiply the vector composed of the function values {f (xj )} N j=0 by the differen- tiation matrix D, whose elements are given by Di,j = g′j (xi).In this way, numerical differentiation through collocation spectral methods can be accomplished through matrix-vector multiplication. Once the collocation points are provided and the differentiation matrix is obtained, the solution of a PDE is very close. For instance, given an initial temperature distribu- tion ~u0 = [ u00, u 0 1, ..., u 0 N ] , with u0j = u (xj , t = 0), one can obtain numerical approxi- mations to future temperature profiles through a first order finite difference temporal 10 Bruno Costa 6, 4(2004) discretization of the heat equation ut − νuxx = 0, conjugated with the collocation spatial derivative: ~un+1 − ~un ∆t − νDxx~un = 0, or ~un+1 = ~un + ∆tDxx~un, (14) where Dxx is the second derivative pseudospectral matrix that can be obtained by squaring the differentiation matrix D. If we suppose, as in Section 3.2, that the diffusion coefficient is given by ν (u) = cu. Thus, scheme (14) becomes ~un+1 = ~un + c∆t~un · Dxx~un, where the product ′·′means pointwise product. Note that at equation (14), all values of u, including boundary values, are com- puted simultaneously. In Section 5, we will see that some of the boundary values must be discarded in order to respect the physics of the problem being solved. 3.1 The Chebyshev Collocation Method We will now present the Chebyshev collocation method starting by its most used set of collocation points, the Chebyshev-Gauss-Lobatto points: xi = cos ( πi N ) i = 0, ..., N, (15) which take this name because they are the nodes of the Gauss-Lobatto quadrature formula for the Chebyshev polynomials (see [8]). They are also the extrema of the N–th order Chebyshev polynomial TN (x) and can be seen as the projection on the x−axis of the equispaced set of points θi = πiN , in the semicircle, as in Figure (2) below Figure 2: The Chebyshev points as projections on the x-axis of an equispaced set of points in the semicircle. Note that points are concentrated close to the boundaries. 6, 4(2004) Spectral Methods for Partial Differential Equations 11 The N–th order cardinal function gj (x) is given in this case by gj (x) = (−1)j+1(1 − x2)T ′N (xj ) cj N 2(x − xj ) cj = 1 + δj,0 + δj,N , (16) and the differentiation matrix in closed form is given by Dkj = ck cj (−1)j+k xk − xj , k 6= j (17) Dkk = −1 2 xk 1 − x2k , k 6= 0, N (18) D00 = −DN N = 2N 2 + 1 6 . (19) We solved the linear heat equation (4) in the interval [−1, 1], using the Chebyshev collocation method, see Figure (3). Since we have imposed temperature zero at both boundaries, we see that the initial condition e−cx 2 converges to a uniform zero tem- perature on the whole interval, as it should be if we heat a metal bar in the center and immerse its ends on ice. Figure 3: Solution of the Heat equation by the Chebyshev pseudospectral method. We end these two introductory sections on spectral methods pointing out that the value of the derivative of a function f at a point xj is a linear combination of all values of f at the remaining points, as it can be seen from (13). This is why spectral methods are called global methods and also means that the differentiation matrices in collocation spectral methods are full, making more expensive in general the numerical differentiation when compared to Finite Differences or Finite elements, where derivatives are computed at each point using only a small number of adjacent points, a characteristic of local methods. 12 Bruno Costa 6, 4(2004) 4 Mapping The Chebyshev and Legendre bases of polynomials are the most used sets of functions for the spectral expansions of solutions of Partial Differential Equations. However, being eigenfunctions of Sturm-Liouville problems, the number of domains which can be modelled with them is limited. In this section we will see how to make use of mappings between two-dimensional domains in order to overcome this barrier and still be able to use valuable numerical tools such as the Fast Fourier Transform. Mappings also have different applications other than only conforming physical to computational domains. For instance, concentration of points on a region of large gradients might improve resolution. On the other hand, increasing the minimum distance of a set of collocation points increases the stability of temporal integration schemes. Below, we deal with such questions when developing a general framework for 1D mappings, aiming the integration of mappings to differentiation operators in order to solve equations on more complex domains without including the metrics terms into the equations. 4.1 Grid Transformation in 1D Consider a set of collocation points {yj} j = 0, ..N, and a mapped set xj = g (yj ) . Given a function u, with values u (xj ) at the new points {xj},the derivative of u can be evaluated as du dx = du dy dy dx , (20) where dy dx = 1 g′(y) . In matricial form, (20) can be rewritten as ~u′ = M D~u, (21) where M is the diagonal matrix with elements Mjj = 1g′(yj ) and D is the differentiation matrix with respect to the set of points {yj}. Note that we need only O(N ) more operations to apply M , after the differentiation process D~u is concluded. Since fast transforms are only applicable to the Chebyshev-Gauss-Lobato points (15), we see that the above algorithm keeps the possibility of using the FFT for differentiation on sets of collocation points which are mapped from the Chebyshev ones. If we take the example of the most simple and common 1D mapping, which is used to conform the Chebyshev computational domain [−1, 1] to a physical domain [a, b] of interest: x = g (y) = a + (y + 1) 2 (b − a) . (22) Then, the matrix M is the constant diagonal matrix Mjj = b−a2 . In general, 2D or 3D grids are cross-products of 1D sets of points. For instance, one can simulate a periodic channel by considering a Chebyshev grid in the vertical and a Fourier grid in the horizontal direction. One can also define a grid on the circle by using this same arrangement in the radial and angular directions, respectively. 6, 4(2004) Spectral Methods for Partial Differential Equations 13 This is a natural way to proceed since the derivative operators work in each direction separatly. Thus, characteristics of the 1D distribution are transported to the 2D or 3D grid and some properties of the multidimensional method are determined by individual properties of each onedimensional grid component. A relevant practical issue when numerically solving evolutionary PDEs is the size of the maximum timestep allowed by the numerical method. It is empirically shown that this size is in inverse proportion to the minimum separation between the points, which is directly proportional to the number of points. The Chebyshev points are naturally agglomerated at the boundaries and this causes the timestep to be exceedingly small in comparison to other methods. For instance, Finite Differences schemes allow a O ( N−1 ) time step when solving the first order wave equation, while a O ( N−2 ) size must be respected to obtain stability with spectral methods. In [26], Kosloff and Tal-Ezer proposed the following mapping in order to mitigate this problem: x = g (y, α) = arcsin (αy) arcsin (α) , (23) The mapping above stretches the grid spacing on the boundary pushing the points to the center of the domain, generating a quasi-uniform grid, where the strength of the endpoints separation depends on α ∈ (0, 1). Kosloff and Tal-Ezer claimed that the mapping decreases the spectral radius of the derivative operator from O ( N 2 ) to O (N ) , increasing the allowed time step from O ( N−2 ) to O ( N−1 ) , in the case of the one-way wave equation. However, it can be easily seen that the mapping has singularities at y = ± 1 α , which introduce some approximation error and make necessary some expertise on choosing the correct value of α. See [26] and [17] for detailed theoretical and numerical information on the choice of α, which might depend on N. Besides the increase of the timestep, the mapping (23) also increases the resolution of the high modes because the maximum grid spacing ∆xmax is diminished. Resolution of a mode is attained when exponential decay is observed on the spectral coefficients. Thus, the improvement of resolution can be measured by looking to the number of points necessary to achieve this order of decay for a fixed frequency. Figure(4) shows the size of the spectral coefficients of cos (mx) for varying m and N , where we see that resolution is attained exactly at the point where N > mg′ (0, α) = mα arcsin α , where it should be noticed that ∆xmax = g′ (0, α) ∆ymax increases when α → 1. For a more up to date discussion on the Tal-Ezer mapping see [4],[9], [13],[30],[31] and [28]. 14 Bruno Costa 6, 4(2004) Figure 4: Decay of the spectral coefficients of cos(mx) for several values of m and number of gridpoints N showing the dependence of spatial resolution on the parameter α. 4.2 Transfinite Mapping In this section we define a simple and useful bidimensional mapping which consists of mapping the square [0, 1]2 onto any quadrilateral with curved boundaries. The Transfinite mapping is defined by four functions from the segment [0, 1] onto the curves defining the sides of the quadrilateral. Let ~πi (ξ) = (xi (ξ) , yi (ξ)) , ξ ∈ [0, 1] , i = 1, 3 and ~πj (η) = (xj (η) , yj (η)) , η ∈ [0, 1] , j = 2, 4 be the functions defining the sides of the quadrilateral, then the transfinite mapping is expressed by ~Ω (ξ, η) = 1 − η 2 ~π1 (ξ) + 1 + η 2 ~π3 (ξ) (24) + 1 + ξ 2 ( ~π2 (η) − 1 + η 2 ~π2 (1) − 1 − η 2 ~π2 (−1) ) + 1 + ξ 2 ( ~π4 (η) − 1 + η 2 ~π4 (1) − 1 − η 2 ~π4 (−1) ) . This mapping was first presented in [19] and in the figures below we show some domains which can be generated. Two very interesting cases are the ones of the circle and the ellipse. Below, in Figure (5), we show a circular grid which is obtained by dividing the boundary of the circle in four curves and mapping them to the interval [0, 1] as above. This is an easy way to avoid the singularities of polar coordinates which appears at the center of the circle due to the vanishing of the radius. 6, 4(2004) Spectral Methods for Partial Differential Equations 15 Figure 5: Circular grid generated by the Transfinite Mapping. Figure (6) below shows the solution of the second-order wave equation utt − c2∆u = 0, on a elliptical domain, generated in a similar way as the circular one above, simulating the travelling of sound waves between the foci of an elipse. Figure 6: Acoustic wave equation simulation on the elliptical grid generated by the transfinite mapping. 4.3 Two-Dimensional Mappings Consider the classical case of the polar coordinates transformation (x (r, θ) , y (r, θ)) 7→ (r cos θ, r sin θ) , (r, θ) ∈ [r1, r2] × [0, 2π] , 16 Bruno Costa 6, 4(2004) where the radial direction is discretized with the Chebyshev collocation points and the angular, with the Fourier points. If we need to solve a problem in the annular domain of Figure (7), this is clearly the transformation to be used. Figure 7: Grid obtained from the polar mapping. For instance, the problem one wants to solve might be the case of a fluid between two concentric and rotating cylinders. Suppose we are only interested in the dynamics on a cross section perpendicular to the axis and, to simplify even more, let us look only to the diffusion of heat on the space between the cylinders. Thus, the equation to be discretized is the two-dimensional heat equation ut − ν (uxx + uyy) = 0. (25) We have two alternatives: The first one is to apply the change of variables to the equation and solve it in the new variables r and θ. For instance, the new equation would become ut − ν ( urr + ur r + uθθ r2 ) = 0. (26) This approach has the drawback of requiring some analytical work that sometimes can get cumbersome if the mapping is a little bit more involved. Take a look at the transfinite mapping (24) above. The second approach is to express only the derivatives in (x, y) of the original equation in terms of the computational variables (r, θ) . For instance, differentiation with respect to x is carried out using the chain rule: ∂u ∂x = ∂u ∂r ∂r ∂x + ∂u ∂θ ∂θ ∂x . The derivatives with respect to r and θ are respectively the Chebyshev and Fourier derivatives, but we still need to spend some work on the computation of the quantities 6, 4(2004) Spectral Methods for Partial Differential Equations 17 ∂r ∂x , ∂r ∂y ∂θ ∂x and ∂θ ∂y . No problem in the polar mapping case, however, look again to the transfinite mapping and see that it can get very troublesome. Thus, it would be handy to have a formula for ∂u ∂x and ∂u ∂y not involving the metrics ∂r ∂x , ∂r ∂y , ∂θ ∂x and ∂θ ∂y , but the inverse metrics ∂x ∂r , ∂x ∂θ , ∂y ∂r and ∂y ∂θ , which can be computed on a straightforward way. We can resort to some very simple linear algebra and write ( ∂u ∂r ∂u ∂θ ) = [ ∂x ∂r ∂y ∂r ∂x ∂θ ∂y ∂θ ]( ∂u ∂x ∂u ∂y ) . We solve the system above for [ ∂u ∂x , ∂u ∂y ] by simply using the straightforward formula for the inverse of a 2 × 2 matrix and obtain:( ∂u ∂x ∂u ∂y ) = 1 J [ ∂y ∂θ −∂y ∂r −∂x ∂θ ∂x ∂r ]( ∂u ∂r ∂u ∂θ ) , (27) where J = ( ∂x ∂r ∂y ∂θ − ∂y ∂r ∂x ∂θ ) is the Jacobian of the mapping. At this point, we would like to advocate a computational procedure which makes simpler the solution of PDEs on mapped domains. Looking at formulae, we see that after computing the inverse metrics and the Jacobian, which can be done at the begining of the computations, we can define the operators ∂ ∂x and ∂ ∂y and code the equation in its original form, with no need to include the metrics. The library PseudoPack contains routines for the computation of the metrics and the Jacobian, which makes the computation of a term like uxx as simpler as Call PS Diff (D t, D r, u, uxx, 2, radx, rady, thx, thy, Jacob) We will see below, in Section 5, that this procedure is of great help on clarifying the coding and implementation of more involved equations as the Navier-Stokes and enables the clear use of more complex differential operators, like the ones of Vector Calculus. 5 Applications 5.1 The Wave Equation The first-order wave equation:{ ut + cux = 0, t > 0, x ∈ [−1, 1], u (x, 0) = u0 (x) , x ∈ [−1, 1], (28) provides us with a typical situation on the study of boundary conditions: The one where some quantity is carried through the domain towards its boundary. It can be seen as the modelling of a travelling wave on a unidimensional medium as a material pulse on a rope or an acoustic wave on a piece of wire. The constant c is the speed of the wave. One can easily check that u (x, t) = u0 (x − ct) (29) 18 Bruno Costa 6, 4(2004) is the solution to equation (28), and if c > 0, it means that information is being propagated to the right, or that values of the solution at a point x depends on values of the function at points to the left of x and at previous times. If one is interested in the propagation of the initial pulse u0 (x) with no other influence than the pure advection phenomenon, no boundary condition must be set at the point x = 1, since the solution values are going to be ”brought” there by the numerical integration itself (remember that spectral differentiation also computes values for the boundary points, see Section 3). The imposition of a boundary condition such as u (1) = 0, as it might be naively proposed, would simulate a numerical barrier on which the travelling pulse would hit and reflect backwards with inverse values. Just think about swinging a rope with one end fixed at a wall. On the other hand, we necessarily have to impose a boundary condition at x = −1 and discard the value computed by the spectral differentiation, otherwise, we would be going against the restriction expressed by equation (29) that information must always come from the left. If we use the Gauss-Lobato set of points {xj}, then x0 = 1 and xN = −1 are boundary points and both are updated at the temporal integration process, since the differentiation matrix includes lines for these points. However, the value of u at −1 has to be thrown out and substituted by a convenient boundary condition, such as u (−1) = 0, (30) at the end of every step (at the end of every stage, when using a Runge-Kutta scheme, see [11]). Remark 5.1 In the case above, the inclusion of the boundary point x = 1 is not essential for the simulation. In these situations one could use instead the Gauss- Radau set of points xj = cos 2πj 2N + 1 , 1 ≤ j ≤ N. In the case of the second-order wave equation, however, both endpoints are needed, since we will have waves moving in both directions, as we shall see below. Remark 5.2 Statement (30) is equivalent to zeroing the last row of the differentiation matrix. This can be generalized to the case of a second order PDE with homogeneous boundary conditions, where one would fill the first and last lines with zeroes. This does not apply if one has nonzero Dirichlet boundary conditions. The typical situation provided by the one-way wave equation (28) will help us to design boundary conditions for problems of greater interest like the Navier-Stokes equations. Equations whose solutions involve propagation of information at a finite speed are called hyperbolic equations and (28) is their very simplest example. In what follows, we will try to reduce more complicated hyperbolic equations to a form which is similar to the one-way wave equation: Ut + AUx = 0, (31) 6, 4(2004) Spectral Methods for Partial Differential Equations 19 where U is a vector with the dependent variables and A is a matrix that does not depend on the derivatives of U , but may depend of U itself. For this reason, we denote (31) a quasilinear system of partial differential equations. Another characteristic of hyperbolic equations in the form of the system (31), is that all the eigenvalues of A are nonzero real numbers, in this way we are able to diagonalize A = ΓΛΓ−1(for the sake of simplicity, let us disregard the case of defective eigenvalues), where Γ and Λ are the matrices of eigenvectors and eigenvalues of A, respectively. Applying the variables transformation V = Γ−1U to equation (31) yields the uncoupled system of equations Vt + ΛVx = 0. (32) Since each equation of the above system is similar to the one-way wave equation, we are able to construct boundary conditions for all the transformed variables. The eigenvalues of A are the velocities with which the new dependent variables in V are being propagated. For practical applications, the diagonalization of A is performed ahead of the temporal integration process and the system (32) is used only to impose the boundary conditions to the transformed variables. The process ends by the computation of the boundary values of the original variables by an inverse transformation. Below, we illustrate this technique with the case of the second order wave equation utt − c2uxx = 0 by reducing it to the first order linear system Ut = AUx, where U = ( u v ) and A = [ 0 c c 0 ] . The diagonalization of the matrix A = ΓΛΓ−1 makes possible to write the system above in its characteristic form:( u + v u − v ) t = c ( u + v u − v ) x , (33) since Γ = [ 1 1 1 −1 ] and Λ = [ c 0 0 −c ] . Defining the characteristic functions R1 = u + v and R2 = u − v, we easily obtain the solution of (33) as R1 (x, t) = R1 (x − ct, 0) and R2 = R2(x + ct, 0). Thus, R1 has to have its right boundary condition imposed and R2, its left one. This reasoning is particularly useful when simulating infinite domains, since they impose no barrier for the travelling quantity to leave the computational domain. Due to its derivation, these are called characteristic boundary conditions. Another name for the characteristic functions above is ”Riemann Quantities ” and they are the main subject of study when designing numerical schemes for hyperbolic conservation laws. 20 Bruno Costa 6, 4(2004) For us, they will be of great use in Section 6 when performing this same characteristic decomposition for the more complex (and more interesting) system of the Navier Stokes equations. They will also be useful for the setting of interface conditions for multi-domain implementation, in Section 6. Figure 8: Second-order wave equation with characteristic boundary conditions. In Figure (8) above, we show the solution of the second-order wave equation using characteristic boundary conditions with an initial pulse at the middle of the domain. This initial pulse is decomposed into two smaller ones travelling to opposite directions. Note that when passing through the boundaries, little wiggles start to appear, but nothing is reflected inside the domain and once the waves are gone, nothing is ”sounding”. 5.2 Nonlinearity: Burgers Equation The viscous Burgers equation in 1D is given by: ut + uux − νuxx = 0, (34) and is the simplest equation where the phenomena of diffusion and nonlinear advection are present. It is therefore very used on testing numerical methods for fluid dynamics and their proposed improvements. Another important feature of the Burgers equation is the development of disconti- nuities, or shocks, in finite time. In Figure (9) below, we used the Fourier collocation method with 64 points in order to show the numerical solution of (34), with a Gaus- sian function centered at the middle of the periodic domain as the initial condition. The advective speed is given by u itself, which means that the initial profile is going to move to the right, however, since the velocity start to decrease at center of the domain, particles departing from the center will colide with slower particles on their right, causing the discontinuity observed in the numerical solution. The value of the diffusion constant ν is the measure of how much the solution is going to behave as the heat equation. Thus, smaller the value of ν, sharper is the discontinuity. In this example, a value of v = 0.001 was used. The occurrence of a discontinuity in the solution poses an extra difficulty for spectral methods. The osccilations emanating from the discontinuity are known as the Gibbs phenomenon and are caused by the fact that a discontinuous solution is being approximated by an osccilatory set of smooth functions. 6, 4(2004) Spectral Methods for Partial Differential Equations 21 Filtering is an effective tool to get rid of the osccilations and recover exponential accuracy away from the discontinuity. The central idea of using filters is to increase the order of convergence away from the discontinuity by attenuating the high order coefficients, which decay only linearly. Care must be taken in order to maintain the important information contained in these coefficients, otherwise we might obtain a strongly smeared function. The Gibbs phenomenon is a relevant issue due to the contamination of the solution with the spurious osccilations. The use of filters is recomended in these situations and in Figure (5.2) below, we see the attenuation of the osccilations after the use of a filter of the exponential type, also paying the price of an extra smoothing of the discontinuity. More information on the several types of filters and their properties can be found in [8]. Figure 9: Burgers Equation in 1D. In the two-dimensional case, the viscous Burgers equation assumes the form:   Ut − ν∆U + (U · ∇U ) = f, in (0, 1)2, t in [0, T ] , U = g on ∂Ω, U (x, y, 0) = U0(x, y), ∀ (x, y) ∈ Ω, (35) where U = (u, v). The initial condition U0(x, y) = ( cos (πx) cos (πy) , ( x+y 4 )) and forcing function f (x, y) = (sin(3πx) sin(4πy + 3 sin(16πx) sin(8πy)), sin(4πx) sin(2πy + 2 sin(10πx) sin(16πy))) are chosen in such a way to promote enough dynamics in the fluid, forcing an inter- action between the several frequencies involved, and generate some ¨turbulence¨. As it can be seen in Figure (10), we have a discontinuity front close to x = 0.6. In this example, taken from [7], we used 64 points in each direction and a timestep of size ∆t = 10−4. Without the use of the exponential filter, the code became unstable. 22 Bruno Costa 6, 4(2004) Figure 10: 5.3 2D Flow through an Obstacle -The Compressible Navier Stokes Equation We now have all the necessary ingredients (mapping, filtering, boundary conditions) to simulate a more realistic example. We chose the classical case of the flow past a circular cylinder. We will use the two–dimensional compressible viscous Navier– Stokes equations in strong conservation form (see [2] and [32]). These equations describe the laws of conservation of mass, momentum and energy when the system is isolated from the influence of external forces. They are much more involved equations than the ones we have been dealing, nevertheless, their derivation can be found in many basic books on Fluid Dynamics, for instance, Anderson [2] presents a clear and thorough construction of these equations in two and three dimensions. The non-dimensional form of these equations in the Cartesian coordinates (x , y) is ∂Q ∂t + ∂F ∂x + ∂G ∂y = 1 Re ( ∂Fν ∂x + ∂Gν ∂y ) (36) where Q =   ρ ρu ρv E   , 6, 4(2004) Spectral Methods for Partial Differential Equations 23 and F =   ρu ρu2 + P ρuv (E + P ) u   , Fν =   0 τ xx τ xy γκ Pr−1 ∂T ∂x + uτ xx + vτ xy   , (37) G =   ρv ρuv ρv2 + P (E + P ) v   , Gν =   0 τ xy τ yy γκ Pr−1 ∂T ∂y + uτ xy + vτ yy   . The variables are the density ρ, the velocity field (u, v) and the total energy E. The parameter γ = cρ cv is the ratio of specific heats at constant density and volume, Pr is the Prandtl number and κ is the coefficient of thermal conductivity. The stress tensor elements are given by τ xx = 2 3 µ ( 2 ∂u ∂x − ∂v ∂y ) , τ xy = µ ( ∂u ∂y + ∂v ∂x ) , τ yy = 2 3 µ ( 2 ∂v ∂y − ∂u ∂x ) . and the viscosity µ is related to the temperature by the Sutherland viscosity law µ = (1 + S)T 3 2 S + T , (38) where S is a thermodynamic non-dimensionalized constant. The Navier Stokes equations are coupled with the non-dimensional equations of state and internal energy, respectively: P = (γ − 1) ρ T, E = ρ [ T + 1 2 (u2 + v2) ] . The Navier–Stokes equations (36) and their necessary supplementary relations such as the Sutherland viscosity law are characterized by four parameters [32]. These parameters are the free–stream Mach number M∞, the free–stream Reynolds number Re, the diameter of the cylinder D, and the dimensional temperature T0. The free– stream Reynolds number Re∞ = ρ∞U∞D/µ∞ of the flow is based on the diameter of the cylinder D, the free–stream velocity U∞, the free–stream density ρ∞ and the dynamic viscosity µ∞, which are the values of these quantities away from the cylinder. The physical domain will be represented by the ring-shaped region shown in Figure (7) below, which is the image of the computational domain (ξ , η) ∈ [0, 2π] × [−1, 1] through the polar mapping. The Navier–Stokes equations (36) written in the coordinates (ξ , η) assume the form ∂Q̄ ∂t + ∂F̄ ∂ξ + ∂Ḡ ∂η = 1 Re ( ∂F̄ν ∂ξ + ∂Ḡν ∂η ) , (39) 24 Bruno Costa 6, 4(2004) where Q̄ = 1 J   ρ ρu ρv E   , F̄ = F J ∂ξ ∂x + G J ∂ξ ∂y , F̄ν = Fν J ∂ξ ∂x + Gν J ∂ξ ∂y , (40) Ḡ = F J ∂η ∂x + G J ∂η ∂y , Ḡν = Fν J ∂η ∂x + Gν J ∂η ∂y ; and J is the Jacobian of the transformation (ξ, η) → (x, y), given by J = ξxηy −ξyηx. All the derivative terms in (37) have to be computed with respect to the coordi- nates (ξ, η). Therefore, the elements of the stress tensor are given by τ xx = 2 3 µ [ 2 ( ∂u ∂ξ ∂ξ ∂x + ∂u ∂η ∂η ∂x ) − ( ∂v ∂ξ ∂ξ ∂y + ∂v ∂η ∂η ∂y )] , τ xy = µ [ ( ∂u ∂ξ ∂ξ ∂y + ∂u ∂η ∂η ∂y ) + ( ∂v ∂ξ ∂ξ ∂x + ∂v ∂η ∂η ∂x )] , τ yy = 2 3 µ [ 2 ( ∂v ∂ξ ∂ξ ∂y + ∂v ∂η ∂η ∂y ) − ( ∂u ∂ξ ∂ξ ∂x + ∂u ∂η ∂η ∂x )] and the gradient of the temperature T : ∂T ∂x = ∂T ∂ξ ∂ξ ∂x + ∂T ∂η ∂η ∂x , ∂T ∂y = ∂T ∂ξ ∂ξ ∂y + ∂T ∂η ∂η ∂y . In this problem, we take γ = 1.4. and Pr = 0.72 Figure (11) shows a contour plot of the density function ρ after the trail of vortices started to develop, as the streamlines of the velocity field show. Go to www.labma. ufrj.br/ ˜bcosta to see the movie. We used 64 points in both directions and applied the Kosloff-Tal-Ezer mapping in the radial direction. An exponential filter of order 16 was also applied. The boundary conditions were imposed through the characteristic decomposition as above, in order to simulate a free boundary. For more information on the characteristic variables decomposition and further results on the flow through an obstacle problem, the interested reader can check [22]. 6, 4(2004) Spectral Methods for Partial Differential Equations 25 Figure 11: Contour density plot for the flow through an obstacle problem. Equations (36) and (39) can be written in the general form ∂Q ∂t + ∇x,y · (F, G) = 1 Re ∇x,y · (Fv, Gv) (41) and one only has to note that ∇ξ,η · (F̄ , Ḡ) is the conservative form of the Divergence ∇x,y · (F, G) in the general Cartesian coordinates (x, y)(see [2]). It is now worth pointing that the software PseudoPack2000 has all the tools needed for the trans- formation from a classical set of computational variables like Fourier, Chebyshev or Legendre to a general set of curvilinear coordinates like polar, cylindrical, spherical or user defined, including the computation of the metrics and the Jacobian. Below, we show the one-page code in Fortran 90 necessary to implement the main part of the Navier-Stokes program to generate Figure (11) when using the subroutines of PseudoPack2000. 26 Bruno Costa 6, 4(2004) !********************************************************************** ! Compute the velocities u and v and the Temperature !********************************************************************** u = Q(:,:,2)/Q(:,:,1) v = Q(:,:,3)/Q(:,:,1) T = Q(:,:,4)/Q(:,:,1) T = T - HALF*(u**2+v**2) !********************************************************************** ! Compute the viscosity array mu !********************************************************************** mu = (ONE+S)*T*SQRT(T)/(S+T) !********************************************************************** ! Compute the derivatives of u and v with respect to x and y !********************************************************************** Call PSGrad(Dt, Dr, u, ux, uy, radx, rady, thx, thy, Jacob) Call PSGrad(Dt, Dr, v, vx, vy, radx, rady, thx, thy, Jacob) !********************************************************************** ! Compute the Stress Terms !********************************************************************** Tauxx = rre*mu*(c43*ux-c23*vy) Tauxy = rre*mu*( uy+ vx) Tauyy = rre*mu*(c43*vy-c23*ux) !********************************************************************** ! Compute Tx and Ty !********************************************************************** Call PSGrad(Dt, Dr, T, T x, T y, radx, rady, thx, thy,Jacob) Tx = cq*rre*mu*Tx Ty = cq*rre*mu*Ty !********************************************************************** ! Compute the Pressure term !********************************************************************** P = Q(:,:,1)*T*gm1 !********************************************************************** ! Compute the Divergence of (F-Fv, G-Gv) in the conservative form !********************************************************************** F(:,:,1) = Q(:,:,2) F(:,:,2) = Q(:,:,2) *u - Tauxx + P F(:,:,3) = Q(:,:,3) *u - Tauxy F(:,:,4) = (Q(:,:,4)+P)*u - Tauxx*u - Tauxy*v - Tx G(:,:,1) = Q(:,:,3) G(:,:,2) = Q(:,:,2) *v - Tauxy G(:,:,3) = Q(:,:,3) *v - Tauyy + P G(:,:,4) = (Q(:,:,4)+P)*v - Tauxy*u - Tauyy*v - Ty do i = 1,4 Call PSDiv(Dt, Dr, F (:, :, i), G(:, :, i), DF lux(:, :, i), radx, rady, thx, thy,Jacob, 1) enddo 6, 4(2004) Spectral Methods for Partial Differential Equations 27 5.4 Multidomain As we saw in Section 4, the bases used in spectral methods are sets of eigenfunctions of Sturm-Liouville problems, having limitations with the types of domain that they can be used at. A common way to represent a solution in complex geometries is to use a multidomain setting, where each domain element is of a simple type, or can be mapped to a simple one. For instance, the figure below shows a multidomain representation of the circular obstacle problem with an inner tube through the obstacle. Figure 12: Multidomain grid for the cylinder with an inner tube. This setting came from the general problem of Vortex Induced Vibrations (VIV) common to submarine raisers in oil extraction. It was noticed experimentally that positioning the tube in order to allow passage of the fluid through the obstacle would eliminate the ressonance of the obstacle induced by the vortex generation. In this simple representation, we are only extending the problem of last section to a mul- tidomain configuration, therefore, the equations are being solved separatly in each subdomain. For pasting information through the domain interfaces, we are simply averaging the solution, ensuring its continuity. While being a very simple setting, it shows the potentiallity of the transfinite mapping associated to a multidomain spec- tral scheme. It is also worth pointing that heavy filtering is necessary in order to sustain stability on the corners of the tube. In this experiment we applied a global filtering of the exponential type, however, this can be improved by also applying a heavier filter only to the solution at the corners. There is a whole new area of research for the multidomain setting, the spectral element method, where the idea is to use much smaller elements than the ones above ( see [25]). 28 Bruno Costa 6, 4(2004) 5.5 The Fortran 90 Library PseudoPack2000 While several software tools for the solution of partial differential equations (PDEs) exist in the commercial (e.g. DiffPack) as well as the public domain (e.g. PETSc), they are almost exclusively based on the use of low-order finite difference, finite ele- ment or finite volume methods. Geometric flexibility is one of their main advantages. For most PDE solvers employing pseudospectral (collocation) methods, one major component of the computational kernel is the differentiation operator. The differen- tiation must be done accurately and efficiently on a given computational platform for a successful numerical simulation. It is not an easy task given the number of choice of algorithms for each new and existing computational platform. Issues involving the complexity of the coding, efficient implementation, geometric restriction and lack of high quality software library tended to discourage the general use of pseudospectral methods in scientific research and practical applications. In particular, the lack of standard high quality library for pseudospectral methods forced individual researchers to build codes that were not optimal in efficiency and accuracy. Furthermore, while pseudospectral methods are at a fairly mature level, many critical issues regarding efficiency and accuracy have only recently been addressed and resolved. The knowl- edge of these solutions is not widely known and appears to restrict a more general usage of such schemes. When developing the software PseudoPack2000, the authors aimed to develop a numerical software to make available to the user, in a high performance computing environment, a library of subroutines providing an accurate, versatile, optimal and efficient implementation of the basic components of global pseudospectral methods on which to address a variety of applications of interest to scientists. This library provides subroutines for computing the derivative of the Fourier collo- cation methods for periodical domain and Chebyshev and Legendre collocation meth- ods for the non-periodical domain. State-of-the-art numerical techniques such as Even-Odd Decomposition [33], and specialized fast algorithms are employed to in- crease the efficiency of the library. Kosloff-Tal-Ezer Mapping is used for reducing roundoff error in the Chebyshev and Legendre collocation methods [17]. Moreover, highly accurate Lagrange polynomial interpolation is used when forming the differ- entiation matrices [12], decreasing solution contamination by roundoff error when dealing with a large number of grid points. Other routines for filtering and grid map- ping are also included. Since the user is shielded from any coding errors in the main computational kernels, reliability of the solution is enhanced. PseudoPack will speed up code development, increase productivity and enhance re-usability. The library contains several simple user callable subroutines that return the deriva- tives and/or filtering (smoothing) of, possibly multi-dimensional, data sets. In term of flexibility and user interaction, any aspect of the library can be modified by a simple change of a small set of input parameters. All source codes are written in Fortran 90. Using the macro and conditional capability of the C Preprocessor, this software package can be compiled into several versions with several different computational platforms. Several popular computational platforms (IBM RS6000, SGI Cray, SGI, SUN) are supported to take advantages of any existing optimized native library such 6, 4(2004) Spectral Methods for Partial Differential Equations 29 as General Matrix-Matrix Multiply (GEMM) from Basic Linear Algebra Level 3 Sub- routine (BLAS 3), Fast Fourier Transform (FFT) and Fast Cosine/Sine Transform (CFT/SFT).The Fortran 90 Library PseudoPack2000 is available by request at the webpage of the authors and it comes with a manual. 6 Conclusions The main goal of this article was to provide a more than rough introduction to the subject of spectral methods and their applications in the study of PDEs and the field of the Dynamics of Fluids. We followed the approach of starting from the very begining and only builting the necessary elements to arrive at a meaningful application, like the flow through an obstacle. Thus, many points of relevant interest were ignored and, sometimes, only quickly mentioned through the text. Among them we can cite time integration processes and their consequences towards numerical stability. Nevertheless, the reader with some experience on FD or FE noticed that the mere substitution of the discrete differentiation operator by the spectral one will make most of the temporal integration formulae work with higher precision, with a much smaller time step, however. The interested reader might now check the cited literature and a very effective starting point is the book of Trefethen [35], containing the essentials of spectral col- location methods and a large range of applications. The book is rich of practical examples using a very powerful tool for scientific computing, the Matlab language. Those interested on a more complete overview will be very pleased. For more exten- sive insights, the work of Canuto et al [8] is a good encyclopedia in spectral methods directed towards Fluid Dynamics applications, as well as the more recent monography by Fornberg [18], where a close relationship between collocation spectral methods and high order Finite Diference schemes is developed. The book by Boyd [6] also con- tains plenty information on spectral methods and several areas of research, as well as a substantious discussion on spherical harmonics and other bases for problems in the sphere, with applications to Meteorology. The reader interested in Multidomain spectral methods can found state-of-the-art information in the book by Karniadakis. Finally, the work of Gottlieb and Orszag [20] will satisfy those looking for the math- ematical foundations of spectral methods in a modern setting. References [1] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972. [2] Anderson, J. D.; Degrez, G.; Dick, E.; Grundmann, R. Computational fluid dy- namics. An introduction. Edited and with a preface by John F. Wendt. A von Karman Institute Book. Springer-Verlag, Berlin, 1992. 30 Bruno Costa 6, 4(2004) [3] J. Augenbaum, An adaptive pseudospectral method for discontinuous problems, ICASE report No. 88-54. Appl. Numer. Math. (1988) [4] M. R. Abril-Raimundo and B. Garćia-Archilla, Approximation properties of a mapped Chebyshev method, Appl. Numer. Math., Vol. 32, (2000) pp. 119-136. [5] C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers, McGraw Hill Book, New York (1978) [6] J.P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd ed., Dover, New York, 2000. [7] Calgaro, C.; Laminie, J.; Temam, R. Dynamical multilevel schemes for the solution of evolution equations by hierarchical finite element discretization. Appl. Numer. Math. 23 (1997), no. 4, 403–442. [8] C. Canuto, M. Y. Hussaini, A. Quarteroni and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics, Springer, Berlin, 1988. [9] J. M. Carcione, A 2D Chebyshev Differential Operator for the Elastic Wave Equa- tion, Comput. Meth. Appl. Mech. Engr., Vol. 130, (1996), pp. 33-45. [10] J. W. Cooley, P. A. Lewis and P. D. Welch,The Fast Fourier Transform Al- gorithm: Programming considerations in the Calculations of Sine, Cosine and Laplace Transforms,J. Sound Vib, Vol 12, No. 3, pp 315-337 (1970) [11] B. Costa, L. Dettori, D. Gottlieb and R. Temam, Time Marching Techniques for the Nonlinear Galerkin Method, SIAM Journal of Sci. Comput., Vol. 23 (2001), No. 1, pp. 46-65. [12] B. Costa and W. S. Don, On the Computation of Higher Order Pseudospectral Derivatives, Appl. Numer. Math., Vol. 33, No. 1, (2000), pp. 151-159. [13] B. Costa, W. S. Don and A.Simas, Spectral Convergence of Mapped Chebyshev Methods, submitted. [14] B. Costa and W.S. Don, General Mappings for the Chebyshev Collocation Method, in preparation. [15] T. J. Dekker, A Floating-Point Technique for Extending the Available Precision, Numer. Math. (18), pp. 224-244 (1971) [16] W. S. Don and A. Solomonoff, Accuracy and speed in computing the Chebyshev collocation derivative, SIAM J. Sci. Comp, Vol.16, No. 6,pp. 1253-1268 [17] W. S. Don and A. Solomonoff, Accuracy Enhancement for Higher Derivatives Using Chebyshev Collocation and a Mapping Technique, SIAM J. Sci. Comput., Vol. 18, (1997) pp. 1040-1055. 6, 4(2004) Spectral Methods for Partial Differential Equations 31 [18] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, New York, 1996. [19] W. J. Gordon and C. A. Hall, Transfinite Element Methods: Blending-Function Interpolation over Arbitrary Curved Elements Domains, Numer. Math. 21, pp. 109-129, 1973. [20] D. Gottlieb, S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, SIAM, Philadelphia, PA, 1977. [21] D. Gottlieb and W. S. Don, The Chebyshev-Legendre Method: Implementing Legendre Methods on Chebyshev Points, SIAM J. Numer. Anal., Vol. 40, No. 5, pp. 1666-1682. [22] D. Gottlieb and W. S. Don, Spectral Simulation of Unsteady Compressible Flow Past a Circular Cylinder [23] D. Gottlieb and J. S. Hesthaven, Spectral Methods for Hyperbolic Problems, Jour- nal of Computational and Applied Mathematics, 128 (2001), 83-131. [24] J. S. Hesthaven, P. G. Dinesen and J. P. Lynov, Spectral Collocation Time- Domain Modeling of Diffractive Optical Elements, J. Comput. Phys., Vol. 155, (1999), pp. 287-306. [25] Karniadakis, George Em; Sherwin, Spencer J. Spectral/$hp$ element methods for CFD. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, 1999. [26] D. Kosloff and H. Tal-Ezer, Modified Chebyshev Pseudospectral Method With O(N−1) Time Step Restriction, J. Comput. Phys., Vol. 104, (1993), pp. 457-469. [27] A. I. Markushevich, Theory of Functions of a Complex Variable, Chelsea, New York (1977). [28] J. L. Mead and R. A. Renaut, Accuracy, Resolution and Stability Properties of a Modified Chebyshev Method, SIAM J. Sci. Comput., Vol. 24, No. 1, (2002), pp. 143-160. [29] S. A. Orszag, Comparison of Pseudospectral and Spectral Approximation, Studies in Applied Mathematics, Vol. LI, No. 3, 1972, pp. 253-259. [30] S. C. Reddy, J. A. C. Weideman and G. F. Norris, On a Modified Chebyshev Pseudospectral Method, (unpublished, Oregon State University, Corvalis, OR). [31] R. A. Renaut and Y. Su, Evaluation of Chebyshev Pseudospectral Methods for Third Order Differential Equations, Numer. Algo., Vol. 16, No. 3, (1997), pp. 255-281. [32] P. J. Roache, Computational Fluid Dynamics, Hermosa Publishers, pp 209-223, 1976. 32 Bruno Costa 6, 4(2004) [33] Solomonoff, Alex A fast algorithm for spectral differentiation. J. Comput. Phys. 98 (1992), no. 1, 174–177. [34] P. Swarztrauber, Symmetric FFT’s, Math. Comp., Vol 47, No. 175,pp. 323-346 (1986) [35] L. N. Trefethen, Spectral Methods in Matlab, SIAM, Philadelphia, 2000. [36] L. N. Trefethen and M. R. Trummer, An Instability Phenomenon in Spectral Methods, SIAM J. Numer. Anal., Vol. 24, No. 5, (1987) [37] J. A. C. Weideman and L. N. Trefethen, The Eigenvalues of Second-Order Spec- tral Differentiation Matrices, SIAM J. Numer. Anal., Vol. 23, No. 6, (1988)