10 CHAPTER 1.INTRODUCTION Although we shall discuss many types of basis functions,the best choice for 95%of all applications is an ordinary Fourier series,or a Fourier series in disguise.By "disguise"we mean a change of variable which turns the sines and cosines of a Fourier series into differ- ent functions.The most important disguise is the one worn by the Chebyshev polynomials, which are defined by Tn(cosf)≡cos(nd) (1.16) Although the Tn(z)are polynomials in z,and are therefore usually considered a sepa- rate and distinct species of basis functions,a Chebyshev series is really just a Fourier cosine expansion with a change of variable.This brings us to the first of our proverbial sayings: MORAL PRINCIPLE 1: (i)When in doubt,use Chebyshev polynomials unless the solution is spatially periodic,in which case an ordinary Fourier series is better. (ii)Unless you're sure another set of basis functions is better,use Chebyshev polynomials. (iii)Unless you're really,really sure that another set of basis functions is better,use Chebyshev polynomials. There are exceptions:on the surface of a sphere,it is more efficient to use spherical harmonics than Chebyshev polynomials.Similarly,if the domain is infinite or semi-infinite, it is better to use basis sets tailored to those domains than Chebyshev polynomials,which in theory and practice are associated with a finite interval. The general rule is:Geometry chooses the basis set.The engineer never has to make a choice.Table A-1 in Appendix A and Figure 1.6 summarize the main cases.When multiple basis sets are listed for a single geometry or type of domain,there is little to choose between them. It must be noted,however,that the non-Chebyshev cases in the table only strengthen the case for our first Moral Principle.Though not quite as good as spherical harmonics, Chebyshev polynomials in latitude and longtitude work just fine on the sphere (Boyd, 1978b).The rational Chebyshev basis sets are actually just the images of the usual Cheby- shev polynomials under a change of coordinate that stretches the interval [-1,1]into an infinite or semi-infinite domain.Chebyshev polynomials are,as it were,almost idiot-proof. Consequently,our analysis will concentrate almost exclusively upon Fourier series and Chebyshev polynomials.Because these two basis sets are the same except for a change of variable,the theorems for one are usually trivial generalizations of those for the other. The formal convergence theory for Legendre polynomials is essentially the same as for Chebyshev polynomials except for a couple of minor items noted in Chapter 2.Thus, understanding Fourier series is the key to all spectral methods. 1.7 Boundary conditions Normally,boundary and initial conditions are not a major complication for spectral meth- ods.For example,when the boundary conditions require the solution to be spatially pe- riodic.the sines and cosines of a Fourier series (which are the natural basis functions for all periodic problems)automatically and individually satisfy the boundary conditions.Con- sequently,our only remaining task is to choose the coefficients of the Fourier series to minimize the residual function
10 CHAPTER 1. INTRODUCTION Although we shall discuss many types of basis functions, the best choice for 95% of all applications is an ordinary Fourier series, or a Fourier series in disguise. By “disguise” we mean a change of variable which turns the sines and cosines of a Fourier series into different functions. The most important disguise is the one worn by the Chebyshev polynomials, which are defined by Tn(cosθ) ≡ cos(nθ) (1.16) Although the Tn(x) are polynomials in x, and are therefore usually considered a separate and distinct species of basis functions, a Chebyshev series is really just a Fourier cosine expansion with a change of variable. This brings us to the first of our proverbial sayings: MORAL PRINCIPLE 1: (i) When in doubt, use Chebyshev polynomials unless the solution is spatially periodic, in which case an ordinary Fourier series is better. (ii) Unless you’re sure another set of basis functions is better, use Chebyshev polynomials. (iii) Unless you’re really, really sure that another set of basis functions is better, use Chebyshev polynomials. There are exceptions: on the surface of a sphere, it is more efficient to use spherical harmonics than Chebyshev polynomials. Similarly, if the domain is infinite or semi-infinite, it is better to use basis sets tailored to those domains than Chebyshev polynomials, which in theory and practice are associated with a finite interval. The general rule is: Geometry chooses the basis set. The engineer never has to make a choice. Table A-1 in Appendix A and Figure 1.6 summarize the main cases. When multiple basis sets are listed for a single geometry or type of domain, there is little to choose between them. It must be noted, however, that the non-Chebyshev cases in the table only strengthen the case for our first Moral Principle. Though not quite as good as spherical harmonics, Chebyshev polynomials in latitude and longtitude work just fine on the sphere (Boyd, 1978b). The rational Chebyshev basis sets are actually just the images of the usual Chebyshev polynomials under a change of coordinate that stretches the interval [-1, 1] into an infinite or semi-infinite domain. Chebyshev polynomials are, as it were, almost idiot-proof. Consequently, our analysis will concentrate almost exclusively upon Fourier series and Chebyshev polynomials. Because these two basis sets are the same except for a change of variable, the theorems for one are usually trivial generalizations of those for the other. The formal convergence theory for Legendre polynomials is essentially the same as for Chebyshev polynomials except for a couple of minor items noted in Chapter 2. Thus, understanding Fourier series is the key to all spectral methods. 1.7 Boundary conditions Normally, boundary and initial conditions are not a major complication for spectral methods. For example, when the boundary conditions require the solution to be spatially periodic, the sines and cosines of a Fourier series (which are the natural basis functions for all periodic problems) automatically and individually satisfy the boundary conditions. Consequently, our only remaining task is to choose the coefficients of the Fourier series to minimize the residual function
1.7.BOUNDARY CONDITIONS 11 Periodic Non-Periodic Fourier Chebyshev or Legendre θ∈[0,2] ×∈[-1,1] Semi-Infinite Infinite Rational Cheby.TL Rational Cheby.TB or Laguerre or Hermite or Sinc X∈[0，∞] X∈[-∞，∞] Figure 1.6:Choice of basis functions.Upper left:on a periodic interval,use sines and cosines.This case is symbolized by a ring because the dependence on an angular coordi- nate,such as longitude,is always periodic.Upper right:a finite interval,which can always be rescaled and translated to zE[-1,1].Chebyshev or Legendre polynomials are optimal. Lower left:semi-infinite interval z E [0,oo],symbolized by a one-sided arrow.Rational Chebyshev functions TLn(z)are the generic choice,but Laguerre functions are sometimes more convenient for particular problems.Lower right:xE [-oo,oo](double-ended ar- row).Rational Chebyshev functions TBn(z)are the most general,but sinc and Hermite functions are widely used,and have similar convergence properties. For non-periodic problems,Chebyshev polynomials are the natural choice as explained in the next chapter.They do not satisfy the appropriate boundary conditions,but it is easy to add explicit constraints such as anon(1)=a (1.17) n=0 to the algebraic equations obtained from minimizing R(z;ao,a1,...,an)so that u(1)=a is satisfied by the approximate solution
1.7. BOUNDARY CONDITIONS 11 Periodic Non-Periodic Semi-Infinite Infinite x ∈ [0, ∞] x ∈ [−1, 1] x ∈ [− ∞, ∞] θ ∈ [0, 2 π] Fourier Rational Cheby. TL or Laguerre Chebyshev or Legendre Rational Cheby. TB or Hermite or Sinc Figure 1.6: Choice of basis functions. Upper left: on a periodic interval, use sines and cosines. This case is symbolized by a ring because the dependence on an angular coordinate, such as longitude, is always periodic. Upper right: a finite interval, which can always be rescaled and translated to x ∈ [−1, 1]. Chebyshev or Legendre polynomials are optimal. Lower left: semi-infinite interval x ∈ [0,∞], symbolized by a one-sided arrow. Rational Chebyshev functions T Ln(x) are the generic choice, but Laguerre functions are sometimes more convenient for particular problems. Lower right: x ∈ [−∞,∞] (double-ended arrow). Rational Chebyshev functions T Bn(x) are the most general, but sinc and Hermite functions are widely used, and have similar convergence properties. For non-periodic problems, Chebyshev polynomials are the natural choice as explained in the next chapter. They do not satisfy the appropriate boundary conditions, but it is easy to add explicit constraints such as X N n=0 an φn(1) = α (1.17) to the algebraic equations obtained from minimizing R(x; a0, a1,... ,aN ) so that u(1) = α is satisfied by the approximate solution
12 CHAPTER 1.INTRODUCTION Alternatively,one may avoid explicit constraints like(1.17)by writing the solution as u(x)≡v(x)+w(x) (1.18) where w(z)is a known function chosen to satisfy the inhomogeneous boundary conditions. The new unknown,v(z),satisfies homogeneous boundary conditions.For(1.17),for ex- ample,w(1)=a and v(1)=0.The advantage of homogenizing the boundary conditions is that we may combine functions of the original basis,such as the Chebyshev polynomials, into new basis functions that individually satisfy the homogeneous boundary conditions. This is surprisingly easy to do;for example,to satisfy v(-1)=v(1)=0,we expand v(x)in terms of the basis functions p2n(x）≡T2n(x)-1, n=1,2, p2m+1(x)三T2n+1(x)-x, n=1,2,… (1.19) where the Tn(z)are the usual Chebyshev polynomials whose properties(including bound- ary values)are listed in Appendix A.This basis is complete for functions which vanish at the ends of the interval.The reward for the switch of basis set is that it is unnecessary. when using basis recombination,to waste rows of the discretization matrix on the bound- ary conditions:All algebraic equations come from minimizing the residual of the differen- tial equation. 1.8 The Two Kingdoms:Non-Interpolating and Pseudospec- tral Families of Methods Spectral methods fall into two broad categories.In the same way that all of life was once divided into the "plant"and "animal"kingdoms2,most spectral methods may be classed as either "interpolating"or "non-interpolating".Of course,the biological classification may be ambiguous-is a virus a plant or animal?How about a sulfur-eating bacteria?The mathematical classification may be ambiguous,too,because some algorithms mix ideas from both the interpolating and non-interpolating kingdoms.Nonetheless,the notion of two exclusive "kingdoms"is a useful taxonomical starting point for both biology and nu- merical analysis. The "interpolating"or "pseudospectral"methods associate a grid of points with each basis set.The coefficients of a known function f()are found by requiring that the trun- cated series agree with f(x)at each point of the grid.Similarly,the coefficients an of a pseudospectral approximation to the solution of a differential equation are found by re- quiring that the residual function interpolate f =0: R(x;ao,a1,.,aw）=0,i=0,1,,N (1.20) In words,the pseudospectral method demands that the differential equation be exactly satisfied at a set of points known as the "collocation"or "interpolation"points.Presum- ably,as R(;an)is forced to vanish at an increasingly large number of discrete points,it will be smaller and smaller in the gaps between the collocation points so that Rz every- where in the domain,and therefore un(z)will converge to u(z)as N increases.Methods in this kingdom of algorithms are also called "orthogonal collocation"or "method of selected points”. The "non-interpolating"kingdom of algorithms includes Galerkin's method and the Lanczos tau-method.There is no grid of interpolation points.Instead,the coefficients of 2Modern classification schemes use three to five kingdoms,but this doesn't change the argument
12 CHAPTER 1. INTRODUCTION Alternatively, one may avoid explicit constraints like (1.17) by writing the solution as u(x) ≡ v(x) + w(x) (1.18) where w(x) is a known function chosen to satisfy the inhomogeneous boundary conditions. The new unknown, v(x), satisfies homogeneous boundary conditions. For (1.17), for example, w(1) = α and v(1) = 0. The advantage of homogenizing the boundary conditions is that we may combine functions of the original basis, such as the Chebyshev polynomials, into new basis functions that individually satisfy the homogeneous boundary conditions. This is surprisingly easy to do; for example, to satisfy v(−1) = v(1) = 0, we expand v(x) in terms of the basis functions φ2n(x) ≡ T2n(x) − 1, n = 1, 2,... φ2n+1(x) ≡ T2n+1(x) − x, n = 1, 2,... (1.19) where the Tn(x) are the usual Chebyshev polynomials whose properties (including boundary values) are listed in Appendix A. This basis is complete for functions which vanish at the ends of the interval. The reward for the switch of basis set is that it is unnecessary, when using basis recombination, to waste rows of the discretization matrix on the boundary conditions: All algebraic equations come from minimizing the residual of the differential equation. 1.8 The Two Kingdoms: Non-Interpolating and Pseudospectral Families of Methods Spectral methods fall into two broad categories. In the same way that all of life was once divided into the “plant” and “animal” kingdoms2, most spectral methods may be classed as either “interpolating” or “non–interpolating”. Of course, the biological classification may be ambiguous – is a virus a plant or animal? How about a sulfur-eating bacteria? The mathematical classification may be ambiguous, too, because some algorithms mix ideas from both the interpolating and non-interpolating kingdoms. Nonetheless, the notion of two exclusive “kingdoms” is a useful taxonomical starting point for both biology and numerical analysis. The “interpolating” or “pseudospectral” methods associate a grid of points with each basis set. The coefficients of a known function f(x) are found by requiring that the truncated series agree with f(x) at each point of the grid. Similarly, the coefficients an of a pseudospectral approximation to the solution of a differential equation are found by requiring that the residual function interpolate f ≡ 0: R(xi; a0, a1,... ,aN )=0, i = 0, 1, ..., N (1.20) In words, the pseudospectral method demands that the differential equation be exactly satisfied at a set of points known as the “collocation” or “interpolation” points. Presumably, as R(x; an) is forced to vanish at an increasingly large number of discrete points, it will be smaller and smaller in the gaps between the collocation points so that R ≈ x everywhere in the domain, and therefore uN (x) will converge to u(x) as N increases. Methods in this kingdom of algorithms are also called “orthogonal collocation” or “method of selected points”. The “non–interpolating” kingdom of algorithms includes Galerkin’s method and the Lanczos tau-method. There is no grid of interpolation points. Instead, the coefficients of 2Modern classification schemes use three to five kingdoms, but this doesn’t change the argument
1.9.NONLINEARITY 13 a known function f(z)are computed by multiplying f(z)by a given basis function and integrating.It is tempting to describe the difference between the two algorithmic kingdoms as "integration"versus "interpolation",but unfortunately this is a little simplistic.Many older books,such as Fox and Parker(1968),show how one can use the properties of the basis functions-recurrence relations,trigonometric identities,and such -to calculate co- efficients without explicitly performing any integrations.Even though the end product is identical with that obtained by integration,it is a little confusing to label a calculation as an "integration-type"spectral method when there is not an integral sign in sight!Therefore, we shall use the blander label of "non-interpolating". Historically,the "non-interpolating"methods were developed first.For this reason, the label "spectral"is sometimes used in a narrow sense as a collective tag for the "non- interpolating"methods.In these notes,we shall use "spectral"only as catch-all for global expansion methods in general,but the reader should be aware of its other,narrower usage. (Actually,there are several other uses because "spectral"has other meanings in time series analysis and functional analysis-ugh!). Many spectral models of time-dependent hydrodynamics split the calculation into sev- eral subproblems and apply different techniques to different subproblems.To continue with our biological metaphor,the computer code becomes a little ecology of interacting “interpolating'”and“non-interpolating”parts.Each algorithm(or algorithmic kingdom) has ecological niches where it is superior to all competition,so we must master both non- interpolating and pseudospectral methods. At first glance,there is no obvious relation between the pseudospectral method and the alternatives that use weighted integrals of R(z;an)to choose the [an).Worse still,we now have the further burden of choosing the interpolation points,{ri).Fortunately,there is a natural choice of interpolation points for each of the common basis sets.These points are the Gaussian quadrature points for the integrals of Galerkin's method.The pseudospectral method is therefore equivalent to the spectral if we evaluate the integrals of the latter by numerical quadrature with(N+1)points.This is the reason why the interpolation-based methods are now commonly called "pseudospectral". Better yet,we shall show in later chapters that the accuracy of pseudospectral methods is only a little bit poorer than that of the non-interpolating kingdom-too little to outweigh the much greater simplicity and computational efficiency of the pseudospectral algorithms. Consequently,we shall emphasize pseudospectral methods in this book.Nonetheless,the justification for the pseudospectral kingdom is derived from that for the non-interpolating methods,and the latter are still superior to interpolation for specialized but important applications.We cannot understand the high efficiency of either kingdom of spectral algo- rithms without first reviewing the theory of Fourier series(Chapter 2). 1.9 Nonlinearity Nonlinearity is not a major complication for spectral methods per se.For expository sim- plicity,we shall usually concentrate on linear algorithms,especially in explaining the basic ideas.The extension to nonlinear problems usually only requires minor modifications To illustrate this,we shall do a very simple nonlinear boundary value problem here. Such equations normally are solved by Newton's iteration.If we set up the iteration by first linearizing the differential equation about the current iterate(the "Newton-Kantorovich" method of Appendix C),then we solve a linear differential equation at each step. In this example,we shall instead apply the spectral method first.The only difference from a linear problem is that the system of algebraic equations for the coefficients is nonlin- ear.It is usually irrelevant whether the Newton iteration is created before or after applying
1.9. NONLINEARITY 13 a known function f(x) are computed by multiplying f(x) by a given basis function and integrating. It is tempting to describe the difference between the two algorithmic kingdoms as ”integration” versus ”interpolation”, but unfortunately this is a little simplistic. Many older books, such as Fox and Parker (1968), show how one can use the properties of the basis functions – recurrence relations, trigonometric identities, and such – to calculate coefficients without explicitly performing any integrations. Even though the end product is identical with that obtained by integration, it is a little confusing to label a calculation as an ”integration-type” spectral method when there is not an integral sign in sight! Therefore, we shall use the blander label of ”non-interpolating”. Historically, the “non–interpolating” methods were developed first. For this reason, the label “spectral” is sometimes used in a narrow sense as a collective tag for the “non– interpolating” methods. In these notes, we shall use “spectral” only as catch–all for global expansion methods in general, but the reader should be aware of its other, narrower usage. (Actually, there are several other uses because “spectral” has other meanings in time series analysis and functional analysis – ugh!). Many spectral models of time-dependent hydrodynamics split the calculation into several subproblems and apply different techniques to different subproblems. To continue with our biological metaphor, the computer code becomes a little ecology of interacting “interpolating” and “non–interpolating” parts. Each algorithm (or algorithmic kingdom) has ecological niches where it is superior to all competition, so we must master both non– interpolating and pseudospectral methods. At first glance, there is no obvious relation between the pseudospectral method and the alternatives that use weighted integrals of R(x; an) to choose the {an}. Worse still, we now have the further burden of choosing the interpolation points, {xi}. Fortunately, there is a natural choice of interpolation points for each of the common basis sets. These points are the Gaussian quadrature points for the integrals of Galerkin’s method. The pseudospectral method is therefore equivalent to the spectral if we evaluate the integrals of the latter by numerical quadrature with (N + 1) points. This is the reason why the interpolation-based methods are now commonly called “pseudospectral”. Better yet, we shall show in later chapters that the accuracy of pseudospectral methods is only a little bit poorer than that of the non-interpolating kingdom – too little to outweigh the much greater simplicity and computational efficiency of the pseudospectral algorithms. Consequently, we shall emphasize pseudospectral methods in this book. Nonetheless, the justification for the pseudospectral kingdom is derived from that for the non-interpolating methods, and the latter are still superior to interpolation for specialized but important applications. We cannot understand the high efficiency of either kingdom of spectral algorithms without first reviewing the theory of Fourier series (Chapter 2). 1.9 Nonlinearity Nonlinearity is not a major complication for spectral methods per se. For expository simplicity, we shall usually concentrate on linear algorithms, especially in explaining the basic ideas. The extension to nonlinear problems usually only requires minor modifications. To illustrate this, we shall do a very simple nonlinear boundary value problem here. Such equations normally are solved by Newton’s iteration. If we set up the iteration by first linearizing the differential equation about the current iterate (the “Newton-Kantorovich” method of Appendix C), then we solve a linear differential equation at each step. In this example, we shall instead apply the spectral method first. The only difference from a linear problem is that the system of algebraic equations for the coefficients is nonlinear. It is usually irrelevant whether the Newton iteration is created before or after applying
14 CHAPTER 1.INTRODUCTION the spectral method;take your choice! 09 08 07- 06 3o 04 03 02 01 0 02 04 608 Figure 1.7:Comparison of exact and approximate solutions:Nonlinear diffusion equation The nonlinear boundary value problem is uzz +a(uz)2+quuzz =0 (1.21) subject to the boundary conditions that u(0)=0:u(1)=1 (1.22) We will take the approximate solution u2(z)to be a quadratic polynomial.The most gen- eral quadratic polynomial which satisfies the boundary conditions is 2=x+a2(x2-x) (1.23) Since there is only one undetermined coefficient a2,only a single collocation point is needed. The obvious choice,the midpoint of the interval,is best. The residual function is R(x;a2)=aa[6x2-6x+1]+2a2[3ax+1-a+a (1.24) The condition that R(z =1/2;a2)=0 then gives the quadratic equation aa【-1/2+2a2[a/2+1刂+a=0 (1.25) We note an amusing fact:although pseudospectral methods are usually considered only as numerical techniques,we have in fact obtained an analytical solution to this nonlinear problem.To see how accurate it is,let us specialize to a =1 for which the exact solution is u(xa=1)=-1+(1+3x)1/2 (1.26) There are two roots to the quadratic,of course,but one gives an unphysical heat flux to- wards the boundary source at z =1,so it can be rejected.3 3The ambiguity of multiple solutions is a difficulty raised by the nonlinearity of the differential equation,not by the method used to solve it.All algorithms for solving nonlinear boundary value problems have the drawback that the algebraic equations that are the discretization of the differential equation have multiple solutions.Most are unphysical and must be rejected on various grounds including(i)imaginary parts(ii)unrealistic behavior such as the heat flux for this example or(iii)failure to converge as N is varied
14 CHAPTER 1. INTRODUCTION the spectral method; take your choice! Figure 1.7: Comparison of exact and approximate solutions: Nonlinear diffusion equation The nonlinear boundary value problem is uxx + α(ux) 2 + αuuxx = 0 (1.21) subject to the boundary conditions that u(0) = 0; u(1) = 1 (1.22) We will take the approximate solution u2(x) to be a quadratic polynomial. The most general quadratic polynomial which satisfies the boundary conditions is u2 = x + a2(x2 − x) (1.23) Since there is only one undetermined coefficient a2, only a single collocation point is needed. The obvious choice, the midpoint of the interval, is best. The residual function is R(x; a2) = a2 2 α [6x2 − 6x + 1] + 2a2 [3αx + 1 − α] + α (1.24) The condition that R(x = 1/2; a2)=0 then gives the quadratic equation a2 2 α [−1/2] + 2a2 [α/2 + 1] + α = 0 (1.25) We note an amusing fact: although pseudospectral methods are usually considered only as numerical techniques, we have in fact obtained an analytical solution to this nonlinear problem. To see how accurate it is, let us specialize to α = 1 for which the exact solution is u(x; α = 1) = −1 + (1 + 3x) 1/2 (1.26) There are two roots to the quadratic, of course, but one gives an unphysical heat flux towards the boundary source at x = 1, so it can be rejected.3 3The ambiguity of multiple solutions is a difficulty raised by the nonlinearity of the differential equation, not by the method used to solve it. All algorithms for solving nonlinear boundary value problems have the drawback that the algebraic equations that are the discretization of the differential equation have multiple solutions. Most are unphysical and must be rejected on various grounds including (i) imaginary parts (ii) unrealistic behavior such as the heat flux for this example or (iii) failure to converge as N is varied