A sort of diary for Math 421:02, spring 2004
Part 4: Fourier series and some classical PDE's

BACK TO THE
SECOND HALF OF
LINEAR ALGEBRA
Fourier series and
some classical PDE's

DateTopics discussed
Thursday,
April 29
I will begin by asking part 1 of the QotD. The indicated heat distribution is given initially, with boundary conditions u(0,t)=6 and u(10,t)=2. Sketch the limit of u(x,t) as t-->infinity.

As we discussed this, I tried to indicate how heat would flow "intuitively". But of course this intuition would be guided by our past thought experiements and computations. Notice that the function -(4/10)t+6 is a steady state heat flow: that is, it has the desired boundary conditions and is only a function of t. Therefore the difference between this function and the u(x,t) in the QotD is a steady state heat flow for the 0 boundary conditions, which must be 0.

I also remarked that this initial temperature distribution could not be the result of the time evolution of an earlier temperature distribution. The principle reason for that is because u(x,0) is not differentiable at x=1/2. Temperature distributions which have evolved using the heat equation must be smooth.

Then I discussed temperature distributions in a rod with insulated ends (as well as insulated sides). It took some convincing but the appropriate boundary conditions are ux(0,t)=0 and ux(L,t)=0: u doesn't change in the spatial directions across the ends of the bar. I then asked the second part of the QotD: if the ends of the bar are also insulated, what is the limiting heat distribution? Here we want heat not to pile up (I tried to indicate that heat really dislikes being at a local max or min): we want it to distribute itself as well as possible ("entropy", disorganization, should be increasing at all times). Therefore the limiting steady state heat distribution in this case has temperature constant, equal to the average of the initial temperature. In this case, the limiting (average) temperature is 7, and the limiting temperature distribution is the constant 7.

I wanted to justify this limiting fact, and to analyze the general boundary/initial value problem for the heat equation:
ut= uxx; ux(0,t)=0 and ux(L,t)=0; u(x,0)=f(x) given.

So for the fourth time (I think) we analyzed this problem with separation of variables. Assume u(x,t)=X(x)T(t). Then, as before, T'(t)X(x)=X''(x)T(t), so that T'(t)/T(t)=X''(x)/X(x). Again this must be constant, as before. We call it, as before, -.

Start with X(x): we know X''(x)/X(x)=-. Again, X''(x)+X(x)=0. The boundary conditions for this equation are from ux(0,t)=0 and ux(L,t)=0. Therefore X'(0)=0 and X'(L)=0. If <0, solutions are linear combinations of cosh(sqrt()x) and sinh(sqrt(x). X'(0)=0 shows that there's nothing coming from the cosh part, and X'(L)=0 shows that the sinh part must vanish, also. Therefore must be positive.

Solutions are X(x)=C1sin(sqrt()x)+C2cos(sqrt()x). X'(0)=0 means that C1=0 (the derivative of sine is cosine). And, finally, X'(L)=0 means sin(sqrt()L)=0. Therefore sqrt()L is an integer (usually written, n) times Pi, and =[n2Pi2]/L.

What about T(t)? Since T'(t)/T(t)=- =-[n2Pi2]/L, we know that T(t) must must scalar multiples of e-([n2Pi2]/L)t. Notice that if n>0, this -->0 very rapidly, and therefore the only term which persists is the n=0 term.

You can check in the book (pp.921-2) to see that the general solution for the time evolution of a heat equation in an insulated bar with insulated ends is
u(x,t)=(c0/2+SUMn=1infinitycncos([nPi x]/L)e-([n2Pi2]/L)t
where cn=(2/L)0Lf(x)cos([nPi x]/L) dx for n>=0.
c0/2 is the average value of f(x) over the interval [0,L], which agrees with the intuition we may have developed giving this constant function as the steady state solution for this problem as t-->infinity.

Steady state heat flows in one dimension are solutions of ut=kuxx which are functions of x alone. That is, uxx=0. These functions are just Ax+B where A and B are constants. So these steady state heat flows whose graphs are straight lines have no bumps and no bending. These are the only functions with those qualities. In two dimensions and more, things get more complicated. There steady state heat flows must satisfy uxx+uyy=0. Such functions are called harmonic functions, there are infinitely many linearly independent harmonic functions, and their behavior can be complicated.

The room was had numerous gnats in this last class, which plagued bothstudents and instructor. The class ended slightly early.

Study for the final exam!

Tuesday,
April 27
Derivation of the Heat Equation
I tried to "derive" the one-dimensional heat equation, also known as the diffusion equation. You should realize that the heat equation is also used to model a very wide variety of phenomena, such as drugs in the bloodstream diffusing in the body or gossip (!) or spread of disease in epidemics or ... The text discusses the temperature of a thin homogeneous straight rod, whose long side is insulated (no heat flow through that allowed) but heat can escape (or enter) through the ends.

I mostly followed the discussion in this web page. Another derivation is given on this page. and yet another is here.

If you are someone who learns more easily with pictures, as I do, you may find the various animations of solutions of the initial value problem for the heat equation on this web page useful. Two questions which appear on a web page linked to that reference are interesting to keep in mind:

  1. Describe in words the time evolution of a one-dimensional heat flow for which the ends of the bar are both kept at zero degrees. How does your answer depend on the initial temperature distribution?
  2. Describe in words the time evolution of a one-dimensional heat flow for which the ends of the bar are kept at different temperatures. How does your answer depend on the initial temperature distribution?
Separating ...
I separated variables for the initial value problem u(x,0)=f(x) with boundary conditions u(0,t)=0 and u(Pi,t)=0. (The text does this with a more general L instead of Pi.) I assumed a solution to ut=uxx looked like u(x,t)=X(x)T(t), as in the text. The development is very similar to what was done for the wave equation on Tuesday, April 20. I got X''(x)/X(x)=T'(t)/T(t). Therefore each side must be constant. Look at the x side first. The boundary conditions imply that we need to solve the ODE problem X''(x)=(constant)X(x) with X(0)=0 and X(Pi)=0. Therefore to get a non-trivial solution, we must have constant=-(nPi)2 for some positive integer n, and then X(x)=sin(nx). Now here is where we differ from the wave equation: we need to solve T'(t)=-n2T(t). Solutions of this are multiplies of T(t)=e-n2t. Even for small n's this exponential function rapidly decreases from 1 to nearly 0. A general solution would be
u(x,t)=SUMn=1infinitycnsin(nx)e-n2t
where for the initial value problem u(x,0)=f(x) the cn's would be the coefficients of the Fourier sine series of f(x).

What should be expected about heat flow if the ends are kept at temperature 0? We discussed this, and then I gave out this which showed a good approximation to a solution. The heat does flow out, and the heat profile gets smoother. Notice this initial condition is one which I analyzed for the wave equation, and the solutions there were very different. Also notice that the original square profile can't be the result of a heat diffusion (otherwise it would be smooth!) so we can't run the heat equation backwards. The wave equation is time-reversible.

Comparing qualitative aspects of solutions of the wave and heat equations
Wave equationHeat equation
  • Conservation of energy: in what we studied, the sum of potential and kinetic energy is constant (other models allow friction, etc.)
  • Corners and shocks allowed to persist
  • Finite propagation speed: model has a top speed (speed of sound or light or ...)
  • Time can be reversed: waves just reverse direction, etc.
  • Diffusion: heat oozes away
  • Model is smoothing: even if initial temperature distribution is kinky, for later time, the temperature is smooth
  • Propagation speed turns out to be infinite, but far-off effects for small time are exponentially small.
  • Time can't generally be reversed: you can't make heat flow backwards.

I asked what would happen if instead of the boundary condition u(0,t)=0 and u(Pi,t)=0 we substituted u(0,t)=0 and u(Pi,t)=5. First I asked people to think about the physical situation. The general reply is that the temperature distribution would eventually tend to a steady state where one end (at 0) was 0 and the other end, at Pi, was 5. So what I did was write uss(x,t)=(5/Pi)x (this is a guess, but based on the physical suggestion!). It turns out that this uSS ("ss" for steady-state) does solve the heat equation: check that the t derivative is 0 and so is the second x derivative. If we take a solution of the boundary value with u(0,t)=0 and u(Pi,t)=5 (call this solution unew) then unew(x,t)-uss(x,t) does solve the heat equation (since the heat equation is linear) and this solution satisfies the old boundary conditions, u(0,t)=0 and u(Pi,t)=0.But we know because of the )e-n2t factors that this solution must -->0 as t-->infinity. So the difference between unew(x,t) and uss(x,t) must -->0 as t-->infinity. Therefore in fact the steady state is the limit for a solution of this boundary value problem.

I should have asked a QotD. It could have been this: suppose one end of a thin homogeneous bar is kept at temperature 7 and the other is kept at temperature -2. Suppose some initial temperature distribution is given (something really complicated). Graph the resulting solution of the heat equation for t=10,000.

Think about this: when should I have office hours and possible review sessions for the final? I'll try to do some problems from section 17.2 at the next (last) class.

Thursday,
April 22
The spring, a one-dimensional system
Everybody knows the vibrating spring, one of the simplest idealized mechanical systems. Here the dynamics are shaped by Hooke's law, that displacement from equilibrium results in a force of magnitude directly proportional to the absolute value of the displacement and directed opposite to the displacement. This together with Mr. Newton's F=mA gives us D''(t)=-(k/m)D(t) where D(t) is the displacement from equilibrium. The minus sign comes from the word "opposite". Solutions are sines and cosines, yielding the model's solution as vibrating up and down. (A typical vibrating string is the ammonia molecule, NH3.)

Looking at energy
This model conserves energy, where energy is here the sum of kinetic energy and potential energy. I won't worry about constants, since I can almost never get them right. The kinetic energy is (1/2)m(velocity)2. But velocity is D'(t). What about the potential energy? This is essentially the amount of work needed to move the spring from equilibrium to another position. Well, moving the spring against the varying resistant force of kD(t) means we need to integrate it with respect to D. This gets me (1/2)kD(t)2. The total energy in this model is therefore (1/2)m(D'(t))2+(1/2)k(D(t))2. I can try to check that energy is conserved by differentiating this total energy with respect to time. The result (using the chain rule correctly) is (1/2)m2(D'(t)D''(t))+(1/2)k2(D(t)D'(t)). This is D'(t)(mD''(t)+kD(t)). But the differential equation of motion says exactly that mD''(t)+kD(t) is always 0! So energy is conserved. Both D(t) and D'(t) are sines/cosines, and the sum of their squares is constant (sin2+cos2 is a well-known constant).

Taking a static picture at a fixed time can be rather deceptive. For example, if the spring is wiggling (no friction: an ideal spring) it always keeps wiggling. If we happened to snap a picture when the spring is at the equilibrium point it is still moving. The kinetic energy is largest at that position, and the potential energy is smallest. When the spring is highest or lowest, the potential energy is large, and the kinetic energy is 0 (the spring isn't moving, in the sense of calculus: D'(t)=0).

I hope that energy considerations will help your understanding of the vibrating string as well as the spring. In fact, we could almost understand the vibrating string as a bunch of springs tied together. This isn't totally correct, but if it helps ... In a very simple case (Warning! big file!) you can see the string in equilibrium position momentarily, but it has relatively high velocity (and therefore large kinetic energy) at that time.

The D'Alembert solution: the two waves
Here the vibrations of the string are written as a sum of left and right moving waves. Again, y(x,t) is supposed to be the height of the string at time t and position x. We write y(x,t)=K(x-ct)+J(x+ct). Mr. Marchitello helped me to understand that the K part is moving to the right and the J part is moving to the left. Let me use the chain rule:
yt(x,t)=K'(x-ct)(-c)+J'(x+ct)c
ytt(x,t)=K''(x-ct)(-c)2+J''(x+ct)c2
yx(x,t)=K'(x-ct)+J'(x+ct)
yxx(x,t)=K''(x-ct)+J''(x+ct)
Now we see that ytt(x,t)=c2yxx(x,t), so indeed this satisfies the wave equation.

Comparing initial conditions and the wavea
We worked hard to solve (via separation of variables and Fourier series) the initial value problem for the wave equation. For this initial value problem, we considered y(x,0)=f(x) and yt(x,0)=g(x) as given. (Again Mr. Marchitello commented that I seemed to be ignoring the boundary conditions, and I agreed that I was. These details, which are not so trivial, are discussed in the text on pp. 871-872.) So how can we go from the K and J description to the f and g description, and back, and does this help? It turns out that it will help a great deal with the understand of the physical model. Since y(x,t)=K(x-ct)+J(x+ct) we know that y(x,0)=K(x)+J(x)=f(x), and then since yt(x,t)=K'(x-ct)(-c)+J'(x+ct)c we know that yt(x,0)=K'(x)(-c)+J'(x)c=g(x). So this is how we could get f and g from K and J.

A short intermission
The other example we discussed last time had f(x) equal to 4 on [Pi/3,Pi/2] and 0 elsewhere on [0,Pi], while g(x) was 0 everywhere on [0,Pi]. I produced a picture (Warning! big file!) of a Fourier approximation to this string and its motion through time (for a short period of time). What the heck is going on? Well, we know that K(x)+J(x)=f(x) and (c=1) K'(x)-J'(x)=0, so that K'(x)=J'(x). Can you think of what K(x) and J(x) could be? You are allowed to guess. I can guess, especially when helped by that useful moving picture. K(x) and J(x) are both equal to a bump half the size of f(x): they are both equal to 2 on [Pi/3,Pi/2] and both 0 elsewhere. But one travels left and one travels right. We just happened at initial time to take a picture when the bumps reinforced. When the bumps reach the end of the interval, the boundary conditions (the string is pinned down) cause the bumps to be reflected negatively. Several engineering students seemed to think that this was correct physically. To me much of what is happening in this example is unexpected.

Back to the general case
We'll try now to get K and J from f and g in general. Let's integrate yt(x,0)=K'(x)(-c)+J'(x)c=g(x) from 0 to x. But first think: should we really write 0xK'(x) dx? There are too many x's. The x's inside the integral could easily get confused with the x on the upper limit of integration. So, following a suggestion of Ms. Sirak, I'll use w for the integration variable. Now we compute:
0xK'(w) dw=K(w)|w=0w=x=K(x)-K(0), and 0xJ'(w) dw=J(x)-J(0). So if I do the arithmetic correctly, the equation K'(x)(-c)+J'(x)c=g(x) leads to [K(x)-K(0)](-c)+[J(x)-J(0)]c=0xg(w) dw or (divide by c, and move J(0) and K(0)):
-K(x)+J(x)=(1/c)0xg(w) dw+J(0)-K(0).
This together with K(x)+J(x)=f(x) allows me to solve for K(x) and J(x) in terms of stuff involving f and g.

Here is what you get:
K(x)=-[1/(2c)]0xg(w) ds+(1/2)f(x)-(1/2)J(0)+(1/2)K(0)
and
J(x)=[1/(2c)]0xg(w) ds+(1/2)f(x)-(1/2)K(0)+(1/2)J(0)
I remarked in class and I will repeat now that this is one of the derivations in the whole course that I think is very useful to know about. It really does lead to very useful ideas in the "real world". Let's see: the minus sign in front of the integral for K(x) can be lost if we interchange the top and bottom limits:
K(x)=[1/(2c)]x0g(w) ds+(1/2)f(x)-(1/2)J(0)+(1/2)K(0)
and also I am really interested in y(x,t). But y(x,t) is K(x-ct)+J(x+ct). So in the formulas for K and J, make the substitutions x-->x-ct in K and x-->x+ct in J. Here is the result for y(x,t):
[1/(2c)]x-ct0g(w) ds+(1/2)f(x-ct)-(1/2)J(0)+(1/2)K(0)+ [1/(2c)]0x+ctg(w) ds+(1/2)f(x+ct)-(1/2)K(0)+(1/2)J(0)
Now look closely, and see things collapse. Indeed (!) the K(0) and J(0) just cancel. And the integrals? Well, the intervals are adjoining, so the integrals can actually be combined into one. Here is the resulting formula:
y(x,t)=[1/(2c)]x-ctx+ctg(w) ds+(1/2)f(x-ct)+(1/2)f(x+ct)
By the way, it is natural (to me) that g is integrated. g is a velocity, and y(x,t) is a position. So g's information has to be integrated to get the units to agree.

Using the formula to understand the physical situation
Look to the right. What's displayed is a different way to look at solutions of the wave equation. The function y(x,t) is a function of two variables, one which locates a position on the string, x, and one which indicates time, t. Let's look down on (x,t) space.

I would like to know about y(x,t) and what part of the initial data influence y(x,t). Certainly f's values at x+/-ct influence y(x,t) (look at the formula above!) and g's values in the interval [x-ct,x+ct] influence y(x,t). Mr. Tahir bravely came to the board and completed a diagram which looks much like what is draw on the left. The dotted lines have slopes +/-(1/c). (I got this wrong in class!) So the influence of the initial conditions is as shown.

Here is yet another way to show what's going on. I could imagine a short interval on the time=0 line, say from x to x+(delta)x (as usual, (delta)x is supposed to be small). I could imagine changing f or g or both on that interval, wiggling these initial conditions. Then what parts of the string could be affected at a later time, t=t0? Well, if you think about it and look at the formula, the only parts could be in the upward part of the trapezoid sketched. The wave equation models phenomena which have a finite propagation speed (here the speed is c). Such phenomena include sound and light. The slanted edges of the trapezoid sketched are usually called the light cone and it represents the spreading the maximum speed that the model allows. This name reflects the physical theory that light travels at the maximum universal speed, and any other "signal" travels more slowly. I remark that Mr. Woodrow courageously drew a picture similar to this on the board.

The QotD was ... well ... student evaluations were requested. The local supply of explectives was exhausted. I remarked that the welfare of my family members and pets (alligator, boa constrictor, etc.) depended on these reviews. Ms. Lu kindly agreed to take charge of the evaluations. I wanted to bribe her to fill in the unused ones in a way beneficial to me. This was not done.

16.2: #3 (try n=20 and t=0, t=.5, t=1), #16 a and c, 16.4: #1, #22 (try t=0 and 1 and 2). Please read the text. I remarked that I wanted to spend next week discussing the heat equation, and that I would therefore not cover section 16.7: we don't have enough time. Technically 16.7 covered the wave equation in a rectangle. What's used are several variable Fourier series, but most of the qualitative aspects (about speed of influence, etc.) are much the same as in one variable.

Tuesday,
April 20
Copying the book
The details of the material in section 16.2 are enormous. To me this material is probably the most daunting (dictionary: discouraging, intimidating) of the whole semester. What I wrote during the beginning of this lecture was essentially copied from the text. Although I have my own opinions about how the material should be organized and presented, my belief is that right now is the wrong time to display the results of these opinions. So here we go:

The setup
The text starts with a function y(x,t) which is supposed to describe the motion of a string at position x (x between 0 and L) and time t (t non-negative). The string is fastened at the ends, and also has an initial position and velocity. Physical intuition (?!) is supposed to tell you that later positions of the string are determined by these initial conditions and boundary conditions, if y(x,t) obeys the wave equation, ytt=c2yxx.

Boundary conditions
Here y(0,t)=0 and y(L,t)=0 for all t: the string is fastened down at the ends.

Initial conditions
y(x,0) is some function f(x) and yt(x,0) is some function g(x). f(x) is supposed to be initial position and g(x) is supposed to be initial velocity.

The importance of linearity
We will first solve the wave equation with initial data any f(x) and the velocity equal to 0. Then we will solve it with velocity g(x) and initial position equal to 0. The general solution will be the sum of these two solutions. Also important in all this stuff (this is very involved) is that the boundary conditions are homogeneous, so the sum of solutions satisfying the boundary conditions still satisfies the boundary conditions.

In what follows I first will ask that y(x,0)=f(x) and yt(x,0)=0 for all x.

The boundary conditions lead to eigenvalues via separation of variables
The word eigenvalue is also used with differential equations. Please read on. We make the completely wild assumption that the solution y(x,t) can be written as X(x)T(t). This assumption is only being taught to you here because it is useful when considering this equation and a whole bunch of other equations. So learn it!
Then (pp.862-863 of the text): XT''=c2TX''. This becomes X''(x)/X(x)=c2T''(t)/T(t). This is sort of a ridiculous equation. The left-hand side is a function of x,and the right-hand side, a function of t. Only constants are such functions. The constant is called -. The reason for the minus sign is to make other computations easier.

Since X''(x)/X(x)=-, we get X''(x)+X(x)=0. The boundary condition y(0,t)=0 means that X(0)T(t)=0 for all t. Well, if the T(t) isn't supposed to be 0 always, we must conclude X(0)=0. Similarly the boundary y(L,t)=0 yields X(L)=0. So we are led to the problem:
X''(x)+X(x)=0 with X(0)=0 and X(L)=0.
What can we say about the boundary value problem for this ordinary differential equation? Well, I know a solution for it. The solution X(x)=0 for all x is a solution. This resembles what happens in linear homogeneous systems. There's always the trivial solution. Are there other solutions? Consider the following.

Stupid example
Suppose we are looking for a solution to X''(x)+56X(x)=0 and X(0)=0 and X(48)=0. Well, solutions to X''(x)+56X(x)=0 are C1sin(sqrt(56)x)+C2cos(sqrt(56)x). If X(0)=0, then since cos(0)=1 and sin(0)=0, this means C2=0. So the solution is now C1sin(sqrt(56)x). But if X(48)=0, we plug in x=48: C1sin(48sqrt(56))=0. When is sin(x)=0? It is only when x is an integer multiple of Pi. But I don't think that 48sqrt(56) is an integer multiple of Pi. (Check me, darn it: compute [48sqrt(56)]/Pi and see if this is an integer!) This means that if the equation C1sin(48sqrt(56))=0 is true, then C1=0. Therefore (blow horns and pound drums!!!!) the boundary value problem X''(x)+56X(x)=0 and X(0)=0 and X(48)=0 has only the trivial solution.

O.k., back to the general case: X''(x)+X(x)=0 and X(0)=0 and X(L)=0. The solutions are C1sin(sqrt()x)+C2cos(sqrt()x). The X(0)=0 condition gives C2=0. What about inserting x=0 in the other piece? C1sin(sqrt()L)=0 has a solution with C1 not equal to 0 exactly when sqrt()L=nPi for an integer n. Then is [n2Pi2]/L2. Selecting as one of these values (with n=1 and n=2 and n=3 and ...) allows the possibility of non-trivial solutions. If is not only of these values then the only solutions are the trivial solutions. These special values of are called the eigenvalues of this boundary value problem.

Now let's consider the T(t) function. Well, (1/c2)T''(t)/T(t)=- (where is now identified as above!). So T''(t)+[n2Pi2c2]/[L2]T(t)=0. We haven't used the initial condition yt(x,0)=0 for all x. This leads to T'(0)X(x)=0 for all x. Again, either we have a really silly (and trivial) solution X(x)=0 always or we require that T'(0)=0. So if we know that T''(t)+[n2Pi2c2]/[L2T(t)=0 then T(t)=C1sin([nPi c/L]t)+C2cos([nPi c/L]t). What's T'(t)? It has a sine and a cosine, and requiring that T'(0)=0 menas that the coefficient in front of cosine (in the derivative!) must be 0. Then C1 must be 0. So T(t) is C2cos([nPi c/L]t).
We then get the solution y(x,t)=X(x)T(t) is (CONSTANT)sin([nPi/L]x)cos([nPi c/L]t). Here again n is a positive integer. We have done lots of work to get to this place.

Linearity of the wave equation
When t=0 in the solution listed above we get (CONSTANT)sin([nPi/L]x). So if you (we?) happen to be lucky, and your c=1 and your L=Pi and your initial position is 36sin(5x) then you can just "guess" a corresponding solution to the wave equation: 36sin(5x)cos(5t). Below to the left is the initial shape of this wave form. Then to the right is an animation of what happens as t goes from 0 to 4.

The pictures are linked here.
Warning!
Large file (~1 meg) to be transmitted.

Remember this is for an idealized vibrating string with no friction, etc., so it just keeps going back and forth endlessly (remember Hooke's law and its validity: the ideal spring compared to a spring vibrating in thick honey, etc.) Energy is constant in this system. When the string at max amplitutde, it is not moving: all the energy is potential. When the string overlays to x-axis, its energy is all kinetic. This is indeed just like an idealized vibrating spring.

But what happens if we want a more complicated initial wave? Well, we take advantage of linearity. We can add up the simple waves, and we get y(x,t)=SUMn=1infinitycnsin([nPi/L]x)cos(nPi c/L]t) (this is on p.864 of the text). When t=0, this becomes y(x,0)=SUMn=1infinitycnsin([nPi/L]x) since cosine of 0 is 1. Let me do an example.

Here I will suppose that c=1 and L=Pi. Then y(x,0)=SUMn=1infinitycnsin(nx). Let me ask that the initial shape of the string be like this: 0 if x<Pi/3 or if x>Pi/2, and 4 for x in the interval [Pi/3,Pi/2]. This is a rectangular bump, and please ignore the fact that it is not continuous. I would like to represent this function as a sum of sines. We have the technology: we can extend the function as an odd function and compute the Fourier coefficients. Here is the result graphically. The picture on the left shows the original rectangular bump, and compares it with the sum of the first 40 terms of the Fourier sine series (obtained by extending the function in an odd way to [-Pi,Pi]). This static picture shows the usual phenomena, including the Gibbs over- and undershoot. The next picture shows the resulting vibration, if you can follow it, during the interval from t=0 to t=4. Notice that fixing the boundary points at 0 "reflects" the waves in a reversed fashion off the ends. This is quite fascinating to me.

The pictures are linked here.
Warning!
Large file (~1 meg) to be transmitted.

The initial value problem (velocity)
Here I suppose that y(x,0)=0 for all x, and yt(x,0)=g(x), for x in [0,L]. Almost the same computations are used. Look at p. 868. In this case, the boundary value problem for the t variable makes the cosine term vanish, because y(x,0)=X(x)T(0)=0. So here y(x,t)=SUMn=1infinitycnsin([nPi/L]x)sin(nPi c/L]t) and you find the cn's from the derivative:
yt(x,0)=SUMn=1infinitycn[nPi/L]sin([nPi/L]x)=g(x)
where the coefficients are from the sine series of g(x). Wow. I don't want to make more pictures because the storage space needed for the animated gifs is enormous.

So what can you do ...
Not much is easy to do by hand. I would sincerely hope that anytime you need to work with these things, you will have a program like Maple or Matlab or Mathematica at hand.

Waves, moving right and left
I think it was Mr. Sosa who kindly let me use the string from the hood of his jacket to show a rather different representation of the wave shapes. I wiggled one end, and Mr. Sosa wiggled the other. The result was something like what I show below.

Here is an example of a solution to the wave equation which I checked in class:
y(x,t)=(x-t)342+(x+t)221
To check this you must remember the Chain Rule carefully. The precise numbers don't matter. It turns out that the dependency on x+t and x-t are the only things that matter. As t grows, the x-t part represents a wave moving to the right and the x+t part represents a wave moving to the left. This is called the D'Alembert representation of the solution, and I will try to understand it more on Thursday and connect it with the work done earlier. Please look at section 16.4 of the text.

The QotD was to look at two bumps, one moving right and one moving left, and draw a picture some time later of what the profile would look like. It was not a wonderful QotD. Please read the text!

Thursday,
April 15
Homework: DONE
13.2 #9, was presented on the board (thanks to Mr. Araujo), 13.3 #5, was presented on the board (thanks to Ms. Mudrak), and 13.5 #2, was presented on the board (thanks to Mr. McNair). These problems all had their technical difficulties (especially 13.3, #5, when done "by hand") but are all straightforward applications of formulas in the tex.

PDE's
In 244 you learned about ordinary differential equations and initial value problems. The basic theory can be stated easily: a nice ODE with the correct initial conditions (order n equation with n initial conditions) has a solution. Of course, "solving" an ODE, in the sense of getting a neat symbolic solution or even approximating the solution well, required tricks and lots of thought. But the theory is there.

By contrast, partial differential equations are horrible. The equations we will analyze in the short time remaining in this course (the classical wave and heat equations) have been studied a great deal for two hundred years, but even these PDE's have properties which are difficult to understand, both theoretically and computationally. More general PDE's have very little theory supporting their study. This is a pity, because many "real" phenomena are best modeled by PDE's.

The Wave Equation
The wave equation is an attempt give a mathematical model of vibrations of a string (space dimension=1) and a thin plate (space dimension=2). Please read the assigned sections in the text in chapter 16. It will be very helpful if you read these sections in advance, since the computational details are intricate!

In one dimension, the wave equation described the dynamics of a string imagined drawn over the x-axis, whose height at position x at time t is y(x,t), and gets a condition which is satisfied by y(x,t): ytt=c2yxx. Here the subscripts tt and xx refer to partial derivatives with respect to time (t) and position (x). To me, a really neat part of the subject is discovering why this equation is true, based on simple physical assumptions. Here and here are two derivations of the wave equation I found on the web. Please note that the wave equation describes small vibrations of an idealized homogeneous string. If your string is lumpy (?) or if you want to pull it a large amount, then such a simple description is not valid. The wave equation will model fairly well the motion of a guitar string, say.

One nice thing to notice is that ytt=c2yxx is linear: if y1 and y2 are solutions, then the sum y1+y2 is a solution, and so is ky1 for any constant k. So the solutions of the wave equation form a subspace. This is frequently called the principle of superposition (yeah, there are over 86,000 web pages with this phrase!). A random function is not likely to solve the wave equation, which essentially asks that yxx and ytt be constantly proportional: that's quite a stiff requirement. Here are some candidates and how I might check them:

y(x,t)ytytt yxyxxConstantly
proportional?
x5+t3 3t26t 5x420x3 No!
t3ex 3t2ex6tex t3ext3ex No!
sin(xt) x cos(xt)-x2sin(xt) t cos(xt)-t2sin(xt) No!
cos(x2+5t) -5sin(x2+5t)-25cos(x2+5t) -2x sin(x2+5t)-2sin(x2+5t)-4x2cos(x2+5t) No!
sin(3x)cos(5t) -5sin(3x)sin(5t)-25sin(3x)cos(5t) 3cos(3x)sin(5t)-9sin(3x)cos(5t)Yes! c2=9/25
e4tcos(5x) 4e4tcos(5x)e4tcos(5x) -5e4tsin(5x)-25e4tcos(5x) No!
 The constant of proportionality, c2, must be positive!
log(6x+7y) 7/(6x+7y)-49/(6x+7y)2 6/(6x+7y)-36/(6x+7y)2 Yes! c2=49/36

My reasons for compiling this table are to warm up by computing some partial derivatives in case you have forgotten how, and, more seriously, to show you how rarely the equation is satisfied.

Setting up a boundary value problem
The physical aspect of the wave equation is important in what we should expect of solutions. It turns out that there are lots of solutions, but which ones will make physical sense? Here maybe it does help to think of a guitar string. Again, I'll follow the text closely (please read the text!). Usually the guitar string is fastened at the ends, so right now we will assume y(0,t)=0 and y(L,t)=0. But we need to know more information: we need the position of the string at time t=0, y(x,0) for x in [0,L], and, further, we need to know the velocity of the string at time t=0, yt(x,0) for x in [0,L]. To the right is my effort to show you a picture of what is happening. There is an x-axis on which the string is "vibrating". The other horizontal axis is t, time. Time is progressing futurewards (?) as t increases. The height from the (x,t) plane is y(x,t), the string's height at time t and position x. Then y(x,t) actually describes a surface, and each slice with a fixed time, t, is a frozen picture of where the string is at time t. I specify that the string is fixed at height 0 (blue line and dashed line) at x=0 and x=L. But I also specify the initial position when t=0 (that's the green curve) and the initial velocity when t=0 (that's the red curve). The red curve declares that the string is moving up when t=0, and moving up faster when x is small than when x is large. I also tried to show a picture of the string for a little bit later time, and, indeed, wanted to make clear that the part of the string for x small has moved up more than the part for x large.

This is a complicated and maybe silly picture. You can look at it, but maybe you should also try creating your own pictures. Your own pictures will likely make more sense to you. The boundary conditions y(0,t) and y(L,t) and the initial conditions y(x,0) and yt(x,0) for x in [0,L] are generally what people expect with the wave equation. Recognizing that these conditions are appropriate comes from the physical situation.

Trying to solve the wave equation
There are several standard techniques for solving the wave equation. Here is one of them, called separation of variables. Look at ytt=c2yxx and assume that y(x,t)=X(x)T(t): that is, y(x,t) is a product a function of x alone and a function of t alone. Then ytt is X(x)T''(t) (you may need to think about this a bit!) and yxx is X''(x)T(t). If this weird function satisfies the wave equation, then X(x)T''(t)=c2X''(x)T(t) so that T''(t)/T(t)=c2X''(x)/X(x). On the left side of this equation we have a function of t alone, and on the right side of the equation, we have a function of x alone. The only function that is a function of both x alone and t alone is a constant. This constant is usually called - (sorry about the minus sign, but everyone uses it, and you will see why, very soon). Therefore, if X(x)T(t) solves the wave equation, I know that T''(t)/T(t)=- so that T''(t)+T(t)=0 and I know a similar statement for X(x). So now I'll use what I know about ODE's to find eligible T(t) and X(x). Please look at the textbook. I will continue this next time.

A proposal for the Final Exam
I suggested that instead of having a comprehensive final covering all parts of the course that I do the following:

  1. Write an exam, whose length and difficulty will be appropriate (similar to the past two exams) covering the Fourier series part of the course. Everyone would take this, and I would grade it, and count it as one of three similar exams in the course.
  2. Optionally offer a few problems on linear algebra and a few problems on Laplace transforms for students who wish to improve their grades in the exams.

After some discussion and e-mail with students, I would like to specify the possible improvement so that no one can get more than 100. Here is one strategy, which is historically and physically interesting: if a student got GLT on the Laplace transform grade, and then got PLT on the optional few Laplace transform problems then (GLT is in the interval [0,100) and PLT, in the interval [0,{20 or 25 or 30}], I would replace (for term grade purposes) GLT by (GLT+PLT)/(1+GLT·PLT/104). Some of you may recognize this as Lorentz addition of grades! So the Lorentz sum of 55 and 15 is 64.67 (I'll use approximations here). Using this formula, students could never score more than 100.

Comment and information
Here is a quote from a web page entitled How Do You Add Velocities in Special Relativity?:

w = (u + v)/(1 + uv/c2)
If u and v are both small compared to the speed of light c, then the answer is approximately the same as the non-relativistic theory. In the limit where u is equal to c (because C is a massless particle moving to the left at the speed of light), the sum gives c. This confirms that anything going at the speed of light does so in all reference frames.

I think this penalizes students who got a fairly good grade already too much, though. For example, the Lorentz sum of 75 and 15 is only 80.90. So maybe I should find something better.

There was no QotD! But I did ask people to: Begin reading the sections of chapter 16. Please hand in 13.5, #3 and #4, and 16.1, #4.

Tuesday,
April 13
I had prepared a sequence of handouts, which are reproduced here as gifs. I basically wanted to discuss a sequence of examples to illustrate
  • The definition and computation of Fourier coefficients.
  • The behavior of the partial sums of Fourier series.
  • The (theoretical?) limiting behavior of the partial sums of Fourier series: what is the sum of the whole Fourier series?
  • The qualitative behavior of the Fourier series.
  • The Parseval equality.
  • Odd and even extensions of functions and the Fourier sine and cosine series which result.
This is a lot to think about but we can only complete this by beginning, so let's begin.

I remarked that Taylor series are good for "local" computations near a fixed point, in some "tiny" interval around the point. Fourier series are a wonderful thing to use on a whole interval. I recalled some general facts: the definitions of Fourier coefficients and Fourier series and how the partial sums and the whole sum of the Fourier series behaves compared to f(x). I also supplied a (very short) table of integrals we would need today:

FunctionAntiderivative
cos(Kx)(1/K)sin(Kx)+C
sin(Kx)-(1/K)cos(Kx)+C
x cos(Kx)(1/K)x sin(Kx)+(1/K2)cos(Kx)+C
x sin(Kx)-(1/K)x cos(Kx)+(1/K2)sin(Kx)+C

A first example (an old QotD)
My first example analyzed the QotD from the previous regular class meeting. So f(x) is 0 for x<0 and is 1 for x>=0. In everything I am doing today, I will only care about the interval [-Pi,Pi]. If you run across a situation where you "care" about the interval [-5,11], say, probably you should first take x and change it into y=(2Pi/16)x-(3Pi/8). Then as x runs from -5 to 11, y would go from -Pi to Pi. Anyway, I used the following Maple commands to create the pictures I will show and also the resulting computations.

Maple commands
Here is how to write f(x), a piecewise defined function:

f:=x->piecewise(x>0,1,0);
This command creates a function q(n) which produces the cosine Fourier coefficients:
q:=n->(1/Pi)*int(f(x)*cos(n*x),x=-Pi..Pi);
and here is r(n) which produces the sine Fourier coefficients:
r:=n->(1/Pi)*int(f(x)*sin(n*x),x=-Pi..Pi);
Here I have a function Q(k) which creates a pair of functions. The first element of the pair is f(x), and the second element of the pair is the partial sum of the Fourier series up to the terms sin(kx) and cos(kx):
Q:=k->{f(x),(1/2)*q(0)+sum(q(n)*cos(n*x)+r(n)*sin(n*x),n=1..k)};
And then this command produces a graph of the previous pair.
QQ:=n->plot(Q(n),x=-Pi..Pi,scaling=constrained);
I try to write short and simple Maple commands rather than write big programs. I can't understand big programs easily. So all four of the pictures below were produced by the one command line (one: only one!)
QQ(0);QQ(2);QQ(5);QQ(10);
f(x) and the partial sums of the Fourier series for f(x) up to
n=0n=2n=5n=10
The hard part of this is figuring out what we can learn from such pictures. I mention that these pictures and the others I'll show below would be really irritating to obtain (at least for me!) without technological help.

Observation 1
Please realize that the partial sums (which are linear combinations of sin(mx) and cos(nx)) will always have the same values at -Pi and Pi since sin(mx) and cos(nx) are the same at -Pi and Pi. So these functions really "wrap around" on the [-Pi,Pi] interval. Look at the pictures!

Observation 2
The first picture tells me that (1/2)a0, that horizontal line, which here has height 1/2, is the mean square best constant approximating our f(x). That is, if I sort of take values at random, then the horizontal line will have the least error, on average, from the function's value. I hope you can see this here, since the function involved is quite simple.

Observation 3
Go away from the jumps in the function. That is, imagine you restrict your attention to x between -Pi and 0 or x between 0 and Pi. What's happening there? Here's a picture of a piece of the graph of f(x) and the partial sum, with n=20, in the interval for x between 1 and 2, away from the discontinuities.
You should notice that the error between the function and the partial sum is very very small. In addition, the partial sum, in this restricted interval, really does not wiggle very much. That is, its derivative is quite close to 0. The derivative of f(x) is 0 oin this interval, of course, since it is a horizontal line segment.

Observation 4
Here by contrast is a picture of what happens to the 40th partial sum (I'm very glad I'm not doing these by hand!) and f(x) on the interval [Pi-.5,Pi]. Of course (observation 0) the partial sum hops "down" towards 1/2 at x=Pi. But it is sort of interesting (as Mr. Dupersoy noticed!) that there is sort of an overshoot. The overshoot is very very near the jump. I will come back and discuss this a bit later. The overshoot (which has an accompanying "undershoot" for this f(x) at -Pi) is called Gibbs phenomenon, and is rather subtle.

The sum of the whole Fourier series
Theory tells me what the sum of the whole Fourier series should be. This function f(x) is "piecewise smooth". Theory then says that the sum of the whole Fourier series is f(x) if x is a point where f(x) is continuous. If the left- and right-hand limits of f(x) exist at a point and if these limits are different, then the Fourier series attempts to compromise (!) by converging to the average of the two values: graphically, a point halfway in between the limiting values of f(x). In our specific case, there are three (well, two, if you consider -Pi and Pi to be the "same" which they are for sine and cosine) points of discontinuity: -Pi, 0, and Pi. At each of those, the whole Fourier series of f(x) converges to the values shown.

Maple code for Parseval
This command

s:=k->(1/2)*q(0)^2+sum(q(n)^2+r(n)^2,n=1..k);
adds up the Fourier coefficients (properly squared) and is a partial sum of the infinite series appearing on one side of the Parseval formula.
pars:=k->[s(k),evalf(s(k)),(1/Pi)*int(f(x)^2,x=-Pi..Pi),evalf((1/Pi)*int(f(x)^2,x=-Pi..Pi))];
asks Maple to produce four quantities: the first one is the partial sum of the Parseval series symbolically evaluated (if possible), the second is a numerical approximation to that sum, the third is the integral symbolically evaluated (if possible), and the fourth is a numerical approximation to that integral. The command pars(10) gives the response
                        469876
                [1/2 + ---------, 0.9798023932, 1, 1.]
                               2
                       99225 Pi
and pars(20) gives the response
                   204698374253288
           [1/2 + ------------------, 0.9898762955, 1, 1.]
                                   2
                  42337793743245 Pi
Observation 5 (and something sneaky!)
Well, at least .989... is closer to 1 than .979..., so maybe the partial sums of the Parseval series are equal to one. Well, to tell you the truth, very few people in the real world use the Parseval series to compute the integral of f(x)2 on the interval [-Pi,Pi]. Instead, things like the following are done.

What's an? It is (1/Pi)-PiPif(x) cos(nx) dx, which is (for our specific f(x)!) (1/Pi)0Picos(nx) dx=(1/(nPi))sin(nx)|0Pi. But sine at integer multiples of Pi is 0, so an is 0. Errrr .... except for n=0, when that formula is not valid. We know that a0=1 (compute the integral of f(x)2, which is 0 on the left and 1 on the right!).

What's bn? It is (1/Pi)-PiPif(x) sin(nx) dx, which is (for our specific f(x)!) (1/Pi)0Pisin(nx) dx=-(1/(nPi))cos(nx)|0Pi. Now cos(n0) is cos(0) which is 1. What about cos(nPi)? It is either +1 or -1, depending on the parity (even/odd) of n. In fact, it is (-1)n. Therefore we have a formula: bn=(1/(nPi))((-1)n+1).

(I screwed up this calculation in class!)
What the heck does the Fourier series of f(x) look like? If you try the first few terms (all of the cosine terms are 0 except for the constant term!) you will get:
(1/2)+(2/Pi)sin(x)+(2/(3Pi))sin(3*x)+(2/(5Pi))sin(5*x)+(2/(7Pi))sin(7*x)+(2/(9Pi))sin(9*x)+...
Actually, what I did was just type the command sum(r(n)*sin(n*x),n=1..10); into Maple, copy the result and add on (1/2) in front. This allowed me to check all of my computations.

Then the infinite series in Parseval becomes
(1/2)+(4/Pi2)SUMn=0infinity1/(2n+1)2
where I factored out the 4/Pi2 from each term (uhhhh: 4/Pi2 is the square of 2/Pi). So we actually have (now put in (1/Pi)-PiPif(x)2 dx, which in our case is just 1):
1=(1/2)+(4/Pi2)(the big SUM).
We can solve for the sum. And thus (non-obvious fact coming!) The sum of 1/(all the odd squares) is Pi2/8.

Maple evidence for this silly assertion
I tried the following:

evalf(sum(1/(2*n+1)^2,n=0..100));
                             1.231225323
evalf(Pi^2/8);
                             1.233700550
which compares the sum of the first 100 terms with Pi2/8. And then, heck, what's time to a computer?, I tried the first 1000 terms:
evalf(sum(1/(2*n+1)^2,n=0..1000));
                             1.233450800
But the convincing command is this:
sum(1/(2*n+1)^2,n=0..infinity);
                                   2
                                 Pi
                                 ---
                                  8
because Maple can recognize certain infinite series and "remember" their values!

The second example
Here f(x) is the function 3x for x>0, and f(x)=0 for x<0. Here are some pictures:
f(x) and the partial sums of the Fourier series for f(x) up to
n=0n=2n=5n=10
I hope you can see some of the same qualitative features in these pictures as in the previous collection. So you should see in the first picture the best constant that approximates this f(x). In the succeeding pictures, you should see better and better approximations, away from the points of discontinuity (Pi, -Pi, which the trig functions regard as the "same"). Even the derivatives of the partial sums (for x not 0 or +/-Pi) are approaching the derivatives of f(x). There's also a Gibbs phenomenon, the bumps at +/-Pi. Here is:
         A graph of f(x)                       A graph of the sum of
                       the Fourier series for f(x)

Some Fourier coefficients for this f(x)
I tried to compute a11. This is (1/Pi)-PiPif(x) cos(11x) dx=(1/Pi)0Pi3x cos(11x) dx (because f(x) is 0 for x<0). I can antidifferentiate this (especially with the table above!). So it is (3/(11))x sin(11x)+(3/(11)2)cos(11x)|0Pi. I know that sin(11Pi) and sin(0) are 0. Also cos(0)=1 and cos(11Pi)=-1, so after computing, -6/(121Pi).

How about b10? This is (1/Pi)-PiPif(x) sin(10x) dx=(1/Pi)0Pi3x sin(10x) dx and this is -(3/10)x cos(10x)+(3/(10)2)sin(10x))|0Pi, and, evaluating carefully, this will be -3/(10).

Hey: I just wanted to show you that some of each kind of Fourier coefficient is not 0.

Cheap cosines?
We have seen that 3x, for x between 0 and Pi, is some crazy infinite linear combination of sines and cosines. What if evaluating cosines was very very cheap and easy, and evaluating sines was really expensive. I will show you a clever way to get 3x, for x between 0 and Pi, as a sum only of cosines.

The even extension
Fold over 3x on [0,Pi]. That is, define f(x) to be 3x for x>0 and f(x) is -3x for x<0. Then the resulting new f(x) is even, so that multiplying it by sin(mx), an odd function, and integrating the product from -Pi to Pi must result in 0. The Fourier series of this new function has only cosine terms.
f(x) and the partial sums of the Fourier series for f(x) up to
n=0n=2n=5n=10
We're using only cosines, and the two graphs look really close for n=10. The function and the sum of the whole Fourier series are identical in this case (the function is continuous in [-Pi,Pi] and its values at +/-Pi agree).

The odd extension
Change the game, so that now computing sines is cheap, and cosines are impossibly expensive. So we take 3x on [0,Pi] and extend it to an f(x) on [-Pi,Pi] with f(-x)=f(x). This has a simple algebraic definition (no longer piecewise): f(x)=3x for x in [-Pi,Pi]. Nows the integrals defining the cosine coefficients are 0 (cosine is even, this f(x) odd, so the product is odd). Here are the pictures which result:
f(x) and the partial sums of the Fourier series for f(x) up to
n=0n=2n=5n=10
Again, I hope you see some neat qualitative aspects.
         A graph of f(x)                       A graph of the sum of
                       the Fourier series for f(x)
These graphs differ at the endpoints.

What are these darn things?
I asked Maple to compute the beginnings of each Fourier series and I display the results. Please note (it is almost unbelievable!) that each of these series has the same sum, 3x, for x in the interval [0,Pi).
PictureBeginning of the Fourier series

Original f(x)
3 Pi   5 cos(x)              3           2 cos(3x)          
---- - -------- + 3 sin(x) - - sin(2x) - -------- + sin(3x) 
  4       Pi                 2             3 Pi             

Even extension from [0,Pi]
3Pi   12 cos(x)   4 cos(3x)   12 cos(5x)
--- - --------- - --------- - ----------
 2      Pi         3 Pi        25 Pi

Odd extension from [0,Pi]
                                   3            6
6 sin(x) - 3 sin(2x) + 2 sin(3x) - - sin(4 x) + - sin(5 x)
                                   2            5

The Maple command I used here was (1/2)*q(0)+sum(q(n)*cos(n*x)+r(n)*sin(n*x),n=1..5); separately for each function.

The Gibbs phenomenon: (maybe) a simpler version
Here is a simple model situation to help you understand the Gibbs phenomenon. I don't know any particular "real" application that needs to pay attention to this overshoot idea, but lots of people (about 44,500 web pages mention it, according to Google!) are curious about it.

In 1863, J. Willard Gibbs received the first U.S. doctorate in engineering. He saw that at a jump discontinuity, there is always an overshoot of about 9% in the Fourier series. The limit of the overshoot's behavior vanishes! This is a bit difficult to understand at first. Here are some pictures of a similar situation which may help you.


This is a picture of a moving bump. The bump is always non-negative, and is always 1 unit high. It gets narrower and narrower as n increases, and moves toward, but never actually reaches, 0. The point pn is at 1/n, and qn is at 2/n, and rn is at 3/n. The graph is made of four line segments. Notice that if you "stand" at 0, you never feel the bump. Of course, you also never will feel the bump if you are at any negative x. But, and this is most subtle, if x is positive, then for n large enough, the bump eventually passes you, goes toward 0, and your value goes to 0 and stays at 0. So the limit function of this sequence of moving bumps is 0 at every point, but each bump has height 1! This is similar to the bumps found by Gibbs. The partial sums display a sort of moving bump effect, but the limit function (the sum of the whole series) doesn't show any additional bump.


The QotD was: here is a function given graphically on the interval [0,Pi]. The answer(s) will be a series of four graphs. First, graph the odd extension of this function to [-Pi,Pi]. Then draw the graph of the sum of the whole Fourier series of the odd extension on the interval [-Pi,Pi]. For the third graph, draw the even extension of this function to [-Pi,Pi]. Then draw the graph of the sum of the whole Fourier series of the even extension on the interval [-Pi,Pi].

Thursday,
April 8
We had fun with the linear algebra exam, and only about a quarter of the students stayed past 10 PM.
Tuesday,
April 6
Real-life linear algebra [ADVT]
One step in the algorithm used by Google to develop PageRank (a way of deciding which pages are more reputable or useful) is finding the largest eigenvalue of a matrix which is about three billion by three billion in size. This is a ludicrous but real application of linear algebra. I have no idea how to effectively compute such a thing.
End [ADVT]

I tried again to discuss closest approximation, but I am willing to concede that this approach to Fourier series does not seem to resonate with students. Perhaps it is my fault. I may not have prepared the territory enough in linear algebra.

We will start with a function f(x) defined on the interval [-Pi,Pi]. In Fourier series, we try to write f(x) as a sum of sines and cosines. That is, I want to write f(x) as a linear combination of cos(nx) and sin(nx), where n and m are non-negative integers. As to cos(nx), the term corresponding to n=0 is somewhat special: cos(0x) is the constant 1. We can discard sin(0x) since sin(0x) is 0. There are various constants that come up which occur so that orthogonal functions are orthonormal. Here's a general approach which is useful in many integrals.

Integrating certain exponentials
The interval of interest is [-Pi,Pi]. I would like to evaluate -PiPiein x dx, where n is an integer. Please remember Euler's formula:
ein x=cos(nx)+i sin(nx), so that
e-in x=cos(nx)-i sin(nx), since cosine is even and sine is odd. Then adding, subtracting, dividing, etc., gives us:
cos(nx)=(1/2)(ein x+e-in x) and sin(nx)=(1/(2i))(ein x-e-in x).
Now -PiPiein x dx=(1/(in)) ein x|Pi-Pi. This is (1/(in))(einPi-e-inPi). Apparently this is sin(nPi) with some multiplicative constant, if n is an integer, this is 0. Hold it! There is exactly one integer for which the given antiderivative is not valid: n=0. It is impossible and forbidden to divide by 0. When n=0, -PiPiein x dx=-PiPie0x dx=-PiPi1 dx=2Pi.
-PiPiein x dx= 2Pi if n=0.
0 if n is an
integer and not 0.
Computing some trig integrals
Example #1 Let's compute -PiPisin(5x)cos(7x) dx. Certainly this integral can be done using integration by parts, twice. I will compute it the complex because it is easier.

We know sin(5x)=(1/{2i})(e5i x-e-5i x) and cos(7x)=(1/2)(e7i x+e-7i x). Then the product sin(5x)cos(7x) is (1/{4i})(e12i x+e2i x+e-2i x-e-12i x. What is the integral of this over the interval [-Pi,Pi}? Since the function is a linear combination of functions eAi x with A not equal to 0, and integration is linear, the integral is 0! Isn't this easier than several integration by parts?

Example #2 Let's compute -PiPicos(8x)cos(8x) dx. Again, cos(8x)=(1/2)(e7i x+e-7i x) so that [cos(8x)]2 is (1/4)(e14i x+2e0 x +e-14i x). Now the first and last terms integrate to 0 on [-Pi,Pi]. The middle term gives us (1/4)(2) on the interval [-Pi,Pi], whose integral is Pi.

Examples #1 and #2, revisited in the context of linear algebra
This example, slightly generalized, states that sin(mx) and cos(nx) are orthogonal on the interval [-Pi,Pi], when the inner product is the integral of the product of the two functions. Similarly, if m is not equal to n, the integral of sin(nx) and sin(mx), and the integral of cos(nx) and cos(mx), are also 0, so these pairs of functions are also orthogonal. The collection of functions {sin(nx),cos(mx)} for n>0 and m>=0 is not orthonormal, however, because the integrals of the squares aren't 1. The integrals of [sin(nx)]2 and [cos(mx)]2 over the interval [-Pi,Pi] (when n and m are at least 1) are all Pi. The integral of [cos(0x)]2 is 2Pi since cos(0x) is the constant 1. Notice that we don't need (we can't use!) [sin(0x)]2 because that is the 0 function, also known as the 0 vector in this setup.

Consequences
Here are some formulas which need to be memorized, or at least thought about a lot. The reason for the factors in front of the integrals comes from the computations mentioned just previously. Here it all is:  
Suppose f(x) is a piecewise continuous function on the interval [-Pi,Pi]
The Fourier coefficients of f(x)
     an=(1/[Pi])-PiPif(x)cos(nx) dx for n>=0
     bn=(1/[Pi])-PiPif(x)sin(nx) dx for n>0
The Fourier series of f(x)
This is (1/2)a0+SUMn=1infinityancos(nx)+SUMn=1infinitybnsin(nx).
 
The Fourier series tries very hard to converge to f(x).
  • The partial sums of the Fourier series turn out to be the best mean square approximations. Therefore, if you know some Fourier coefficients, and you want to approximate f(x) at random for x's in the interval [-Pi,Pi], the best that can be done is to compute (1/2)a0+SUMn=1Nancos(nx)+SUMn=1Nbnsin(nx). and your average error will be as small as possible.
  • If f(x) is, say, at least piecewise "smooth", so it is differentiable at all but a few points, then the whole Fourier series will try to converge to f(x).
  • Where f(x) jumps (has a discontinuity) it turns out that the Fourier series will still converge, but it will converge to the average of the left- and right-hand limits of f(x): so the sum of the whole Fourier series at x0 converges to (1/2)(limx-->x0-f(x)+limx-->x0+f(x))
         Original graph of f(x)                Graph of the sum of f(x)'s Fourier series
Please compare closely the sum of the Fourier series and the original function. They are the same where f(x) is continuous, and the sum of the Fourier series attempts to compromise (?) where f(x) is not continuous.

I also mentioned
Parseval's Theorem: (1/Pi)-PiPi(f(x))2 dx=(1/2)a02+SUMn=1infinity(an2+bn2).
This works because the Fourier series is actually sort of a linear combination of an orthonormal basis, and therefore the length squared of the sum ((1/Pi)-PiPi(f(x))2 dx) is the same as the sum of the squares of the coefficients of the linear combination (the other side of Parseval's equality).

I only showed one example. I need to show many more examples and show people how Fourier series are used.

One example
f(x)=x2-Pi2 on [-Pi,Pi]. Notice that this f(x) is even: f(x)=f(-x). Therefore any integral of f(x) multiplied by the odd function sin(nx) on [-Pi,Pi] is 0. We only need to deal with the cosine terms.

I used Maple. Here I defined a function to compute the Fourier coefficients:
q:=n->(1/Pi)*int((Pi^2-x^2)*cos(n*x),x=-Pi..Pi);
For example, the entry q(3) gets the response 4/9.

Here is the first (zeroth) term compared with the original function. This is the constant which has the least mean square average error from f(x) on the whole interval:
plot({Pi^2-x^2,(1/2)*q(0)}, x=-Pi..Pi,thickness=4,scaling=constrained);
Here are the sum of the first three terms, compared with the original function (I broke q(0) out from the others because it needs a different normalization):
plot({Pi^2-x^2,(1/2)*q(0)+sum(q(n)*cos(n*x),n=1..2)}, x=-Pi..Pi,thickness=4,scaling=constrained);
And the sum of the first 6 terms, which I find pretty darn close to f(x) on all of [-Pi,Pi]:
plot({Pi^2-x^2,(1/2)*q(0)+sum(q(n)*cos(n*x),n=1..5)}, x=-Pi..Pi,thickness=4,scaling=constrained);

The QotD was to find the Fourier coefficients of the function f(x) which is 0 for -Pi<x<0 and 1 for 0<x<Pi.

Please hand in 13.2 #9 and 13.3 #5 and 13.5 #2. Please read the applicable sections of chapter 13.

Thursday,
April 1
I tried to show some of the real linkage between linear algebra and Fourier series. I began, though, with the QotD: I asserted several times during the last lecture that matrix multiplication was not necessarily commutative. Find simple matrices A and B so that AB and BA are defined (compute them!) and AB does not equal BA.

I put the QotD first because I didn't know how long I would last in class, since I was miserable and suffering from the effects of an upper respiratory infection. One source declares

What are the risks to others?

Most upper respiratory infections are highly contagious. They are transmitted through respiratory secretions. Sneezing and coughing can spread these droplets. The germ can also be passed on when an individual with an upper respiratory infection touches his or her nose and then handles an object that another person later touches. The second person can then pick up the germ from the object and transfer it into his or her own respiratory tract by touching the face.


I tried diligently not to breath on anyone. Sympathy seemed meager.

Well, let's get to work. i began with an extended discussion of an example of a matrix I had analyzed with Maple. What I did was write a 5 by 5 symmetric matrix with small integer coefficients, and find the characteristic polynomial. I kept making entries equal to 0 to simplify the characteristic polynomial. I finally ended up with A=

( 0 -1 0 1  0)
(-1  0 0 0  0)
( 0  0 1 0  0)
( 1  0 0 0  0)
( 0  0 0 0 -1)
Since A=At, the general theory (almost all of which I didn't verify!) says:
  • All eigenvalues of A are real.
  • Eigenvalues corresponding to distinct eigenvectors are orthogonal.
  • Rn has a basis of eigenvectors of A.
The Maple command charpoly applied to A gives this result: 5-33+2. Then I can ask Maple to factor this polynomial, resulting in (-1)(+1)(2-2). Maple can be told to factor this further, using square roots, but this was enough for me. I then used the command eigenvects to get associated eigenvectors. Here is the result:

EigenvalueAn associated eigenvector
0(0,1,0,1,0)
1(0,0,1,0,0)
-1(0,0,0,0,1)
sqrt(2)(sqrt(2),-1,0,1,0)
-sqrt(2)(-sqrt(2),-1,0,1,0)

You can check (as I did, first by myself before class, and then quickly in class) that the list of associated eigenvectors are orthogonal.

Now I would like to create a matrix which diagonalizes A. We can normalize the eigenvectors and then create an orthogonal matrix. This means get eigenvectors with length=1, and get a matrix C so that C-1=Ct.

EigenvalueAssociated normalized
eigenvector
0(0,1/sqrt(2),0,1/sqrt(2),0)
1(0,0,1,0,0)
-1(0,0,0,0,1)
sqrt(2)(sqrt(2)/2,-1/2,0,1/2,0)
-sqrt(2)(-sqrt(2)/2,-1/2,0,1/2,0)

We must write the normalized eigenvectors as column vectors to create the matrix C:

(   0      0 0 sqrt(2)/2 -sqrt(2)/2)
(1/sqrt(2) 0 0    -1/2        -1/2 )
(   0      1 0      0           0  )
(1/sqrt(2) 0 0     1/2         1/2 )
(   0      0 1      0           0  )
then one can check that C-1=Ct and C-1AC=D, where D is a diagonal matrix, with diagonal entries in the same order as the corresponding eigenvectors of C:
(0 0  0    0        0   )
(0 1  0    0        0   )
(0 0 -1    0        0   )
(0 0  0 sqrt(2)     0   )
(0 0  0    0    -sqrt(2))
But this is all "review", things we've done before. I want to think about a slightly different problem, a problem in a context for which most people have little intuition (just means they haven't computed enough examples!).

Closest point or best approximation
Suppose we choose three of the eigenvectors of A, say v1=(0,1/sqrt(2),0,1/sqrt(2),0) and v2=(0,0,0,0,1) and v3=(-sqrt(2)/2,-1/2,0,1/2,0). Let's call S the collection of all linear combinations of v1 and v2 and v3. That is, S is all vectors in R5 which can be written as av1+bv2+cv3. Remember that the vectors are normalized eigenvectors of a symmetric matrix, so each of them has length 1 and they are mutually perpendicular. They are linearly independent. Therefore S is a 3-dimensional subspace of R5. Now suppose p=(2,3,4,-1,-7) is some "random" point of R5. What is the closest point in S to p? That is, I would like to find the best approximation, q, in S to p. Our criterion for goodness is distance, and closer is better. Can we find q in S so that the distance from q to p is smallest? We're considering distance from a point in R5 to a 3-dimensional subspace. Such problems may be difficult in general. Here, though, I would like to convince you that the problem is very easy.

Let's think about the distance from p to q. First, there's a square root in the formula. Well, if we want to minimize something, and that something is non-negative, we could also minimize its square. So we'll be happy to find the minimum of the distance from p to q squared. The distance is computed by taking the square root of a dot product, so we need to minimize (q-p)·(q-p). Since q=av1+bv2+cv3 we have to consider the function ([av1+bv2+cv3]-p)·([av1+bv2+cv3]-p). In this algebraic mess, v1 and v2 and v3 and p are fixed, and a and b and c are varying. So this is some complicated function of three variables. I think I called it F(a,b,c).

Dot product has some nice properties, like (av)·w=v·(aw)=a(v·w) if v and w are vectors and a is a scalar, and v·w=w·v and v·(wi+w2=(v·wi)+(v·w2). So now I can take the formula for F(a,b,c) and use these properties. Therefore F(a,b,c) is:
a2v1·v1+abv1·v2+acv1·v2+bav2·v1+b2v2·v2+bcv2·v3+cav3·v1+cbv3·v2+c2v3·v3
-av1·p-bv2·p-cv3·p-ap·v1-bp·v2-cp·v3 -p·p.
There are so many many terms. But ... but ... but: the vj's are orthonormal. So dot products between them are either 0 or 1. The first row above simplifies a great deal. And dot product is commutative, so the second row collapses a lot. The formula becomes:
F(a,b,c)=a2+b2+c2-2a v1·p-2b v2·p-2c v3·p-p·p.

That is quite a simplification. F(a,b,c) doesn't appear quite as menacing any more. In fact, because of the algebraic structure, I hope you can see that whenever any of the variables a or b or c gets really large, the value of F(a,b,c) will get large also, since the quadratic terms will eventually dominate any linear term. We can find the minimum (candidates for the minimum, maybe) by taking partial derivatives of F(a,b,c) and setting them equal to 0. Here we go:
The partial derivative with respect to a is 2a-2v1·p, and this is 0 only when a=v1·p.
The partial derivative with respect to b is 2b-2v2·p, and this is 0 only when b=v2·p.
The partial derivative with respect to c is 2a-2v3·p, and this is 0 only when c=v3·p.

We have solved the problem. The closest point, and the best approximation in S, to p is q=(v1·p)v1+(v2·p)v2+(v3·p)v3. The dot products are very easy to compute. Notice, though, that this problem was made easy exactly by the orthonormal nature of the original collection of vectors v1 and v2 and v3.

Do the same thing
Now I would like to describe something similar but in a very different setting. Again, this is likely to be unfamiliar territory, and please realize that the similarities I'll show took more than a century to recognize. Suppose f(x) is a function defined on an interval [a,b], and we want to approximate this function by another function, g(x). There could be several reasons for this. Maybe g(x) is easier to compute. Maybe g(x) is easier to analyze theoretically. Maybe f(x) might be "real-world" data, and it might be "noisy" and we might want to clean up the noise. So there could be many reasons. I would like the "closest" g(x) or the best approximation. Of course, I need to specify how we will test our goodness of appproximation. I'd like one number which could specify closeness. The integral of f(x)-g(x) is not used very much, because (look at the picture) the positive area (f above g) could easily cancel with the negative area (f below g) and give a deceptively small result for "closeness". What's another candidate? One is ab|f(x)-g(x)| dx. Here there won't be cancellation because of the absolute value. But absolute value has a corner (o.k.: it is not differentiable) and people perhaps don't like to use it for that reason. Instead, what is frequently used is something like ab(f(x)-g(x))2 dx. Sometimes this is divided by b-a, and it maybe is called the mean square error. This is a good definition of error, and leads to some really neat computation.

The textbook has formulas for more general intervals, but I think I will use the interval [-Pi,Pi]. The analog of inner product is this: take two functions f(x) and g(x) and compute -PiPif(x)g(x) dx. The most commonly used analogs of orthonormal vectors are sines and cosines.

Sines and cosines and 0's and 1's
I tried to examine the values of various integrals. For example, what is -PiPisin(x)cos(2x) dx? I used only ideas ab out even and odd functions. A function is even if f(-x)=f(x). So cos(2x) is even. A function is odd if f(-x)=-f(x). Here sin(x) is certainly odd. And [-Pi,Pi] is a balanced interval, as much negative as positive. It is easy to forget the following fact: the integral of an odd function over such an interval ([-Pi,Pi]) is 0. The product of an odd function and an even function is odd, so the integral of sin(x)cos(2x) is 0: these functions are orthogonal. I also tried to use pictures to convince people of this. An example
I gave an example which I will have to go over again on Tuesday. The example was f(x)=Pi2-x2 on the interval [-Pi,Pi].

Please hand in 8.2: 1,7 8.3: 1 at the next class meeting. Please study for the exam. You may also want to read the beginning of chapter 13 on Fourier series.


Maintained by greenfie@math.rutgers.edu and last modified 4/2/2004.