IV. Calculating Fourier Series

Linear Methods of Applied Mathematics

Evans M. Harrell II and James V. Herod*

*(c) Copyright 1994,1997,2000 by Evans M. Harrell II and James V. Herod. All rights reserved.

version of 8 May 2000

If you wish to print a nicely formatted version of this chapter, you may download the rtf file, which will be interpreted and opened by Microsoft Word or the pdf file, which will be interpreted and opened by Adobe Acrobat Reader.

(Some remarks for the instructor).

IV. Calculating Fourier Series

In this section we calculate several Fourier series. As we know, to find a Fourier series simply means calculating various integrals, which can often be done with software or with integral tables. Since performing integrals is not much more interesting in the modern age than long division, our goals in this section will be to get a visual and analytic impression of what to expect from Fourier series, and to understand the rôle of symmetry in the calculations. Many of the calculations in this chapter are available in a Mathematica notebook or Maple worksheet.

Let's begin by evaluating the Fourier series for the functions:

f(x) = 1 for 0 <= x < L/2, but 0 for L/2 <= x <= L

and

g(x) = x, 0 <=x < L.

The functions have not been defined at the points of discontinuity, but as we know, the Fourier series will converge there to the average of the limit from the left and the limit from the right. The end point L is essentially a jump point, because the periodic extension of the functions make the values x=L and x=0 equivalent.

Here is a graph of the function f, called a "square pulse" or "square wave" (when extended periodically): (graph of square wave)

The length L has been chosen as pi.

We want to represent these functions in the form

a[0] + Sum[a[m] Cos[2 Pi m x/L, {m,1,N}] + Sum[b[n] Sin[2 Pi n x/L, {n,1,N}], beginning with f(x).

Model Problem IV.1. Use software to calculate the full Fourier series for the function f(x) as defined above, and investigate its convergence.

Solution. The formulae for these coefficients were given in Chapter II:

a_0 = (1/L) integral 0 to L of f(x). (i.e., the average of f), and for the other coefficients: a_m = (2/L) integral 0 to L of f(x)*cos(2 Pi m x/L)
and
b_n = (2/L) integral 0 to L of f(x)*sin(2 Pi n x/L). Here are the results of the calculation (see notebook or worksheet): a0 = 1/2 ;
am = 0 for m = 1, 2, ...;
bn = (1 - (-1)n)/(npi) A point to ponder: why are all the am = 0? Surely this is no coincidence!

The sort of expression we got for the bn simplifies in a way which arises often in Fourier series for elementary functions: If n=2k is even, b2k = 0. Thus only odd terms survive, in which case

b2k+1 = 2/(pi(2k+1)). We can plot the Fourier series and watch it converge to the original function as more and more terms are included (see notebook or worksheet for the calculations):

With one sine term:

(graph of square wave and approximation)

With three sine terms:

(graph of square wave and approximation)

With seven sine terms:

(graph of square wave and approximation)

And so on. We call attention to the systematic overshoot that occurs at the edges of the jump. Curiously, the size of this last overshoot does not tend to 0 as we include more terms. (You can see this in an animation prepared for us by Georgia Tech student Ning Wu; the fact that Ning's function has vertical range -1 to 1 is not significant for this purpose.) The bumps next to the jump get thinner rather than shorter. This is known as the Gibbs phenomenon.

Model Problem IV.2. We can now investigate some questions about the excitation of mechanical resonances. Suppose that experiments by K. Battle at the Wiener Staatsoper in Austria show that a crystal goblet will shatter if the intensity of the tone at frequency 1760 Hz (high A) exceeds .01. We use physical units in which the proportionality between the power and the square of the amplitude is 1, i.e., we define the intensity simply as

I = ak2 + bk2. At Georgia Tech, not a particularly musical spot, we can generate a square pulse with amplitude A for 1/704 sec., then amplitude 0 for the same length of time, then A, etc., periodically with period 1/352.

Question: What amplitude A will cause the glass to shatter in our laboratory?

Solution. The issue is essentially to estimate the magnitude of the Fourier coefficient for our function corresponding to the frequency 1760 Hertz.

The mapping from functions to their Fourier coefficients is linear. Hence the coefficients can be obtained from the ones we just calculated.

Step 1. Scale the independent variable by replacing x with t and L with 1/352.

Step 2. Multiply all the coefficients by A.

The harmonic with frequency 1760 is the one with n=5, so the Fourier coefficients corresponding to this frequency are

a5 = 0, and
b5 = 2A/(5 pi) The intensity is 4 A2/(25 pi 2). The amplitude A above which the glass will shatter is pi /4, numerically about 0.78539816.

Model Problem IV.3. For comparison, let us find another Fourier series, namely the one for the periodic extension of g(x) = x, 0 <=x <= 1, sometimes designated x mod 1. Watch it converge.

Solution. (For more details on the calculations, see the Mathematica notebook or the Maple worksheet. For x between 1 and 2, the function is (x-r1L), for x between 2 and 3 it is (x-2), etc. Its graph has a sawtooth shape:

(graph of 'sawtooth' wave)

Because the interval has length 1, the Fourier basis functions are 1, cos(2npi x), and sin(2npi x).

The Fourier coefficients are calculated as:

a0 = 1/2 ;
am = 0 for m = 1, 2, ...;
bn = (-1)/(n pi) The coefficients am are all zero, again. Why?

Let us get an impression of how the series converges, by plotting the contributions of the first two sines, and then of six sines: (graph of sawtooth wave and approximation)
(graph of sawtooth wave and approximation)

We also recall that on a longer interval, the Fourier series produces a periodic function:
(graph of sawtooth wave and approximation)

The Fourier series is converging nicely to the function except at the end-points of the interval, which are places where the full periodic saw-tooth function has jumps, and we saw something similar with the square pulse. In both cases we see numerical evidence for the theorem that the Fourier series converges to f(x) where f(x) is continuous, and where it has a jump, the Fourier series converges to the average of the upper and the lower value at the jump.

Let's try out a different sort of example, where we need to integrate numerically (see notebook).

Model Problem IV.4. Find the Fourier series for the function sin(pi/x).

Solution. The integrals called for are not in the integral tables. Indeed, most  integrals are not in the integral tables, even when they are nowhere near as wild as this function (ask yourself what happens near x=0). In the information age this is no barrier. We simply call on the software to do the integrals numerically. Since the function sin(pi/x) is odd in x, the coefficients associated with even functions, am, will all be zero. With some complaints about the convergence, the software calculates the first six sine coefficients as: 


             -0.3071329778789654123
              0.5062509679699521912
             -0.271449220623610834
             -0.2561291634161970945
              0.1947511859947028479
              0.2192911372159977088
Let us plot the Fourier series with these six terms and compare with the function itself:

(graph of approximation of sin (Pi/x))
(graph of sin (Pi/x))

Notice that the convergence is pretty good away from the nasty singularity at x=0. In effect, truncating the Fourier series has filtered out the high-frequency oscillations.

Model Problem IV.5. Find the Fourier series for a function which is square-integrable but has an infinite jump.

Solution. The example we look at is x-1/3. Like the previous example, the function is odd, and there will only be sine contributions. The first several coefficients are numerically: 


             1.710689953277251171
             0.3736036030974259086
             0.7368949814369842807 
             0.2752161613201822172 
             0.5010701734296802681 
             0.2252029789384463762
Here is a comparison of the function and its Fourier approximation:

(graph and approximation of x^(-1/3))

Think about the series calculation at x=0 and at x=+/-1. Did what you expected happen?

As we have seen in computations, symmetries can be used to simplify the computation of Fourier series. Recall that a function is even if f(-x)=f(x) for all x, and odd if f(-x) = -f(x) for all x. If f is a periodic even function, then all of the sine coefficients bn = 0: bn is the inner product of f with a sine function, but any even function and any odd function are orthogonal on the interval [-A, A]. (If f(x) is periodic with period L, then we can analyze it on any interval of length L, so even if it is originally defined on [0,L], we can just as well look at it on [-L/2, L/2]. Look at a graph of a periodic function to understand why.) For similar reasons, if f is periodic and odd, all the coefficients am = 0.

Any function can be written as the sum of an even function and an odd function, and the Fourier series picks out the two parts. If, for instance, we want to find the Fourier series of a function such as pi x - x2, - pi <= x <= pi , we can save some work by thinking about the symmetries.

Model Problem IV.6. Use symmetries to efficiently find the Fourier series for the function pi x - x2 on the interval -pi < x < pi .

Solution. This function is neither even nor odd when we change x into -x, but its "odd part" is pix, and its "even part" is -x2. The sine coefficients in the Fourier series will come entirely from the odd part and the even coefficients entirely from the even part; -x2 is orthogonal to all the Fourier cosine functions, and pix is orthogonal to all the Fourier sine functions. Moreover, the symmetry can be used to change the integration range so that it begins at 0:

The coefficients bn are calculated as follows:
    2 (2/(2 Pi)) Integrate[Pi x Sin[m x], {x,0,Pi}]
The product of two odd functions is even, so this integral is:
    -2 Pi (-1)^m /m
The final result was calculated by software. In a calculus class you would have probably calculated it by integration by parts. The calculation of the coefficients an can also be simplified. First, notice that only -x2 contributes:
    2 (2/(2 Pi)) Integrate[-x^2 Cos[m x], {x,0,Pi}]
The value of this integral is
    -4 (-1)^m /m^2

When you combine these, you have the Fourier coefficients for pix - x2 on this interval.

Whoever was the marketing expert consulted when they invented "complex numbers" should be fired! Complex numbers almost always make things much simpler, and this is true in particular for Fourier series. The key to the simplification is Euler's formula: ei x := exp(i x) := cos(x) + i sin(x),                           (4.1) It is useful to know several special cases, especially, exp(i n pi) = (-1)n.
and
exp(i (2k+1) pi/2) = (-1)k i, Euler's formula allows us to discuss complex numbers not only in the Cartesian system, z = x+iy, but also in the polar system: z = r exp(i theta), where theta, the argument of z, is the polar angle in the complex plans, and r = |z| (the modulus of z) is the distance from the origin. One of the truly great things about this formula is that is allows us to replace calculations with trigonometric functions with much easier calculations with exponential functions. As an example, let us think of how to calculate the Fourier series for the function f(x) = x, 0 <=x < L (repeated periodically, if you wish, outside this interval - notice that this is not a centered interval, so this calculation is different from the one done above.). a0 is the average of the function, so, obviously, a0 = L/2. For the other am, we have to calculate integrals of the form

    x * sin or cos(2 Pi x/L)

If you don't have software or integral tables handy, you can do these integrals with integration by parts. Or you can reason as follows:

    a_m = Re(integral of x by exp(2 pi i x/L)

Model Problem IV.7. Evaluate this integral by a method different from integration by parts.

Solution. We can use the following useful trick:

    (differentiate the integral of exp(A x) and then fix A=2 pi i/L           (4.2)

After the calculation, we set A=2 pi i m/L. Calling the integral in (4.2) INT, we find

    (continuing as above)

which starts looking nice after we substitute A=2 pi im/L, which reveals that these exponential terms are both equal to 1:

    (same answer as before)             (4.3)

This is purely imaginary, so all am = 0.

Model Problem IV.8. Evaluate the coefficients bn with a similar method.

Solution. STOP! Don't go back to the beginning of the calculation. We don't have to start all over, since

    the integral of x sin(a x) = Im(integral x exp(i a x))

which is the integral we already did. From (4.3) we see that this imaginary part is just

    b_m = -L^2 /2 Pi m^2

Actually, we can be much more systematic if we simply replace the complete set of functions (2.8) by (2.9):

    exp(2 \pi i n x/L)/Sqrt(L), - infinity < n < infinity

These are simply related to the sines and cosines by Euler's formula (4.1) and its easy consequence:

    cos(z) = (exp(iz)+exp(-iz))/2, sin(z) = (exp(iz)-exp(-iz))/2i

Definition IV.9. The Fourier exponential series is an expansion (for an arbitrary square-integrable function):

    f(x) = Sum(c_k exp(2 Pi k i x/L)                          (4.4)

Since the exponential functions are an orthonormal set, a familiar kind of calculation shows us that the formula for ck is:

    c_k = (1/L) integral of f(x)*exp(-2 Pi i k x/L)    (4.5)

Notice the tricky minus sign - this is a place where the complex conjugate in the inner product is important.

Let's observe how much more convenient this formula is than the one without complex numbers. There is only one sum, not two sums and a constant term. There is only one formula for ck, not separate ones for a0, am, and bn. The constant term c0 fits in with all the others, for instance, and is not put aside as a special case. The Parseval formula (3.4) is now much simpler:

    <f,g> = L Sum (c-bar_k * c-tilde_k                                     (4.6)

Here, the coefficients for g have the tildes. In particular,

    ||f||^2 = sum |c_k|^2.                                     (4.7)

Notice that the Parseval formula is similar to the Pythagorean theorem, since, other than the normalization factor L, it states that a certain length squared is equal to the sum of the squares of its components in an orthogonal basis. We see here, however, that the space in which f lives is infinite-dimensional.

The important thing to know is that the Fourier exponential series is completely equivalent to the usual "full" Fourier series. We will later look at the Fourier sine and Fourier cosine series, which are truly different series, but the exponential series is not. If we substitute with Euler's formula, the full series becomes the exponential series or vice versa. We can recombine ck with c-k as follows:

    c_k exp {2 Pi ikx/ L} + c_{-k}exp{-2 Pi ikx/ L} = (c_k + c_{-k}) cos{2 Pi kx/L} + (c_k - c_{-k}) i sin{2 Pi kx/L}

From this we get the connecting formulae that:

a0 = c0, am = cm+c-m, and bn = i (cn - c-n}).    (4.8)

Another useful fact to know about Fourier series is that you can normally safely differentiate them and integrate them, as long as the functions you are representing are periodic and differentiable. Actually, even if the function is not smooth, these manipulations are still often possible. This topic is explored in the next chapter.

Exercises.
  1. If f is a real-valued function, the coefficients am and bn are real-valued. What do you conclude from (4.6) about the coefficients ck in this case? What can you conclude about the ck if f is an even function? if f is an odd function? (solution by Mathematica)
  2. Calculate the Fourier series for the saw-tooth function you get as the periodic extension of f(x) := 1 - |x| from the basic interval [-1,1]. Differentiate the series term by term and compare with the Fourier series for the derivative of f(x). Notice that the Fourier series is not bothered by the corners in the function at -1,0, and 1. (solution by Mathematica)
  3. Calculate the Fourier series for the saw-tooth function you get as the periodic extension of f(x) := x from the basic interval [-1,1]. Differentiate the series term by term and compare with the Fourier series for the derivative of f(x). Explain any discrepancies you find.
  4. Calculate the second derivatives of the Fourier series of the last two problems. Does the resulting series converge? If so, to what?
  5. Calculate the full and exponential Fourier series for the power functions x2, x3, and x4 on a) the interval [- pi, pi ], and b) the interval [0,1]. (Solution)
  6. Calculate the exponential Fourier series for the power functions of Exercises III.

Link to

  • chapter V
  • chapter III

  • Table of Contents
  • Evans Harrell's home page
  • Jim Herod's home page