Lecture 4: Comparison of Methods for the Wave Equation

Flash and JavaScript are required for this feature.

Download the video from iTunes U or the Internet Archive.

NARRATOR: The following content is provided by MIT OpenCourseWare under a Creative Commons license. Additional information about our license and MIT OpenCourseWare in general is available at ocw.mit.edu.

PROFESSOR: So, ready for the second lecture about -- second to last lecture about the linear one way wave equation, so the model for hyperbolic equation. That word hyperbolic tells me that signals are traveling with finite speed. In this case, the speed is c and in this case, the signals are only going to the last. The two way wave equation will travel left and right, but again, finite speed. I emphasize finite speed, because for the heat equation, the speed is infinite. As soon as you change an initial condition, at any later time, you immediately feel the effect -- probably a small effect way, way out. So diffusion travels with infinite speed. Waves travel with finite speed.

The different sort of is that the true -- the true growth factor g -- you remember, the true growth factor was, for the wave equation. That This is the g for the -- that multiplies the exponential e to the ikx and the difference equation -- and g, the true g, is e to the ikct. So that -- what I was going to say is, that has absolute value 1. Again, a difference from the heat equation. For the heat equation with the fusion, g will be smaller than 1. You're losing energy. For the wave equation, the energy is constant and the wave is just moving.

So it's a little more delicate, because we're trying to -- this g in the difference equation has to stay -- still be somehow close to that one, which is right on the unit circle and yet, we got to stay -- we can't go outside the unit circle. So here are two methods. We had one method last time that worked. That was one sided differences and of course, you had to take the difference on the good side, the upwind side and that was first order. The Courant-Freidrichs-Lewy condition was the same as the von Neumann stability condition and it was r between 0 and 1. r is this ratio, c delta t over delta x, it comes in. It comes into the wave [? limit. ?]

So now we really have to think about better methods, other methods and here are two: Lax-Freidrichs, which will only be first order. So when I -- these are the questions I always ask, and then finally, experiment. I mean, answering these questions will certainly be a guide to what we'll see experimentally, but there's more to learn from actually doing the numbers. The accuracy, Lax-Freidrichs, will be only first order, whereas Lax-Wendroff will be second order. So with respect to accuracy, that's certainly better. Can we see that maybe without too much calculation? The notes would have the more details.

Can you see -- here's Lax-Wendroff, Lax-Freidrichs. So what's the idea of Lax-Freidrichs? You remember the total failure, the complete and absolute failure when we used a time difference there and a space difference there? You remember that that gave the time difference that produced a 1 ng and this space difference produced an imaginary part ng, so we had g as 1 plus something imaginary and we were outside the circle automatically. Freidrichs thought of a way to get around that. Instead of using ujn here, instead of the backward difference, he replaced ujn by the average of that value and that value, so that both sides now involved just two numbers -- uj and uu, this one and this one. This one is not used now.

Again, for a still first order, you won't be surprised that this approximation, this has got first order error. But stability now is back with us and easy to check because now each new value depends on what old values -- just these two, and with what coefficients? So let me just check stability for Lax-Freidrichs. So uj plus 1 n, the new value, is some multiple of this guy, the one on the left and some multiple of this one, the one on the right. What are those multiples? Let's see. I'm going to move things over to the right hand side of the equation, so I think there's a half and then I think there's an r over -- minus, I guess. This is the one that comes in with a minus sign -- it's the one on the left and the one on the right, uj plus 1 n. There's a half plus r over 2. Simple idea. Simple idea and it has one important advantage over upwind, one sided that we saw last time. It's not a great method. It's only the first order here. If I -- and you see when it's going to be stable? The coefficients add to 1. Our coefficients always add to 1, right? Because if we have a constant state, this next step should be that same state. It's not going to change, but instead of having the typical 1 in the g, I have -- it's split into a half, so what is g for this guy? g for this guy is this number: 1 minus r over 2 times e to the minus ik delta x plus this number times e to the plus ik delta x. You're seeing now the pattern.

So I've sort of jumped to g without plugging in the e to the ikx , taking a step and then factoring out e to the ikx. I did it without writing down the e to the ikx. This is what factors out. Maybe a 1 plus r would be more correct. So civility requires this complex number, which I could write in different ways. The real part there would be a cosine. This would actually be the cosine plus r times the sine -- plus ir times the sine. This is really -- I could rewrite that as cos k delta x plus ar sine k delta x. I don't know which way you -- here, this shows you the real and imaginary parts. That makes it nice. This shows you the two coefficients. That makes it nice. Can you see stability -- so what is Courant-Freidrichs-Lewy? The characteristic discussion, the argument based on characteristics, says what? Stability is going mean that delta t -- stability is going to happen when r is below 1. It's just like upwind, except that r could also be coming -- c could have the opposite sine. Like Freidrichs is prepared for c to be positive or negative, because it's not purely upwind or purely downwind. It's both winds.

You see that if r is anywhere -- the condition is going to be that minus 1 is -- but r is between minus 1 and 1 -- see what's nice about Freidrichs in that case. This number is positive. This number is positive, because the 1 is bigger than the r. So each step -- and they add to 1. So each step takes the new value as a combination of the two old values and of course, the magnitude can't be bigger than 1. So under this condition, I have positive coefficients. That's positive and that's positive. They multiply things of absolute value 1, so they can't -- the sum can't get beyond 1 -- or you might like it here, that magnitude squared is cosined squared plus r squared sine squared, and if r is below 1 where smaller then cosine squared plus sine squared, which is 1.

So if I draw the -- just for a moment, if I draw the g in the complex plane -- so there's the unit circle. k equals 0 is always there. That's at 1, then these points -- I think probably those points fall on an elipse here, where we go up as high as r or magnitude of r and we go left, right as far as -- this is k delta x equal pi as far as minus 1. So it's stable if this is less than 1. Nothing hard. In fact, easy. Two positive numbers that add to 1.

Now I've added a problem in the problem set. Not requiring it to it be turned in, but just to mention it here -- that whenever the coefficients are all positive, the accuracy could only be first order, so we can't get too far with positive coefficients. First order accuracy is the best possible, except in the extreme case when we're right on the characteristic. So if r was exactly 1 -- you remember that special case when r is exactly 1 and we're following exactly the line that carries the true information -- r is 1. That's gone and this is just -- so if r is 1, then uj plus 1 n is u -- is this, is that number and we're just carrying that information exactly. I have to say, again, that the idea of trying to follow that characteristic is a totally natural one and I'm not writing down here a message that tries to do that, but it wouldn't be a bad idea in one dimension. You couldn't do it, really, well in higher dimensions where these ideas do continue to apply. So that's Lex-Freidrichs. That's Lex-Freidrichs.

So now let me turn to Lax-Wendroff. These are natural ideas, but Lax and Wendroff did it first. They said, if I use these three values, all three, then certainly I can get the accuracy up one more because I have one more coefficient to choose and to see what they did, the easiest way is to write it as a time difference on the left to see how they match. All they want to do with match that second term, that delta x squared term in the Taylor series. Because this is going to give us a d by dt, this will give us a d by dx centered, so we get extra accuracy there and the matching that we need is this term. Do you recognize that? What does that term approximate? Hope I've got it right. That term approximates uxf, right? So this is -- they've thrown in a term that approximates uxx and its purpose is to catch the next order of accuracy and hopefully not lose stability, hopefully remain stable.

So what Lax-Wendroff are doing is staying even closer to this curve. I guess I have forgotten exactly what their curve will be. If I write down g of, maybe you could ask MATLAB to do that. Write down the g 4 -- this for Lax-Wendroff. It's got that extra term and plot absolute value of g. I guess you might plot it along with the Lax-Freidrichs out. So what am I guessing? I'm pretty sure that the Lax-Wendroff one will hug the unit circle more closely, but it'll stay inside. So can I just -- Lax-Wendroff is g closer to 1, at least for a small k, at least for a small k. That's really where the accuracy is checked. How close does these Taylor series match? Lax-Wendroff match better and they stay inside. That takes a little calculation. Took me awhile yesterday to rewriter the -- write out the g, square the real part, square the imaginary part, add them up and this is what I got. So for Lax-Wendroff, I got g squared is 1 minus r squared minus r to the 4th times 1 minus cosine k delta x squared. I put in the notes eventually, much algebra, but it does the job. Because if r squared -- so what's the stability condition? It's again this same guy, the same condition -- r squared below 1.

Again, we could be -- we could have a left going wave and a right going wave. Now you'll say, what's the good of that? Here we just have one wave, but on Friday, we'll have a wave going both ways. So Lax-Freidrichs and Lax-Wendroff can both deal with that. In fact, their stability conditions are the same and in fact, you may be wondering why would anybody use Lax-Freidrichs when Lax-Wendroff gives higher accuracy at no real significant increase in cost? So Lax-Wendroff is a favorite. Don't forget, still experiment has to be done to see what it really does in practice, because right now, everything's looking in Lax-Wendroff's favor with one exception. Lax-Freidrichs, which only used these two coefficients and had them positive. This is positive. This is positive, where Lax-Wendroff is going to have a negative.

Let's see. Where will our negative -- I think there'll be a negative coefficient in here somewhere. Probably one of the three outside guys. I know that we couldn't have three positive coefficients and get up to the second order actors. Now where do we pay the price for that? Where do we pay the price -- I should finish here. Suppose r squared's bigger than 1. If you accept this calculation, which has a nice result, but took awhile -- if r squared's below 1, then what do I know about g? Suppose r squared's bigger than 1. What can we say from this formula? Does that hit you right away? g -- if r squared -- this is certainly -- absolute g. I've taken the actual g and squared its real parts, squared its imaginary part, did all the algebra I could think of and got this answer. For the magnitude squared -- so the magnitude squared is never going to go below 0. I don't have to think of that. It's just I want to be sure it's below 1.

So all I really want to do is be sure that I'm subtracting something positive. Am I subtracting something positive? That's certainly positive. Is this positive? r squared'll be bigger than r to the fourth. If r squared's less than 1, then r to the fourth is even smaller. So this is positive, this is positive. I'm subtracting something so g is certainly not bigger than 1. So that formula really makes it evident and I guess also evident that we're getting some high accuracy. You see, what's the size of this when k is small?

When k is small, if I'm taking the cosine of a small angle -- what am I looking at there? If this is a small angle theta, I'm looking at 1 minus cosign theta -- that's of what order? Theta squared. Theta squared. Then that's squared again, so that I'm up to, theta to the fourth. That's why I'm getting off to a good start. So stability we've got, two sided we've got, experiment -- can I just try to draw a little graph of what I think the two methods will do for our model problem? Like this model problem of a wave moving along? I want to try to draw it at, let's say, t equal to 1, for example. So at t equal 1, the wave should be here. I'll draw -- you might say this is totally unsatisfactory, and it is, because I'm just drawing by hand what you can see. Maybe you should. This would be a perfect experiment and I really could use for the notes -- the notes could use a proper graph.

I think what you see -- let's see. At t equal 1, both methods are going to be 0 way out here and they'll be 1 way out here. Why's that? Because the signal speed is finite, right? This way out here, it hasn't heard about this wave of water. Way out here, it doesn't know that it's not a constant. The question is, and I guess the thing has moved, so maybe -- can I try roughly? So the true wall of water would be that. I think that of what Lax-Wendroff -- Lax-Freidrichs -- so Lax-Freidrichs with positive coefficients will show a smooth profile, but the perfect signal will get spread out.

The perfect signal will get spread out. Why is that? Let me just draw a spread out signal. So apologies to do this by hand. I think it's probably going to be pretty standard, but it's going to be spread out. So that, I think, would be Lax-Freidrichs. That's very rough, right? But I really -- hearty thanks for anybody who gives me a code that does it and graphs it. But the main point is, we have a lot of smearing here. The accuracy is not going to be terrific. But there isn't any bumpiness. The reason there's no bumpiness is that we're Lax-Freidrichs with -- has these positive coefficients, so that every value always has to stay between 0 and 1. The starting values were between 0 and 1. At every step, the values are between 0 and 1, because we're always taking a combination of a positive multiple of one value and another positive multiple of the other value and the two coefficients adding to one. We're taking 30% of one and 70% of another and we never -- we're always in between the 0 and the 1.

What about Lax-Wendroff? What are you expecting for Lax-Wendroff? Expecting better accuracy, expecting -- if I were on the ball, I'd have another color. Lax-Wendroff will be a much better profile, sharper profile, but you know what the price is going to be? There'll be some oscillation. There'll be some oscillations behind -- I think behind -- oh gosh, I should know. Forgive me for -- I can fix these if you'll tell me how on the next -- because I'm not sure that I should know whether they're at both or whether the oscillations happens at both ends. I doubt it. No, I doubt it. I don't think it will. Why do I not think it well? Because in Lax-Wendroff, the three coefficients -- one of them will be negative and two will be positive. They'll add to 1, of course. So the oscillations, which are coming from sometimes a negative coefficient, sometimes a positive, will all be on one side. So only one of those oscillations [UNINTELLIGIBLE] So this is my picture of -- can I -- of Lax-Wendroff with oscillations behind the wave front. Alright.

If you've done any computing, you, when you see these oscillations, you immediately think of somebody's name, which is Gibbs. Of course, Gibbs -- you remember the Gibbs phenomenon that comes when you're approximating a square wave, like there by pure Fourier series? So the Gibbs, the classical Gibbs phenomenon is, take the 4 a series for a square wave while it's smooth and then truncate that series, cut that series off at 100 terms or 1,000 terms. Gibbs noticed -- and it's always a problem -- that these oscillations appear right near the front. So far away, no problem, but here we really see it only on one side, because we're not doing Gibbs Fourier series experiment. We're doing Lax-Wendroff wave equation. When I say we, I really hope you'll do it and new lines of MATLAB will be greatly appreciated. Maybe Mr. Cho can -- you could email to him, for example, a code and the plot.

So what does that tell us? Tells us Lax-Wendroff has got a lot of potential, but this oscillation isn't good. How are you going to get rid of it? How do you get rid of bumpiness? You introduce diffusion. So diffusion, heat equation stuff smooth things out. Wave equations don't smooth things out. This wall of water in the wave equation stayed a sharp wall. But we'll see very soon in the heat equation, diffusion, viscosity, all those physical things connect numerically to a smoother solution, smoother profile. So that is really a key idea, is -- and we'll see it, but it will make the problem not linear. If we want to stay with a linear equation and increase the accuracy, we're going to get oscillation. But if we locate the oscillation, introduce some diffusion in the right place, we can reduce the oscillation. So that's where -- this is what you really have to do if you want good output from numerical methods. Even on this simple model problem, you'll have to have this idea of numerical viscosity. Artificial viscosity, numerical diffusion popped in at the right spot. So we'll see that. Maybe the right place to see it will be for when the equation itself is already nonlinear. THat's where it has developed a lot. So we'll -- when we get to nonlinear wave equations, which isn't far off, that's what we'll see.

So now I guess I have one more topic in this lecture. that is convergence. How do you show that -- we know, we sort of know that we need stability. So I guess I'm going to have to come to the question of convergence. Alright. Of course, I have to mention -- more than mention -- the key result in the whole theory, the connection -- so we have methods that are first order or second order -- stable or unstable, and then the key idea is this Lax -- again, he enters -- Lax equivalency. So this is a little bit of pure math, comes into this and it gives exactly the result we would expect. Stability implies convergence and convergence implies stability. So that's the two ideas or equivalent, provided we're talking about a decent problem, for a consistent approximation to well-posed problems.

I've got four words on the board, four keywords words here. Stability, we kind of know what that means. Convergence. These are the new words and I'll just -- and actually -- so consistent means that the method is at least first order and that the difference equation approaches the differential equation at a single step. Well-posed means that the differential equations is OK. So, for example, what would be a non well-posed problem? When we get to the heat equation, if you ran the heat equation backward in time, that would give us an equation dudt as minus uxx. That would be not well-posed. That would be not well-posed. If you tried to solve Laplace's equation equation -- say, utt plus uxx equals 0 from initial values. That would not be well-posed. But we're dealing with well-posed problems. The wave equation, which has the right side and the heat equation with the right side. So our problems are going to be well-posed. You can safely think we're only talking about consistent approximations to well-posed problems.

Here is the main point. Stability holds if and only if the convergence holds. So there's actually two things to show there. How does stability give convergence? Why does convergence requires stability? So that instability could not [UNINTELLIGIBLE] So that's sort of maybe the fundamental theorem of numerical analysis. The fact that stability, a balance -- we better say what stability is and what is convergence. So let me capture these x in any quality, in mathematical terms rather than just words.

So those are the ideas -- maybe I'm going to have to write a little above here, but -- so what does -- what's our methods? So I need some notation and then I'll lift the board and we'll see what each thing means. So our method, we're allowing any difference method that moves forward in time. There's time. Let me use r for the true operator that takes the function, follows the differential equation for the next time. So I'm going to say that true u, which is always a small u at delta t is sum r, something r times u of t, where the approximate guy is sum s times u of t. So r is the true solution. So well-posed. So of course to take -- if this is n delta t and we've taken r n times, so that the true solution is r to the nth times the initial value and the approximate guy is sn times the initial value and we'll -- we might as well start out with those. That's what r and s are. So now we got three ideas here. What does well-posed -- well-posed means that the true thing is bound. That r to the nth times u is less than some constant times u.

These inequalities are sort of part of the mathematics of partial differential equations. I think you've got to take them as not -- I mean, take them as sort of a reasonable statement that the -- this is energy this u typically will measure the energy and in our wave equations, c 1, the constant c 1 will be 1. The energy doesn't change. In another problem, we could allow the possibility of c 1 being a little bigger than 1, but think of c 1 as 1. Then stability will be abound on s to the n, say, c 2. And again, in our examples, c 2 has been 1. This is like -- c 2 is like -- I mean, s to the nth on each exponential is g to the nth and g is below 1. So often c 1 and c 2 will be 1. Then finally, what is this ideas of consistency? And this means -- consistent means that -- so what's consistency?

Consistent is comparing the differential equation with the difference equation. It's comparing r with s and at a single step, they should be close. So at a single step, r minus su -- that's the -- starting with the same u, I'm just taking one step in the differential equation and comparing it with one step in the difference equation and I would want that to be c 3, maybe delta t to the p plus 1. Was that -- is that -- times some size of u, I think. Order delta t to the p plus 1. Like delta t squared. This would be the usual situation, that when it's stable. So at the moment then, I want to say, if these are all true, then I get convergence because what I'm now going to do is show that stability implies convergence. It's what we did for ordinary differential equations. This was the local error. This is the local error that I called de, the discreditization error. These are saying -- so this is saying that the -- here's really what happens. So I want to look at the error up here after n steps. So the error after n steps -- typical term is, I follow the differential equation. That's u. Then I have this little error de, which I make at that step, the difference between r and s and then I follow the difference equation up to the remaining time. What I want to show is that I can track the error and bounce it. I don't know if you remember the discussion for ordinary differential equations, but it went the same way. We had -- we looked at the new error de that was produced at each time step k. When we took steps of -- when we got up to this point with the true solution, we plugged the true solution into the difference equation, it had a little error and whatever error there was, that got thrown in with other errors and propogated up to any later time s. Am I bringing back what we did there or just -- you can just see that's what we're going to have here.

Now how will I express that? Somehow I have to estimate the difference between r to the nth and s to the nth. That's really my problem. Actually, it's a nice problem. Suppose I have control of powers of r and s, and I know that r minus s is small somehow. How do I show that r to the nth minus s to the nth is small? That's my job. This is the error when I applied to a starting value u. This is the difference between little u and big u. So how do I -- what I want to do is split it out into this situation where I just have 1 r minus s. I think this'll work. So suppose I have -- I think I could start with r to the n minus 1 times r minus s, and then in r to the n minus 2 times an r minus s times an s. Maybe I don't -- sorry, this guy is going to carry and I'd rather do -- sorry. Got to do it right because I want the rs to come first. I want maybe r minus s times s to the n minus 1 and then -- sorry. I want rs on the right. Good. Finally got it. So that gives me the r of the nth, but I've subtracted it off so I need an sr minus sr to the n minus 2 and so on.

Sorry for this bit of algebra in the middle of practical things, but do you see that this is the case where I followed the differential equation up and at the last step I made an error? This is the case where I followed the differential equation almost to the top, I made an error and then I had one more step to go with the difference. The s doing this part, so a typical term will be a power of s, s to some power, r minus s, r to some power. I'm breaking the error at n steps into n different terms, each from the error at one step. This is the error at the last step. This is the error at the next to last step and so on. We've got it. Because I just estimate every term, every term has a power of r, which gives me a c 1, a power of s which gives me a c 2 and the discreditization error, the one step actual error, r minus s, which gives me a c 3 and this delta t to the p plus 1. You see that I'm just piling up n of these, so I'm getting n of those, delta t to the p plus 1, so the final result will be, I'll have n of these delta t to the p plus 1s and that will be t and delta t, delta t to the p. So that'll be the error. Since I went to all the trouble of defining a c 1, c 2, and c 3, those three numbers. So the error -- there'll be a c 1 times c 2 times c 3 constant there in the bound. Have a look at the notes for that and let me just concludes with Lax's insight on that other part.

The fact that convergence require stability. So you see, what I've shown is that if we're stable, if I bound these, bound these, bound these, then I have a bound on the error. That's convergence. The other way is, if it's unstable, how do you know it couldn't converge? How do we know that convergence requires stability? It's a little delicate, so I just -- it's the last 30 seconds of the lecture, just to mention this point, because it's easy to miss a key point here.

We've tracked a different -- now I'm supposing that the method is unstable and I want to show that it will fail, right? That's my job now. That's Lax's -- other half of Lax's theorem: If it's unstable, the method will fail. There'll be some initial value for which it blows up. Now you might say, that initial value is just e to the iks. Buw we didn't follow this point, but if I take it e to the iks where g is bigger than 1 -- of course, that is getting amplified a bit, but for that particular e to the ikx, when delta t decreases -- so if I have a -- you see, this is a k delta x. What I'm trying to say is that for a fixed k, even an unstable method can work. For a fixed frequency k, as delta t and delta x gets smaller -- these guys, whether they're outside or inside, will be coming down toward 1 and if you're patient, you could see it actually worked. So it's somehow the initial function that doesn't work, that blows up, is a combination of frequencies. A single frequency, actually, you move down here and it's not catastrophic in the end, but a combination of frequencies, so that as one frequency moves down, another one is moving into this position, another one into this -- there is a typical initial function and it is disaster. So that's what Lax had to catch. The notes will use the single word, uniform boundedness.

Anyway, this is one point where pure math -- a kind of neat result in pure math says that if stability fails, then there is some single starting value on which it will fail. That starting value won't be a pure e to the ikx, because it, as we said, actually in the end gets to be OK, but there will be a combination of frequencies that'll fail and so convergence will fail. So that's the key theorem. I probably won't use the word theorem again, but this is one time when it's justified and when a little sort of abstract functional analysis gives the conclusion that convergence requires stability.

So you'll see that in the notes and next time we move on to the true wave equation, second order. Thanks.

Free Downloads

Video


Audio


Subtitle

  • English - US (SRT)