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 equations.
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 left.
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 difference sort of is that the true true growth factor G
-- tiy remember, the true growth factor was,
for the wave equation -- This is the G for the --
that multiplies the exponential e to the i*k*x
in the difference equation -- and G, the true G,
is e to the i*k*c*t.
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 diffusion,
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 --
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-Friedrichs-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 equation.
So now, we really have to think about better methods,
other methods.
And here are two: Lax-Friedrichs,
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-Fridrichs, 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-Friedrichs.
So what's the idea of Lax-Friedrichs?
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 in G and this space difference produced
an imaginary part in G, so we had G as 1 plus something
imaginary and we were outside the circle automatically.
Friedrichs thought of a way to get around that.
Instead of using U_(j, n) here, instead of the backward
difference, he replaced U_(j, n) by the average of that value
and that value, so that both sides now involved just two
numbers -- U_(j, n) -- this one and this one.
This one is not used now.
Again, 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-Friedrichs.
So U_(j+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, U_(j+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 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 i*k delta x plus this number times e to the plus
i*k delta x.
You're seeing, now, the pattern.
So I've sort of jumped to G without plugging in the e
to the i*k*x, taking a step and then factoring out e
to the i*k*x.
I did it without writing down the e to the i*k*x.
This is what factors out.
Maybe a 1 plus r would be more correct.
So stability 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 i*r times the sine.
This is really -- I could rewrite that as cos k delta x
plus i*r 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-Friedrichs-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 sign.
Like Friedrichs 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 -- r is between minus 1 and 1 --
see what's nice about Friedrichs 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,
we're smaller than 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 ellipse 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 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 U_(j+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 method 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-Friedrichs.
That's Lex-Friedrichs.
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 u_xx, right?
So this is -- they've thrown in a term that approximates u_xx,
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 -- Oh, maybe you could ask MATLAB to do
that.
Write down the G for 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-Friedrichs value.
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 small k, at least for 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 a while yesterday to rewrite 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 fourth 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-Friedrichs 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-Friedrichs 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-Friedrichs, which only used these two coefficients,
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 second-order accuracy.
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 a while -- 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 part,
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 cosine
theta -- that's of what order?
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 what Lax-Wendroff -- Lax-Friedrichs --
so Lax-Friedrichs, 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 centered,
but it's going to be spread out.
So that, I think, would be Lax-Friedrichs.
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.
All right.
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 -- do 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 Fourier series for a square wave -- well,
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 a few 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.
All right.
And 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 are 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 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 equation 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 du/dt as minus u_xx.
That would be not well-posed.
Or, if you tried to solve Laplace's equation equation --
say, u_tt plus u_xx 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 sign< and the heat
equation with the right sign.
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 require stability?
So that instability could not [UNINTELLIGIBLE]
So that's sort of maybe the fundamental theorem
of numerical analysis.
The fact that stability, a -- we better say what stability is,
and what is convergence.
So let me capture these in inequalities,
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 method?
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
to the next time.
So I'm going to say that true u, which is always a small u,
at delta t is some R, something R times u of t,
where the approximate guy is some 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,
then we've taken R n times, so the true solution is R
to the n-th times the initial value,
and the approximate guy is S^n times the initial value
and we'll -- we might as well start out with those the same.
So that's what R and S are.
Now, we got three ideas here.
What does well-posed mean?
Well-posed means that the true thing is bounded.
That R to the n-th 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 equation, 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 a bound on S to the n, say, C_2.
And again, in our example, C_2 has been 1.
This is like -- C_2 is like -- I mean,
S to the n-th on each exponential is this G
to the n-th, 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 S u, 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 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 bound 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 k 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 propagated up to any later time n.
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 n-th and S to the n-th.
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 n-th minus S to the n-th is small?
That's my job.
This is the error when I apply it 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 one 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 an 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 R's to come first.
I want maybe R minus S times S to the n minus 1 and then --
sorry.
I want R's on the right.
Good.
Finally got it.
So that gives me the R of the n-th, but I've subtracted off,
so I need an S R minus S R 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 is 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 1's and that will be T -- n 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 conclude with Lax's insight
on the other part.
The fact that convergence requires 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 --
this is 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 i*k*x.
But, we didn't follow this point,
but if I take an e to the i*k*x where G is bigger than 1 --
of course, that is getting amplified a bit,
but for that particular e to the i*k*x,
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 i*k*x,
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.