Fourier Analysis

Chapter 6 Fourier analysis (Historical intro: the heat equation on a square plate or interval.) Fourier’s analysis was tremendously successful in the 19th century for formulating series expansions for solutions of some very simple ODE and PDE. This class shows that in the 20th century, Fourier analysis has established itself as a central tool for numerical computations as well, for vastly more general ODE and PDE when explicit formulas are not available. 6.1 The Fourier transform We will take the Fourier transform of integrable functions of one variable x ∈ R. Definition 13. (Integrability) A function f is called integrable, or absolutely integrable, when Z ∞ |f (x)| dx < ∞, −∞ in the sense of Lebesgue integration. One also writes f ∈ L1 (R) for the space of integrable functions. We denote the physical variable as x, but it is sometimes denoted by t in contexts in which its role is time, and one wants to emphasize that. The frequency, or wavenumber variable is denoted k. Popular alternatives choices for the frequency variable are ω (engineers) or ξ (mathematicians), or p (physicists). 1 CHAPTER 6. FOURIER ANALYSIS Definition 14. The Fourier transform (FT) of an integrable function f (x) is defined as Z ∞ ˆ f (k) = e−ikx f (x) dx. (6.1) −∞ When fˆ(k) is also integrable, f (x) can be recovered from fˆ(k) by means of the inverse Fourier transform (IFT) Z ∞ 1 eikx fˆ(k) dk. (6.2) f (x) = 2π −∞ Intuitively, fˆ(k) is the amplitude density of f at frequency k. The formula for recovering f is a decomposition of f into constituent waves. The justification of the inverse FT formula belongs in a real analysis class (where it is linked to the notion of approximate identity.) We will justify the form of (6.2) heuristically when we see Fourier series in the next section. The precaution of assuming integrability is so that the integrals can be understood in the usual Lebesgue sense. In that context, taking integrals over infinite intervals is perfectly fine. If (6.1) and (6.2) are understood as limits of integrals over finite intervals, it does not matter how the bounds are chosen to tend to ±∞. One may in fact understand the formulas for the FT and IFT for much larger function classes than the integrable functions, namely distributions, but this is also beyond the scope of the class. We will generally not overly worry about these issues. It is good to know where to draw the line: the basic case is that of integrable functions, and anything beyond that requires care and adequate generalizations. Do not be surprised to see alternative formulas for the Fourier transfom in other classes or other contexts. Wikipedia lists them. Here are some important properties of Fourier transforms: • (Differentiation) fb0 (k) = ikfˆ(k). Justification: integration by parts in the integral for the FT. • (Translation) If g(x) = f (x + a), then ĝ(k) = eika fˆ(k). Justification: change of variables in the integral for the FT. 2 6.1. THE FOURIER TRANSFORM Let’s see some examples of FT. Example 17. Let 1 f (x) = χ[−a,a] (x) = 2a Then 1 fˆ(k) = 2a Z  1 2a if x ∈ [−a, a]; otherwise. 0 a e−ikx dx = −a sin(ka) . ka This function is a scaled version of the sinc function, sinc(k) = sin k . k It is easy to check by L’Hospital’s rule that sinc(0) = 1. At k → ∞, sinc(k) decays like 1/k, but does so by alternating between positive and negative values. It is a good exercise to check that sinc is not absolutely integrable. It turns out that the Fourier transform can still be defined for it, so lack of integrability is not a major worry. Example 18. Consider the Gaussian function f (x) = e−x 2 /2 . By completing the square and adequately modifying the contour of integration in the complex plane (not part of the material for this class), it can be shown that √ 2 fˆ(k) = 2π e−k /2 . Example 19. The Dirac delta δ(x) has a FT equal to 1 (why?). Another basic property of Fourier transforms is the convolution theorem. Theorem 11. (The convolution theorem.) Denote convolution as f ? g(x) = R∞ f (y)g(x − y) dy. Then −∞ f[ ? g(k) = fˆ(k) ĝ(k). 3 CHAPTER 6. FOURIER ANALYSIS Proof. Let h = f ? g. We can use Fubini below provided every function is integrable. Z Z −ikx ĥ(k) = e f (y)g(x − y) dydx Z Z = e−iky f (y)e−ik(x−y)) g(x − y) dydx Z  Z  −iky −ikx0 0 0 = e f (y) dy e g(x ) dx = fˆ(k) ĝ(k). The Fourier transform is an important tool in the study of linear differential equations because it turns differential problems into algebraic problems. P For instance, consider a polynomial P (x) = an xn , and the ODE   d P u(x) = f (x), x ∈ R, dx P dn u which means an dxn = f . (Such ODE are not terribly relevant in real life because they are posed over the whole real line.) Upon Fourier transformation, the equation becomes P (ik)û(k) = fˆ(k), which is simply solved as û(k) = fˆ(k) , P (ik) and then Z ∞ 1 fˆ(k) u(x) = eikx dk. 2π −∞ P (ik) Beware the zeros of P when applying this formula! They always carry important physical interpretation. For instance, they could be resonances of a mechanical system. ˆ(k) The formula û(k) = Pf(ik) also lends itself to an application of the convolution theorem. Let K(x) be the inverse Fourier transform of 1/P (ik). Then we have Z u(x) = K(x − y)f (y) dy. The function K is called Green’s function for the original ODE. 4 6.2. SAMPLING AND RESTRICTION 6.2 Sampling and restriction We aim to use Fourier transforms as a concept to help understand the accuracy of representing and manipulating functions on a grid, using a finite number of degrees of freedom. We also aim at using a properly discretized Fourier transform as a numerical tool itself. For this purpose, x ∈ R and k ∈ R must be replaced by x and k on finite grids. Full discretization consists of sampling and restriction. Let us start by sampling x ∈ hZ, i.e., considering xj = jh for j ∈ Z. The important consequence of sampling is that some complex exponential waves eikx for different k will appear to be the same on the grid xj . We call aliases such functions that identify on the grid. Definition 15. (Aliases) The functions eik1 x and eik2 x are aliases on the grid xj = jh if eik1 xj = eik2 xj , ∀j ∈ Z. Aliases happen as soon as k1 jh = k2 jh + 2π × integer(j). Letting j = 1, and calling the integer n, we have k1 − k2 = 2π n, h for some n ∈ Z. To wave numbers k1 , k2 are indistinguishable on the grid if they differ by an integer multiple of 2π/h. For this reason, we restrict without loss of generality the wavenumber to the interval k ∈ [−π/h, π/h]. We also call this interval the fundamental cell in frequency (in reference to a similar concept in crytallography.) Real-life examples of aliases are rotating wheels looking like they go backwards in a movie, Moiré patterns on jackets on TV, and stroboscopy. The proper notion of Fourier transform on a grid is the following. Definition 16. Let xj = hj, fj = f (xj ). Semidiscrete Fourier transform (SFT): ∞ X fˆ(k) = h e−ikxj fj , k ∈ [−π/h, π/h]. (6.3) j=−∞ 5 CHAPTER 6. FOURIER ANALYSIS Inverse semidiscrete Fourier transform (ISFT): Z π/h 1 fj = eikx fˆ(k) dk. 2π −π/h (6.4) As we saw, sampling in x corresponds to a restriction in k. If one still wanted to peek outside [−π/h, π/h] for the SFT, then the SFT would simply repeat by periodicity: 2nπ ) = fˆ(k). fˆ(k + h (why?). That’s why we restrict it to the fundamental cell. We can now define the proper notion of Fourier analysis for functions that are restricted to x in some interval, namely [−π, π] for convention. Unsurprisingly, the frequency is sampled as a result. the following formulas are dual to those for the SFT. Definition 17. Fourier series (FS): Z π ˆ e−ikx f (x) dx. fk = (6.5) −π Inverse Fourier series (IFS) ∞ 1 X ikx ˆ f (x) = e fk , 2π k=−∞ x ∈ [−π, π]. (6.6) If one uses the Fourier series inversion formula (6.6) for x outside of its intended interval [−π, π], then the function simply repeats by periodicity: f (x + 2nπ) = f (x). (again, why?) The twoR formulas (6.5) and (6.6) can be justified quite intuitively. The expression f (x)g(x) dx is an inner product on functions. It is easy to see that the complex exponentials r h −ikxj e = vj (k) 2π form an orthonormal set of functions on [−π/h, π/h], for this inner product. Hence, up to normalization constants, (6.5) is simply calculation of the coefficients in an orthobasis (analysis), and (6.6) is the synthesis of the function 6 6.3. THE DFT AND ITS ALGORITHM, THE FFT back from those coefficients. We’d have to understand more about the peculiarities of infinite-dimensional linear algebra to make this fully rigorous, but this is typically done in a real analysis class. Here’s an example of SFT. Example 20. fj = 1 χ[−a,a] (xj ) 2a Then a h X −ikjh fˆ(k) = e 2a j=−a 2a h ikah X −ikjh = e e 2a j=0 h ikah (e−ikh )2a+1 − 1 e 2a e−ikh − 1 h sin(kh(a + 1/2)) = . 2a sin(kh/2) = (geometric series) This function is called the discrete sinc. It looks like a sinc, but it periodizes smoothly when k = −π/h and k = π/h are identified. Our first slogan is therefore: Sampling in x corresponds to restriction/periodization in k, and restriction/periodization in k corresponds to sampling in x. 6.3 The DFT and its algorithm, the FFT The discrete Fourier transform is what is left of the Fourier transfom when both space and frequency are sampled and restricted to some interval. Consider xj = jh, j = 1, . . . , N. The point j = 0 is identified with j = N by periodicity, so it is not part of the grid. If the endpoints are x0 = 0 and xN = 2π, then N and h relate as h= 2π N ⇒ π N = . h 2 7 CHAPTER 6. FOURIER ANALYSIS For the dual grid in frequency, consider that N points should be equispaced between the bounds [−π/h, π/h]. The resulting grid is k=− k, N N + 1, . . . , . 2 2 We have the following definition. Definition 18. Discrete Fourier transform (DFT): fˆk = h N X e−ikjh fj , N N ,..., . 2 2 (6.7) j = 1, . . . , N. (6.8) k=− j=1 Inverse discrete Fourier transform (IDFT) 1 fj = 2π N/2 X eikjh fˆk , k=−N/2+1 The DFT can be computed as is, by implementing the formula (6.7) directly on a computer. The complexity of this calculation is a O(N 2 ), since there are N values of j, and there are N values of k over which the computation must be repeated. There is, however, a smart algorithm that allows to group the computation of all the fk in complexity O(N log N ). It is called the fast Fourier transform (FFT). It is traditionally due to Tukey and Cooley (1965), but the algorithm had been discovered a few times before that by people who are not usually credited as much: Danielson and Lanczos in 1942, and well as Gauss in 1805. The trick of the FFT is to split the samples of f into even samples (j even) and odd samples (j odd). Assume that N is a power of 2. Then fˆk = h N X e−ikjh fj j=1 N/2 =h X −ik(2j)h e f2j + h j=1 =h N/2 X N/2 X e−ik(2j+1)h f2j+1 j=1 −ikj(2h) e ikh f2j + he j=1 N/2 X j=1 8 e−ikj(2h) f2j+1 . 6.4. SMOOTHNESS AND TRUNCATION The first term on the last line is simply the DFT of length N/2, on a grid of spacing 2h, of the even samples of f . The second term is, besides the multiplicative factor, the DFT of length N/2, on a grid of spacing 2h, of the odd samples of f . Note that those smaller DFTs would normally be only calculated for k = − N4 + 1, . . . , N4 , but we need them for k = − N2 + 1, . . . , N2 . This is not a big problem: we know that the DFT extends by periodicity outside the standard bounds for k, so all there is to do is copy fˆk by periodicity outside of k = − N4 + 1, . . . , N4 . Already, we can see the advantage in this reduction: solving one problem of size N is essentially reduced to solving two problems of size N/2. Even better, the splitting into even and odd samples can be repeated recursively until the DFT are of size 1. When N = 1, the DFT is simply multiplication by a scalar. At each stage, there are O(N ) operations to do to put together the summands in the equation of the last line above. Since there are O(log N ) levels, the overall complexity is O(N log N ). There are variants of the FFT when N is not a power of 2. 6.4 Smoothness and truncation In this section, we study the accuracy of truncation of Fourier transforms to finite intervals. This is an important question not only because reallife numerical Fourier transforms are restricted in k, but also because, as we know, restriction in k serves as a proxy for sampling in x. It will be apparent in Chapter 2, section 2.1 that every claim that we make concerning truncation of Fourier transforms will have an implication in terms of accuracy of sampling a function on a grid, i.e., how much information is lost in the process of sampling a function f (x) at points xj = jh. We will manipulate functions in the spaces L1 , L2 , and L∞ . We have already encountered L1 . Definition 19. Let 1 ≤ p < ∞. A function f of x ∈ R is said to belong to the space Lp (R) when Z ∞ |f (x)|p dx < ∞. −∞ p Then the norm of f in L (R) is kf kp = R 9 ∞ −∞ p |f (x)| dx 1/p . CHAPTER 6. FOURIER ANALYSIS A function f of x ∈ R is said to belong to L∞ when ess sup |f (x)| < ∞. Then the norm of f in L∞ (R) is ess sup |f (x)|. In the definition above, “ess sup” refers to the essential supremum, i.e., the infimum over all dense sets X ⊂ R of the supremum of f over X. A set X is dense when R\X has measure zero. The notions of supremum and infimum correspond to maximum and minimum respectively, when they are not necessarily attained. All these concepts are covered in a real analysis class. For us, it suffices to heuristically understand the L∞ norm as the maximum value of the modulus of the function, except possibly for isolated points of discontinuity which don’t count in calculating the maximum. It is an interesting exercise to relate the L∞ norm to the sequence of Lp norms as p → ∞. We will need the very important Parseval and Plancherel identitites. They express “conservation of energy” from the physical domain to the frequency domain. Theorem 12. (Parseval’s identity). Let f, g ∈ L1 (R) ∩ L2 (R). Then Z ∞ Z ∞ 1 f (x)g(x) dx = fˆ(k)ĝ(k) dk. 2π −∞ −∞ Proof. Let h be the convolution f ? g̃, where g̃(x) = g(−x). It is easy to see that the Fourier transform of g̃ is ĝ(k) (why?). By the convolution theorem (Section 1.1), we have ˆ h(k) = fˆ(k) ĝ(k). If we integrate this relation over k ∈ R, and divide by 2π, we get the IFT at x = 0: Z ∞ 1 h(0) = fˆ(k)ĝ(k) dk. 2π −∞ On the other hand, Z ∞ Z ∞ f (x)g(−(0 − x)) dx = h(0) = −∞ f (x)g(x) dx. −∞ 10 6.4. SMOOTHNESS AND TRUNCATION Theorem 13. (Plancherel’s identity). Let f ∈ L1 (R) ∩ L2 (R). Then Z ∞ Z ∞ 1 2 |fˆ(k)|2 dk. |f (x)| dx = 2π −∞ −∞ Proof. Apply Parseval’s identity with g = f . (With the help of these formulas, it is in fact possible to extend their validity and the validity of the FT to f, g ∈ L2 (R), and not simply f, g ∈ L1 (R) ∩ L2 (R). This is a classical density argument covered in many good analysis texts.) We need one more concept before we get to the study of truncation of Fourier transforms. It is the notion of total variation. We assume that the reader is familiar with the spaces C k (R) of bounded functions which are k times continuously differentiable. Definition 20. (Total variation) Let f ∈ C 1 (R). The total variation of f is the quantity Z ∞ |f 0 (x)| dx. (6.9) kf kT V = −∞ 1 For functions that are not C , the notion of total variation is given by either expression Z ∞ X |f (x) − f (x − h)| kf kT V = lim dx = sup |f (xp+1 )−f (xp )|, h→0 −∞ |h| {xp } finite subset of R p (6.10) These more general expressions reduce to −∞ |f (x)| dx when f ∈ C 1 (R). When a function has finite total variation, we say it is in the space of functions of bounded variation, or BV(R). R∞ 0 The total variation of a piecewise constant function is simply the sum of the absolute value of the jumps it undergoes. This property translates to a useful intuition about the total variation of more general functions if we view them as limits of piecewise constant functions. The important meta-property of the Fourier transform is that decay for large |k| corresponds to smoothness in x. 11 CHAPTER 6. FOURIER ANALYSIS There are various degrees to which a function can be smooth or rates at which it can decay, so therefore there are several ways that this assertion can be made precise. Let us go over a few of them. Each assertion either expresses a decay (in k) to smoothness (in x) implication, or the converse implication. • Let fˆ ∈ L1 (R) (decay), then f ∈ L∞ (R) and f is continuous (smoothness). That’s because |eikx | = 1, so Z Z 1 1 ikx ˆ |f (x)| ≤ |e f (k)| dk = |fˆ(k)| dk, 2π 2π which proves boundedness. As for continuity, consider a sequence yn → 0 and the formula Z 1 f (x − yn ) = eik(x−yn ) fˆ(k) dk. 2π The integrand converges pointwise to eikx fˆ(k), and is uniformly bounded in modulus by the integrable function |fˆ(k)|. Hence Lebesgue’s dominated convergence theorem applies and yields f (x − yn ) → f (x), i.e., continuity in x. • Let fˆ(k)(1 + |k|p ) ∈ L1 (R) (decay). Then f ∈ C p (smoothness). We saw the case p = 0 above; the justification is analogous in the general case. We write Z Z 1 (n) ikx nˆ |f (x)| ≤ |e (ik) f (k)| dk ≤ |k|n |fˆ(k)| dk, 2π which needs to be bounded for all 0 ≤ n ≤ p. This is obviously the case if fˆ(k)(1 + |k|p ) ∈ L1 (R). Continuity of f (p) is proved like before. • Let f ∈ BV(R) (smoothness). Then fˆ(k) ≤ kf kT V |k|−1 (decay). If f ∈ C 1 ∩ BV (R), then this is justified very simply from (6.9), and Z ˆ ikf (k) = e−ikx f 0 (x) dx. Take a modulus on both sides, and get the desired relation Z ˆ |k||f (k)| ≤ |f 0 (x)|dx = kf kT V < ∞. 12 6.4. SMOOTHNESS AND TRUNCATION When f ∈ BV (R), but f ∈ / C 1 , either of the more general formulas (6.10) must be used instead. It is a great practice exercise to articulate a modified proof using the limh→0 formula, and properly pass to the limit. • Let f such that f (k) ∈ L2 (R) for 0 ≤ k < p, and assume that f (p) ∈ BV(R) (smoothness). Then there exists C > 0 such that |fˆ(k)| ≤ |k|−p−1 (decay). This claim is also formulated as point (a) in Theorem 1 on page 30 of Trefethen’s “Spectral methods in Matlab”. The justification is very simple when f ∈ C p+1 : we then get Z p+1 ˆ (ik) f (k) = e−ikx f (p+1) (x) dx, so p+1 |k| |fˆ(k)| ≤ Z |f (p+1) (x)| dx = kf (p) kT V < ∞. Again, it is a good exercise to try and extend this result to functions not in C p+1 . • (This is the one we’ll use later). (Same proof as above.) In summary, let f have p derivatives in L1 . Then |fˆ(k)| ≤ C|k|−p . This is the form we’ll make the most use of in what follows. Example 21. The canonical illustrative example for the two statements involving bounded variation is that of the B-splines. Consider 1 s(x) = χ[−1,1] (x), 2 the rectangle-shaped indicator of [−1, 1] (times one half ). It is a function in BV (R), but it has no derivative in L2 (R). Accordingly, its Fourier transform ŝ(k) is predicted to decay at a rate ∼ |k|−1 . This is precisely the case, for we know sin k ŝ(k) = . k Functions with higher regularity can be obtained by auto-convolution of s; for instance s2 = s ? s is a triangle-shaped function which has one derivative in 13 CHAPTER 6. FOURIER ANALYSIS L2 , and such that s02 ∈ BV (R). We anticipate that sˆ2 (k) would decay like |k|−2 , and this the case since by the convolution theorem 2 sˆ2 (k) = (ŝ(k)) =  sin k k 2 . Any number of autoconvolutions s ? s . . . ? s can thus be considered: that’s the family of B-splines. The parallel between smoothness in x and decay in k goes much further. We have shown that p derivatives in x very roughly corresponds to an inversepolynomial decay |k|−p in k. So C ∞ functions in x have so called superalgebraic decay in k: faster than Cp |k|−p for all p ≥ 0. The gradation of very smooth functions goes beyond C ∞ functions to include, in order: • analytic functions that can be extended to a strip in the complex plane (like f (x) = 1/(1 + x2 )), corresponding to a Fourier transform decaying exponentially (in this case fˆ(k) = πe−|k| ). That’s Paley-Wiener theory, the specifics of which is not material for this course. • analytic functions that can be extended to the whole complex plane 2 with super-exponential growth (like f (x) = e−x /2 ), whosepFourier 2 transform decays faster than exponentially (in this case fˆ(k) = π/2e−k /2 ). • analytic functions that can be extended to the whole complex plane with exponential growth at infinity (like f (x) = sinx x ), whose Fourier transform is compactly supported (in this case fˆ(k) = 2πχ[−1,1] (k)). That’s Paley-Wiener theory in reverse. Such functions are also called bandlimited. More on this in Chapter 4 of Trefethen’s book. An important consequence of a Fourier transform having fast decay is that it can be truncated to some large interval k ∈ [−N, N ] at the expense of an error that decays fast as a function of N . The smoother f , the faster fˆ decays, and the more accurate the truncation to large N in frequency. On the other hand, there may be convergence issues if f (x) is not smooth. To make this quantitative, put fˆN (k) = χ[−N,N ] (k)fˆ(k). 14 6.4. SMOOTHNESS AND TRUNCATION Recall that have 1 2π RN −N eikx dx = sin(N x) . πx By the convolution theorem, we therefore sin N x ? f (x). πx In the presence of f ∈ L2 (R), letting N → ∞ always gives rise to convergence fN → f in L2 (R). This is because by Plancherel’s identity, Z ∞ Z 1 1 2 2 ˆ ˆ kf − fN k2 = |fN (k) − f (k)| dk = |fˆ(k)|2 dk . 2π −∞ 2π |k|>N fN (x) = R This quantity tends to zero as N → ∞ since the integral R |fˆ(k)|2 over the whole line is bounded. The story is quite different if we measure convergence in L∞ instead of 2 L , when f has a discontinuity. Generically, L∞ (called uniform) convergence fails in this setting. This phenomenon is called Gibb’s effect, and manifests itself through ripples in reconstructions from truncated Fourier transforms. Let us illustrate this on the Heaviside step function  1 if x ≥ 0; u(x) = 0 if x < 0. (It is not a function in L2 (R) but this is of no consequence to the argument. It could be modified into a function in L2 (R) by some adequate windowing.) Since u(x) is discontinuous, we expect the Fourier transform to decay quite slowly. Consider the truncated FT ûN (k) = χ[−N,N ] (k)û(k). Back in the x domain, this is sin N x ? u(x) πx Z ∞ sin N (x − y) = dy π(x − y) 0 Z Nx sin y = dy −∞ πy ≡ s(N x). uN (x) = (Draw picture) 15 CHAPTER 6. FOURIER ANALYSIS The function s(N x) is a sigmoid function going from 0 at −∞ to 1 at ∞, going through s(0) = 1/2, and taking on min. value ' −.045 at x = −π/N and max. value ' 1.045 at x = π/N . The oscillations have near-period π/N . The parameter N is just a dilation index for the sigmoid, it does not change the min. and max. values of the function. Even as N → ∞, there is no way that the sigmoid would converge to the Heaviside step in a uniform manner: we always have ku − uN k∞ & .045. This example of Gibb’s effect goes to show that truncating Fourier expansions could be a poor numerical approximation if the function is nonsmooth. However, if the function is smooth, then everything goes well. Let us study the decay of the approximation error 2N = kf − fN k22 for truncation of the FT to [−N, N ]. Again, we use Plancherel to get Z 2 |fˆ(k)|2 dk N = |k|>N Now assume that fˆ(k) ≤ C|k |−p , a scenario already considered earlier. Then Z 2 N ≤ C |k|−2p ≤ C 0 N −2p+1 , |k|>N 00 so, for some constant C (dependent on p but independent of N , hence called constant), N ≤ C 00 N −p+1/2 . The larger p, the faster the decay as N → ∞. When a function is C ∞ , the above decay rate of the approximation error is valid for all p > 0. When the function is analytic and the Fourier transform decays exponentially, the decay rate of the approximation error is even faster, itself exponential (a good exercise). In numerical analysis, either such behavior is called spectral accuracy. 6.5 Chebyshev expansions Chebyshev expansions, interpolation, differentiation, and quadrature rules are not part of the material for 18.330 in 2012. 16 6.5. CHEBYSHEV EXPANSIONS In the previous sections, we have seen that smooth functions on the real line have fast decaying Fourier transforms. On a finite interval, a very similar property hold for Fourier series: if a function is smooth in [−π, π] and connects smoothly by periodicity at x = −π and x = π, then its Fourier series (6.5) decays fast. Periodicity is essential, because the points x = −π and x = π play no particular role in (6.5). They might as well be replaced by x = 0 and x = 2π for the integration bounds. So if a function is to qualify as smooth, it has to be equally smooth at x = 0 as it is at x = −π, identified with x = π by periodicity. 6 For instance if a function f is smooth inside [−π, π], but has f (−π) = f (π), then for the purpose of convergence of Fourier series, f is considered discontinuous. We know what happens in this case: the Gibbs effect takes place, partial inverse Fourier series have unwelcome ripples, and convergence does not occur in L∞ . If f is smooth and periodic, then it is a good exercise to generalize the convergence results of the previous section, from the Fourier transform to Fourier series. How shall we handle smooth functions in intervals [a, b], which do not connect smoothly by periodicity? The answer is not unique, but the most standard tool for this in numerical analysis are the Chebyshev polynomials. For simplicity consider x ∈ [−1, 1], otherwise rescale the problem. Take a C k nonperiodic f (x) in [−1, 1]. The trick is to view it as g(θ) = f (cos θ). Since cos[0, π] = cos[π, 2π] = [−1, 1], all the values of x ∈ [−1, 1] are covered twice by θ ∈ [0, 2π]. Obviously, at any point θ, g inherits the smoothness of f by the chain rule. Furthermore, g is periodic since cos θ is periodic. So g is exactly the kind of function which we expect should have a fast converging Fourier series: Z ĝk = 2π −ikθ e g(θ) dθ, 0 ∞ 1 X ikθ g(θ) = e ĝk , 2π k=−∞ k ∈ Z. Since g(θ) is even in θ, we may drop the sin θ terms in the expansion, as well as the negative k: Z gˆk = 2π cos (kθ)g(θ) dθ, 0 ∞ 1 1X g(θ) = ĝ0 + cos (kθ)ĝk , 2π π k=1 17 k ∈ Z+ . CHAPTER 6. FOURIER ANALYSIS Back to f , we find Z 1 dx cos (k arccos x)f (x) √ , ĝk = 2 1 − x2 −1 ∞ 1 1X f (x) = ĝ0 + cos (k arccos x)ĝk , 2π π k=1 The function cos(k arccos x) happens to be a polynomial in x, of order k, called the Chebyshev polynomial of order k. Switch to the letter n as is usually done: Tn (x) = cos(n arccos x), x ∈ [−1, 1]. The first few Chebyshev polynomials are T0 (x) = 1, T1 (x) = x, T2 (x) = 2x2 − 1. It is a good exercise (involving trigonometric identities) to show that they obey the recurrence relation Tn+1 (x) = 2xTn (x) − Tn−1 (x). The orthogonality properties of Tn (x) over [−1, 1] follows from those of cos(nθ) over [0, π], but watch the special integration weight: Z 1 dx Tm (x)Tn (x) √ = cn δmn , 1 − x2 −1 with cn = π/2 if n = 0 and cn = π otherwise. The Chebyshev expansion of f is therefore Z 1 hf, Tn i = −1 Tn (x)f (x) √ dx , 1 − x2 f (x) = ∞ 1 2X hf, T0 i+ Tn (x)hf, Tn i, π π n=1 Under the hood, this expansion is the FS of a periodic function, so we can apply results pertaining to Fourier series and Fourier transforms to obtain fast decay of Chebyshev expansions. This way, spectral accuracy is restored for nonperiodic functions defined on an interval. 18 k ∈ Z+ . k ∈ Z+