Errata and Extended-Bibliography These may be found on author's homepage,currently at http://www-personal.engin.umich.edu/~jpboyd Errata and comments may be sent to the author at the fol- lowing: jpboyd@umich.edu Thank you! xvi
Errata and Extended-Bibliography These may be found on author’s homepage, currently at http://www-personal.engin.umich.edu/∼jpboyd Errata and comments may be sent to the author at the following: jpboyd@umich.edu Thank you! xvi
Chapter 1 Introduction "I have no satisfaction in formulas unless I feel their numerical magnitude." -Sir William Thomson,1st Lord Kelvin(1824-1907) "It is the increasingly pronounced tendency of modern analysis to substitute ideas for calculation;nevertheless,there are certain branches of mathematics where calculation con- serves its rights." -P.G.L.Dirichlet(1805-1859) 1.1 Series expansions Our topic is a family of methods for solving differential and integral equations.The basic idea is to assume that the unknown u(z)can be approximated by a sum of N+1 "basis functions'”pn(r: N u(a)≈uw(d)=∑ann(a) (1.1) n=0 When this series is substituted into the equation Lu=f(x) (1.2) where L is the operator of the differential or integral equation,the result is the so-called "residual function"defined by R(r;ao,a1,·,aN）=LuN-f (1.3) Since the residual function R(x;an)is identically equal to zero for the exact solution, the challenge is to choose the series coefficients fan}so that the residual function is mini- mized.The different spectral and pseudospectral methods differ mainly in their minimiza- tion strategies. 1
Chapter 1 Introduction “I have no satisfaction in formulas unless I feel their numerical magnitude.” –Sir William Thomson, 1st Lord Kelvin (1824–1907) “It is the increasingly pronounced tendency of modern analysis to substitute ideas for calculation; nevertheless, there are certain branches of mathematics where calculation conserves its rights.” –P. G. L. Dirichlet (1805–1859) 1.1 Series expansions Our topic is a family of methods for solving differential and integral equations. The basic idea is to assume that the unknown u(x) can be approximated by a sum of N + 1 “basis functions” φn(x): u(x) ≈ uN (x) = X N n=0 an φn(x) (1.1) When this series is substituted into the equation Lu = f(x) (1.2) where L is the operator of the differential or integral equation, the result is the so-called “residual function” defined by R(x; a0, a1,... ,aN ) = LuN − f (1.3) Since the residual function R(x; an) is identically equal to zero for the exact solution, the challenge is to choose the series coefficients {an} so that the residual function is minimized. The different spectral and pseudospectral methods differ mainly in their minimization strategies. 1
2 CHAPTER 1.INTRODUCTION 1.2 First Example These abstract ideas can be made concrete by a simple problem.Although large problems are usually programmed in FORTRAN and C,it is very educational to use an algebraic manipulation language like Maple,Mathematica,Macsyma or Reduce.In what follows, Maple statements are shown in bold face.The machine's answers have been converted into standard mathematical notation. The example is the linear,one-dimensional boundary value problem: uzz-(x6+3x2)u=0 (1.4) u(-1)=u(1)=1 (1.5) The exact solution is (Scraton,1965) u(x)=exp([x4-1]/4) (1.6) Polynomial approximations are recommended for most problems,so we shall choose a spectral solution of this form.In order to satisfy the boundary conditions independently of the unknown spectral coefficients,it is convenient to write the approximation as u2:=1 +(1-x*x)*(a0+al*x+a2*x*x): 2=1+(1-x2)(a0+a1x+a2x2) (1.7) where the decision to keep only three degrees of freedom is arbitrary. The residual for this approximation is Resid:=diff(u,x,x)-(x**6+3*x**2)*u; R(x;a0,a1,a2）=u2,xz-(x6+3x2)u2 (1.8) R=(2a2+2a0）-6a1x-(3+3a0+12a2)z2-3a1x3+3(a0-a2)x4 (1.9) +3a1x5+(-1-a0+3a2)z6-a1x7+(a0-a2)x8+a1x9+10a2x10 As error minimization conditions,we choose to make the residual zero at a set of points equal in number to the undetermined coefficients in u2(z).This is called the "collocation" or "pseudospectral"method.If we arbitrarily choose the points i=(-1/2,0,1/2),this gives the three equations:egl:=subs(x=-1/2,Resid); eq2:=subs(x=0,Resid);eq3:=subs(x=1/2,Resid); 659 1683 117149 egl 25600+ 512a1 1024a2- 64 eq2=-2(a0-a2 eg3 6591683117149 256a0 5121 -10242-64 (1.10) The coefficients are then determined by solving eql=eg2 =eg3=0; solutionarray:=solve([eql,eq2,eq3),{a0,a1,a2));yields 784 Q0=- 3807 a1=0, a2=a0 (1.11) Figure 1.1 shows that this low order approximation is quite accurate. However,the example raises a whole host of questions including:
2 CHAPTER 1. INTRODUCTION 1.2 First Example These abstract ideas can be made concrete by a simple problem. Although large problems are usually programmed in FORTRAN and C, it is very educational to use an algebraic manipulation language like Maple, Mathematica, Macsyma or Reduce. In what follows, Maple statements are shown in bold face. The machine’s answers have been converted into standard mathematical notation. The example is the linear, one-dimensional boundary value problem: uxx − (x6 + 3x2)u = 0 (1.4) u(−1) = u(1) = 1 (1.5) The exact solution is (Scraton, 1965) u(x) = exp([x4 − 1]/4) (1.6) Polynomial approximations are recommended for most problems, so we shall choose a spectral solution of this form. In order to satisfy the boundary conditions independently of the unknown spectral coefficients, it is convenient to write the approximation as u2:=1 + (1-x*x)*(a0 + a1*x + a2*x*x); u2 = 1 + (1 − x2)(a0 + a1x + a2x2) (1.7) where the decision to keep only three degrees of freedom is arbitrary. The residual for this approximation is Resid:= diff(u,x,x) - (x**6 + 3*x**2)*u; R(x; a0, a1, a2) = u2,xx − (x6 + 3x2)u2 (1.8) R = (2a2 + 2a0) − 6a1x − (3 + 3a0 + 12a2)x2 − 3a1x3 + 3(a0 − a2)x4 (1.9) +3a1x5 + (−1 − a0 + 3a2)x6 − a1x7 + (a0 − a2)x8 + a1x9 + 10a2x10 As error minimization conditions, we choose to make the residual zero at a set of points equal in number to the undetermined coefficients in u2(x). This is called the “collocation” or “pseudospectral” method. If we arbitrarily choose the points xi = (−1/2, 0, 1/2), this gives the three equations: eq1:=subs(x=-1/2,Resid); eq2:=subs(x=0,Resid); eq3:=subs(x=1/2,Resid); eq1 = −659 256 a0 + 1683 512 a1 − 1171 1024a2 − 49 64 eq2 = −2(a0 − a2) eq3 = −659 256 a0 − 1683 512 a1 − 1171 1024a2 − 49 64 (1.10) The coefficients are then determined by solving eq1 = eq2 = eq3=0; solutionarray:= solve({eq1,eq2,eq3}, {a0,a1,a2}); yields a0 = − 784 3807, a1 = 0, a2 = a0 (1.11) Figure 1.1 shows that this low order approximation is quite accurate. However, the example raises a whole host of questions including:
1.2.FIRST EXAMPLE 3 1.What is an optimum choice of basis functions? 2.Why choose"collocation"as the residual-minimizing condition? 3.What are the optimum collocation points? 4.Why is a zero?Could we have anticipated this,and used a trial solution with just two degrees of freedom for the same answer? 5.How do we solve the algebraic problem for the coefficients when the Maple "solve" function isn't available? The answer to the first question is that choosing powers ofx as a basis is actually rather dangerous unless N,the number of degrees-of-freedom,is small or the calculations are being done in exact arithmetic,as true for the Maple solution here.In the next section,we describe the good choices.In an algebraic manipulation language,different rules apply as explained in Chapter 20. The second answer is:Collocation is the simplest choice which is guaranteed to work, and if done right,nothing else is superior.To understand why.however,we shall have to understand both the standard theory of Fourier and Chebyshev series and Galerkin meth- ods(Chapters 2 and 3)and the theory of interpolation and cardinal functions (Chapters 4 and 5). The third answer is:once the basis set has been chosen,there are only two optimal sets of interpolation points for each basis(the Gauss-Chebyshev points and the Gauss-Lobatto points);both are given by elementary formulas in Appendices A and F,and which one is used is strictly a matter of convenience. The fourth answer is:Yes,the irrelevance of a could have been anticipated.Indeed,one can show that for this problem,all the odd powers of x have zero coefficients.Symmetries of various kinds are extremely important in practical applications(Chapter 8). Exact-u2 Error 0 -0.002 0.95 -0.004 -0.006 0.9 -0.008 -0.01 0.85 -0.012 o 0.8 °o00000o° -0.014 -0.016 0.7 -0.018 0 0 Figure 1.1:Left panel:Exact solution u exp([x4-1]/4)(solid)is compared with the three-coefficient numerical approximation (circles).Right panel:u-u2
1.2. FIRST EXAMPLE 3 1. What is an optimum choice of basis functions? 2. Why choose “collocation” as the residual-minimizing condition? 3. What are the optimum collocation points? 4. Why is a1 zero? Could we have anticipated this, and used a trial solution with just two degrees of freedom for the same answer? 5. How do we solve the algebraic problem for the coefficients when the Maple “solve” function isn’t available? The answer to the first question is that choosing powers of x as a basis is actually rather dangerous unless N, the number of degrees-of-freedom, is small or the calculations are being done in exact arithmetic, as true for the Maple solution here. In the next section, we describe the good choices. In an algebraic manipulation language, different rules apply as explained in Chapter 20. The second answer is: Collocation is the simplest choice which is guaranteed to work, and if done right, nothing else is superior. To understand why, however, we shall have to understand both the standard theory of Fourier and Chebyshev series and Galerkin methods (Chapters 2 and 3) and the theory of interpolation and cardinal functions (Chapters 4 and 5). The third answer is: once the basis set has been chosen, there are only two optimal sets of interpolation points for each basis (the Gauss-Chebyshev points and the Gauss-Lobatto points); both are given by elementary formulas in Appendices A and F, and which one is used is strictly a matter of convenience. The fourth answer is: Yes, the irrelevance of a1 could have been anticipated. Indeed, one can show that for this problem, all the odd powers of x have zero coefficients. Symmetries of various kinds are extremely important in practical applications (Chapter 8). -1 0 1 0.75 0.8 0.85 0.9 0.95 1 Exact - u2 -1 0 1 -0.018 -0.016 -0.014 -0.012 -0.01 -0.008 -0.006 -0.004 -0.002 0 Error Figure 1.1: Left panel: Exact solution u = exp([x4 − 1]/4) (solid) is compared with the three-coefficient numerical approximation (circles). Right panel: u − u2
CHAPTER 1.INTRODUCTION Table 1.1:Maple program to solve linear boundary-value problem u2:=1+(1-x*x)*(a0+al*x+a2*x*x: Resid:=diff(u,x,x)-(x**6+3*x**2)*u; eq1:=subs(x=-1/2,Resid);eq2:=subs(x=0,Resid);eq3:=subs(x=1/2,Resid); solutionarray:=solve(feq1,eq2,eq3),[a0,a1,a2)); The fifth answer is:the algebraic equations can be written(for a linear differential equa- tion)as a matrix equation,which can then be solved by library software in FORTRAN or C. Many other questions will be asked and answered in later chapters.However,some things are already clear. First,the method is not necessarily harder to program than finite difference or finite element algorithms.In Maple,the complete solution of the ODE/BVP takes just five lines (Table 1.1)! Second,spectral methods are not purely numerical.When N is sufficiently small, Chebyshev and Fourier methods yield an analytic answer. 1.3 Comparison with finite element methods Finite element methods are similar in philosophy to spectral algorithms;the major differ- ence is that finite elements chop the interval in z into a number of sub-intervals,and choose the o(z)to be local functions which are polynomials of fixed degree which are non-zero only over a couple of sub-intervals.In contrast,spectral methods use global basis func- tions in which (z)is a polynomial (or trigonometric polynomial)of high degree which is non-zero,except at isolated points,over the entire computational domain. When more accuracy is needed,the finite element method has three different strategies. The first is to subdivide each element so as to improve resolution uniformly over the whole domain.This is usually called "h-refinement"because h is the common symbol for the size or average size of a subdomain.(Figure 1.2). The second alternative is to subdivide only in regions of steep gradients where high resolution is needed.This is "r-refinement". The third option is to keep the subdomains fixed while increasing p,the degree of the polynomials in each subdomain.This strategy of"p-refinement"is precisely that employed by spectral methods.Finite element codes which can quickly change p are far from univer- sal,but those that can are some called "p-type"finite elements. Finite elements have two advantages.First,they convert differential equations into matrix equations that are sparse because only a handful of basis functions are non-zero in a given sub-interval.("Sparse"matrices are discussed in Appendix B;suffice it to say that "sparse"matrix equations can be solved in a fraction of the cost of problems of similar size with "full"matrices.)Second,in multi-dimensional problems,the little sub-intervals become little triangles or tetrahedra which can be fitted to irregularly-shaped bodies like the shell of an automobile.Their disadvantage is low accuracy (for a given number of degrees of freedom N)because each basis function is a polynomial of low degree. Spectral methods generate algebraic equations with full matrices,but in compensation, the high order of the basis functions gives high accuracy for a given N.When fast iterative matrix-solvers are used,spectral methods can be much more efficient than finite element
4 CHAPTER 1. INTRODUCTION Table 1.1: Maple program to solve linear boundary-value problem u2:=1 + (1-x*x)*(a0 + a1*x + a2*x*x); Resid:= diff(u,x,x) - (x**6 + 3*x**2)*u; eq1:=subs(x=-1/2,Resid); eq2:=subs(x=0,Resid); eq3:=subs(x=1/2,Resid); solutionarray:=solve({eq1,eq2,eq3},{a0,a1,a2}); The fifth answer is: the algebraic equations can be written (for a linear differential equation) as a matrix equation, which can then be solved by library software in FORTRAN or C. Many other questions will be asked and answered in later chapters. However, some things are already clear. First, the method is not necessarily harder to program than finite difference or finite element algorithms. In Maple, the complete solution of the ODE/BVP takes just five lines (Table 1.1)! Second, spectral methods are not purely numerical. When N is sufficiently small, Chebyshev and Fourier methods yield an analytic answer. 1.3 Comparison with finite element methods Finite element methods are similar in philosophy to spectral algorithms; the major difference is that finite elements chop the interval in x into a number of sub-intervals, and choose the φn(x) to be local functions which are polynomials of fixed degree which are non-zero only over a couple of sub-intervals. In contrast, spectral methods use global basis functions in which φn(x) is a polynomial (or trigonometric polynomial) of high degree which is non-zero, except at isolated points, over the entire computational domain. When more accuracy is needed, the finite element method has three different strategies. The first is to subdivide each element so as to improve resolution uniformly over the whole domain. This is usually called “h-refinement” because h is the common symbol for the size or average size of a subdomain. (Figure 1.2). The second alternative is to subdivide only in regions of steep gradients where high resolution is needed. This is “r-refinement”. The third option is to keep the subdomains fixed while increasing p, the degree of the polynomials in each subdomain. This strategy of ”p-refinement” is precisely that employed by spectral methods. Finite element codes which can quickly change p are far from universal, but those that can are some called “p-type” finite elements. Finite elements have two advantages. First, they convert differential equations into matrix equations that are sparse because only a handful of basis functions are non-zero in a given sub-interval. (“Sparse” matrices are discussed in Appendix B; suffice it to say that “sparse” matrix equations can be solved in a fraction of the cost of problems of similar size with “full” matrices.) Second, in multi-dimensional problems, the little sub-intervals become little triangles or tetrahedra which can be fitted to irregularly-shaped bodies like the shell of an automobile. Their disadvantage is low accuracy (for a given number of degrees of freedom N) because each basis function is a polynomial of low degree. Spectral methods generate algebraic equations with full matrices, but in compensation, the high order of the basis functions gives high accuracy for a given N. When fast iterative matrix–solvers are used, spectral methods can be much more efficient than finite element