INTRODUCTION: 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: OK.
So I'll start now with the third lecture.
And we're into partial differential equations.
So we had two on ordinary differential equations and the
idea of stability came up, and the order of accuracy came up.
And those will be the key questions today --
accuracy and stability.
We do get beyond that actually.
I mean to look at the contour at the output from a
difference m once you've decided yes it's stable, yes
it's accurate.
You still have to look at the output and see, you know, how
good is it, where is it weak, where is it on target.
But first come the accuracy and stability questions.
And beginning with this equation.
So this will be the simplest initial value problem
we can think of.
I call it a one way wave equation.
Notice the two way wave equation, which we'll get too,
would have second derivatives there.
And that will send waves in both directions.
This we'll see sends a wave only in one direction, so it's
a nice scalar problem.
First order constant
coefficient, that's the velocity.
Perfect.
The model problem for wave equations.
Because it's just first order, I just begin it with the
function at time 0.
So I'm looking for the solution at
time t to this equation.
OK.
So that will not be hard to find.
And I want to do it first for pure exponentials.
When u of x0, the initial function, is a pure
exponential e to the ikx.
Fourier is always around.
And to see what happens suppose this is given, for
example, so this will be my example, e to the ikx.
OK.
Now what's special about e to the ikx?
The special thing is that with constant coefficients and no
boundaries, the solution will be a multiple of e to the ikx.
In other words, we can separate variables.
Maybe I'll put it here.
I can then look for the solution u to be a function
times a number really, but it will depend on the time on the
frequency k and the time t times e to the ikx.
You see how that separated variables? x is separated from
t, and the frequency controls what this
growth factor will be.
So this growth factors is a key quantity to compute.
And we compute it just by substituting that hope for a
solution into the equation.
And getting an equation for g, because e to
the ikx will cancel.
That's the key point, that everything remains a pure
harmonic with that single frequency k.
So if I plug it in, I get dg dt.
I'm taking the time derivative of u, times
e to the ikx, right.
That's du dt.
Because t was separate from x, that's what I get when I
plug in this u.
So I'm always plugging in this u.
What do I get on the right hand side? c, the x
derivative.
So that's just g.
The x derivative will bring down ik, e to the ikx.
No surprise that I can cancel e to the ikx,
which is never 0.
And I have a simple equation for g.
Linear constant coefficients, but the coefficient depends on
k, of course.
The coefficient is ick, and of course the solution is g is
equal to --
we're back to our simplest model of an ordinary
differential equation.
Now the coefficient a is now ikc, so it's e
to the ikct, right.
That shows us what happens to the e to the ikx.
So now I'll put that here, because now we know it.
So we have all together e to the ikct and here's an x.
So let's put those together as x plus ct, because that's a
guide to what's coming.
So that's the solution: separated variables, look for
the growth factor, tried an exponential.
And we'll do exactly the same thing for difference method.
So this was von Neumann's brilliant idea, watch
exponentials, watch every frequency, and see what
happens to the multiple of e to the ikx.
Now here is very special, because all frequencies are
producing this combination x plus ct.
So when we combine frequency, that's what Fourier said, take
combinations of these things.
The solution will be a combination of these things.
And what do we get?
If our original function was u of x and 0, it's a combination
of these, e to the ikx's.
And at later time we have the same combination by linearity
of these guys.
All that's happening is x is getting changed to x plus ct.
So I'm getting the answer.
Here was the answer for one exponential.
But now I'm going to do the answer for a general u, and
it's going to be easy to pick off.
It's just like picking off this x plus ct.
It will be u at x plus ct at time 0.
That's the solution for every u.
Not just for the exponentials, but for all combinations of
them which give us any u of x and 0.
So u of x and 0 was a combination of e to the ikx's,
and then the solution is the same
combination of these guys.
Well, we can understand what that solution is.
It's what I call a one way way.
Let me understand what's happening here.
So here's the key picture now.
I want to a draw a picture in the xt plane.
So I'm not graphing u here.
This is a picture in the xt plane that will show us what
this formula means.
Because this is the solution that we then want to
approximate by difference method.
OK.
So what does this mean?
This means that the value, for example, u of 0, 0.
So there's the point 0, 0.
Here's what's happening, at a later time, say up here some
later time, the value of u, this is the line, x plus ct
equals 0; x plus ct equals 0 along that line.
And what's the point on that line, the solution is the same
all the way along.
So whatever the initial value is here it travels, so that u
at this point is u at this point.
So let me to call that point P, and say what I mean now.
I mean that u of P is the same as u at 0, 0.
Because x plus ct is the same here and here.
Let me take another point.
Let's say x0 another initial value, and again along the
line like this.
This will be the line x plus ct equals capital X, because
it's the line where x plus ct is a constant, and the
constant is chosen so that it starts at this point at t
equals 0, X big X. So now -- let's see --
I better make clear -- that this lines carries this
initial value.
This line carries this initial value.
You see it?
The value of u all along this line is the value it started
with right there.
And you know the name for those lines?
So these are lines, they happen to be straight here
because we got constant coefficients, straight and
even parallel.
And they're the lines on which the information travels.
Whatever the information was at times 0.
By information I mean the value of u at 0, 0, or the
value of u at this point.
That information travels along that line to give this
solution all along the line.
And the name of the line, anybody know it?
Characteristic.
So this is a characteristic line, they all are.
This to.
All these line are characteristic lines.
Can you read characteristic?
I should of have written it a little better.
So this is like a key feature of wave equations.
It's a key feature of wave equations that we will not see
for heat equations.
And so let me just say it again, because it will bear
also on the stability question for difference equations.
The information, the true information, for the true
solution travels along characteristic lines.
We're in one dimension here, and which is a major
simplification.
In two dimensions we'll have characteristic cones.
In three dimension certainly we'd call it a cone.
And actually I guess what I'm saying is that if I speak a
word, you know, my voice or the sound waves solve the 3D
wave equation.
So when I snap my fingers, that sound travels from this
point out along characteristics.
And of course it travels in all directions, so we've got a
more interesting, more complicated picture in higher
dimensions, but we really see it here in one dimension.
OK.
Now I should also include a graph of the solution.
What does the solution look like?
Maybe I'll find a spot here for that graph.
Now I'm going to graph u.
This is still x.
So I'll graph u at maybe two different times.
So the starting time.
Suppose we have a wall of water.
So this will be u of x and 0, a step function, maybe at 0.
So u of x and 0 is 1 and then 0.
OK.
So there's a typical starting initial value.
Well, I don't know if it's typical.
Pretty special to have 1 and then 0, a
perfect step function.
Now I want to ask what's the solution to the equation?
What is u of x and t?
I want to graph u of x and t.
What what happens?
You have to assemble sort of this picture, which is in the
xt plane, into a graph here of u of x and t.
What happens?
Well all these 0 starting values travel along
characteristics, and tell us that u is 0 along all these
characteristics --
u is 0.
This here, the starting value is 1.
So along these characteristics we're starting with a 1, we
end with a 1.
Do you see what happens?
How do I graph it?
The wave is moving.
The wall of water is moving that way, moving to the left.
Because I'm taking c to be a positive number.
So u of x and t is a wall, or a wave you could say, moving
left with speed c.
Every solution to the equation is doing this.
And in this particular case it's a wall of water and it
stays a wall of water.
That the point.
Shapes are not changed as we go.
There's no dispersion would be the right word I think.
We would expect in a general equation that the shape of the
wave would change a little as it travels.
In this equation it doesn't change.
The shape stays the same.
See again I'm thinking of this wave as a combination of pure
exponentials.
And they all travel with the same speed c, so the whole
wave travels with speed c.
And at a later time that's the graph of the solution.
OK.
So that much is all we have to do, all we can do really, in
solving the differential equations.
It's a simple model certainly.
But now I come to the difference equation.
Well when I say the difference equation, I should say
equations plural.
Because here I've written four possibilities.
Three of them are actually used, and one is a disaster.
But none of the four is perfect.
Most people would say Lax-Wendroff is the best of
these four, because it has one extra order of accuracy.
Oh yeah, yeah.
I guess these all have only first order accuracy.
But Lax-Wendroff moves up to second.
So that's a method that everybody knows, because it
gets up this way.
But actually upwind is the first one to start with.
Let me start with that one.
OK.
So I've multiplied up by the delta t that should be in the
denominator there.
I'm taking the standard forward difference in time --
all these methods are going to be explicit , so they're all
going to be solvable in a simple way for
u at the new time.
So a forward difference in time for du, dt and c times a
forward difference in space for du, da.
So that's the first idea that would occur to you, and it's
pretty reasonable.
We can check its accuracy, but you can guess
what it would be.
What's your guess on the accuracy of this method?
You're seeing a forward difference in time, a forward
difference in space.
The odds are those finite differences give you an
accuracy p, if we use p as the measure of accuracy, which
would be probably 1, first order accurate.
So we'll see the accuracy as only p equal to 1.
OK.
What about the stability?
I'm looking at this first method.
So one comment about the method is that this ratio c
delta t over delta x, let me call it r,
that's a key number.
it's called the Courant number by some.
So in the literature you might see the phrase Courant number
for that ratio r.
And maybe my first point is stability is going to put a
bound on r.
We can tell right away for this explicit method and for
any explicit method, that r can't be unlimited.
There's going to be a bound on r, which it means a
bound on delta t.
It's telling us that delta t has to be less than some
constant times delta x.
And stability is the job to find that constant.
How small does delta t have to be?
Now let me sort of understand this forward difference method
in this picture.
So I'm going to do the forward difference, delta t of u
equals that ratio times the forward difference of the
space difference, so time differences is that ratio
times this space difference.
OK.
Method 1.
OK.
So I guess somehow you want to see want to see a picture like
this for the difference equation, So the difference
equation we imagine we've got steps delta x, all these steps
are delta x.
Let's suppose that's 0.
OK.
So what happening now, and time going this way and I'm
solving this equation?
What's happening here?
Let me copy this method, u at a position j, so this might be
the point j delta x, at time n plus 1 is.
I want to write this equation as always in the most
convenient way.
So what am I going to do?
I'm going to take this u of x and t and put it
on the right hand.
So then my equation looks like that ratio r times u at j plus
1 at time n.
That's this.
I'm one position over, but I'm at the given time.
And then what's the coefficient of ujn?
So can you just tell me what do I put in there?
So this will be that nice way to look at
our difference equation.
Each new value is a combination of two old values,
r times the guy here.
So in other words, here's my little molecule.
This is the new value it comes from r times this plus what,
you have to tell me what multiplies
this you at this point.
So here's uj n plus 1, and here's uj plus
1n, and here is ujn.
And what goes in parentheses?
1 minus o.
Because I get the 1 when that flips over to the other side.
Good.
So it's just a combination.
The new value is a combination of those two.
And I guess what I want to show is that it could not
possibly be stable unless r is less or equal to 1.
So this will be the Courant condition.
The condition on the Courant number.
So the Courant condition or the CFL, that's a really a
better condition.
Let me say what it is.
The CFL condition is r less or equal to 1.
That's a condition for stability or convergence that
comes to us from comparing the real characteristics with the
finite difference picture.
OK.
So what's CFL?
C is Courant still.
So the this came from a really early paper, and I think
actually Lewy.
That's Courant, that's Friedrichs,
they were close friends.
That's Lewy, all three were friends actually.
And I think it was actually Lewy who spotted what we're
just going to do right now, the requirement that r had to
be less or equal 1.
Now let me just make clear that this
Courant-Friedrichs-Lewy condition is going to be a
necessary condition.
It has to hold or there is no hope.
But it's not enough, it's not a sufficient condition.
We can't guarantee that a method stable just because the
Courant condition is satisfied.
The Courant condition will tell us something about the
position of characteristics, we'll see it in a second.
For example, this centered one is unstable, this guy when the
right hand side you might say, oh this is a better method;
always better to take a center difference for du
dx than a one sided.
But not now.
The one sided one d OK.
The centered one is going to be totally unstable for all r,
this centered one.
And we'll identify that not from Courant-Friedrichs-Lewy
their paper wouldn't have spotted the difference between
one and two, but van Neumann did; van Neumann quickly
noticed that if you use an exponential, so we'll do one
van Neumann in a moment, if you plug in an exponential the
number 2, this center difference, it takes off.
And it's unstable even if delta t is small.
OK.
Now what is this CFL,
Courant-Friedrichs-Lewy, condition?
Where does that come from?
It comes from just thinking about the information.
I have to write down the thinking now that goes into
the Courant-Friedrichs-Lewy condition, the CFL condition.
OK.
It's straightforward.
I have to think at a time t, let's say t equals n delta t.
Suppose I've taken n times that, then the true u at,
let's say, 0 and that time t is, or x and t because we know
the true solution, the true u at x and t is as we know the
initial function x plus ct at the point x plus ct.
So that's the correct solution.
And Courant-Friedrichs-Lewy are just saying that you
better use this number in finding a difference method or
you don't have a chance.
Right.
That makes sense.
So what do I mean by use this number?
OK.
So let me draw.
So u at x and t is some point here.
This is at the point x and t.
This use these two values.
They use these three values.
This one use that one, that one, that one.
This one use those two.
Those four use these five values.
Right.
See the difference method is not propagating the
information entirely on that line, it's really taking
values over that whole interval with some
combinations of r and 1 minus r.
And in the end after five times step, it's giving us
this value up there.
But the point is that if this point x plus ct were out here,
then refining delta x an refining dealt t, and using
more and more point and filling in and using more of
these, but not getting out here, would
get you know where.
You wouldn't have a chance.
So that's the
Courant-Friedrichs-Lewy condition.
How far out is this?
This is n delta x inside this.
And if this is x, this is n delta x is away, and
this is x plus ct.
And c is n delta t's.
Right?
t is n delta t.
Are you with me?
As I keep delta t over delta x, the ratio r fixed and use
more and more points, I'm filling up this interval with
points that I've used.
But I'm never going to get to a point there.
So the idea is then unstable, couldn't
converge, couldn't work.
If the distance reached here n delta x is smaller than cn
delta t, it couldn't work.
OK.
So I'm taking a negative point of view here.
It can't work in this case.
The reason I take that negative point of view is I'm
not able to say it does work in the case when delta x is
big enough compared to delta t, and I include that point.
I can't be sure.
But I can be sure that if I don't include that
neighborhood of where the true u is produce, it
can't have a chance.
So let me just cancel n divide by delta x.
And I say fail if 1 is smaller then c delta t over delta x
which is r.
It fails if r is bigger than 1 by this reasoning of
characteristics.
There are two more cases that are worth thinking about.
Suppose the method was down wind.
So this is up wind, because I take a forward difference.
The wind is blowing this way.
Forgive me for incomplete use of the word wind.
I don't know what it means.
But one way or another whatever, the wind is
blowing this way.
The values are coming from the up wind direction.
And that's good, because the true value comes from the up
wind direction.
That's what we saw on characteristics, the values
are blowing down wind.
So what would happen if I use a backward difference here?
So I didn't write that bad idea, but it's important to
realizes that's method 1a.
A bad idea is to use a backward difference here.
I should maybe call it 1b for backward.
You see that it would fail?
If I used a backward difference, what will what
happen to this picture?
This value would come from stuff on the left.
It would totally have no chance to use the correct
initial value.
So a backward difference is an immediate failure.
It will also fail if --
well I was going to say r less than 0.
I don't know whether I can say that.
The way r could be less than 0 would be the delta x being
negative, and that's what I mean by the down wind message,
where I'm looking for looking for the value here, and of
course I'm not finding it.
OK.
So I guess I'm saying that Courant-Friedrichs-Lewy tells
us right away that this method will only be possible, And in
fact this is the stability condition.
So let me write that fact, but we haven't proved it yet.
It's stable for r between 0 and 1.
That's the stability condition, which
limits delta t.
Because r is delta t over delta x.
And this says that delta t can't be larger than the
multiple of delta x.
OK.
You might say that sounds like delta t is too small.
Is stability too restrictive?
Well it's restrictive certainly, because it
limits delta t.
But actually for these methods, accuracy also
limits delta t.
If I look for a method that allowed a bigger delta t, then
the error of that in that method proportional to delta t
would be too big.
In other words accuracy as well as stability is keeping
delta t of the order of delta x.
Accuracy as well as stability is keeping r within a bound,
and stability is keeping it within those bounds.
OK.
So I still have to show that the method is stable.
The Courant-Friedrichs-Lewy condition just helps to
eliminate impossible ratios.
OK.
One more point, which of course everybody notices in
solving the difference equation.
Suppose r is exactly 1?
Think about the case where r is exactly 1.
Remember that's c delta t over delta x.
Suppose r is exactly 1, then what does the method do?
The method takes the new value, r is 1.
This is 0 now.
So the new value is then when r is 1 is that old value.
In this picture, when r is 1 this value is this one.
The values are traveling along a line.
I'm not using these you see, because if r is exactly 1,
these are exactly 0.
Now is that good?
Actually it's terrific, because that line when r is 1
is the correct characteristic line.
If r is 1, if c delta t equals delta x, I'm going up exactly
with that slope that the true line does.
In other words, the difference equation is gives exactly the
solution to the differential equation.
Because it's so simple of course.
Couldn't do it.
Well that's the point, we can't do it for big problems,
because in a big problem c is going to be variable.
Maybe the problem we'll receive may depend on the
position or the time, or it might depend on the solution u
in a non-linear problem.
Those are ahead of us.
But there's possibility in 1D.
It's only a really a one dimensional possibility, where
the true information travels on a line, and we could just
hope to stay close to that characteristic line.
OK.
Now I'm ready to tackle for any method.
Let me get these guys up again.
So first I want to see why is up wind stable, and why is
centered unstable.
Should I start with the stability question rather than
the accuracy question?
Yeah, Christie we could guess what's first order until we
get to Lax-Wendroff.
So accuracy is going to come next time.
We've had Courants take on stability, the CFO condition,
but now I'm ready for van Neumann's deeper insight.
OK.
So how does von Neumann study stability?
He watches every e to the ikx.
So von Neumann's is going to be each e to the ikx should
have a growth factor in the difference equation
smaller then 1.
The growth factor in the differential equation of
course was right on.
This has magnitude exactly 1.
So wave equations are not giving us any
space to work in.
Because we want the growth factor and the difference
equation, we want it to be close to the real 1.
And at the same time it better be not the real 1 having
absolute value exactly 1 on the unit circle.
But the difference equation is allowed to go into the unit
circle, but not to go out.
OK.
That will separate the good ones from the bad.
OK.
Let me work here.
Can I copy the equation uj n plus 1 is equal to r uj plus
1n, and 1 minus r ujn.
OK.
Now I'm going to do an exponential again.
I'm going to start from e to the ikx.
Sorry, let me say times 0 and position n will be e to the
ikn delta x.
Right.
That's my initial, e to the ikx, that's my pure frequency.
I plug that in.
Oh not n but j; j is measuring how many, sorry let's get it
right; j is measuring how many delta x's I'm going across; n
is measuring how many delta t's I'm going up.
So this is the key, here's von Neumann.
Try a pure exponential, plug it in.
Again it's going to work, because we have constant
coefficients.
So what do I get on the right hand side?
I get r times this guy at level j plus 1, e to the ik j
plus 1 delta x times the exponential.
Can I save the exponential because it's
going to factor out.
Plus a 1 minus 4 times the exponential itself, and here's
the exponential: e to the ikj delta x.
This will be one step.
I could say uj1 coming from level 0 to level 1.
Sorry e to the ikj delta x, so I don't need it here.
Sorry I did it wrong.
That e to the ikj plus 1 delta x factors into e to the ik
delta x from the 1, and e to the ikj delta x,
which is over there.
So I only wanted this much.
OK.
Sorry that's on the permanent tape but maybe it's OK.
Because it makes us think through what's the result from
a pure exponential?
And again the result is pure exponential in, pure
exponential out, but multiply by some finite difference
growth factor g, that depends on the frequency.
And it depends delta x and it depends on r.
Well actually it depends on k delta x together.
So I could put that together.
What's von Neumann's question now?
He wants to know is this number, it's a complex a
number of course because it's got this cosine of k dealt x
plus i sign of k delta x, he wants to know is it magnitude.
Does the magnitude get bigger than 1?
if so n steps will give the nth power of the thing and it
will blow up, or does the magnitude stay lesser equal to
1, in which case it won't blow up, in which
case we have stability.
Of course it would be great if the magnitude was
always exactly 1.
Because that's what the true solution has.
But that's only going to happen in this special,
special case when r equals 1, and I'm going just right on
the characteristic.
When r is exactly 1, that's the case when that's gone,
this is a 1, and that has magnitude exactly 1.
But normally r is going to be, I'm going to be on the safe
side, r will be less than 1.
And I just want to draw this van Neumann
amplification factor.
I'll often use the word growth factor just because it's
shorter, but another word is amplification factor.
What does an exponential get amplified by?
Can you identify this number in the complex plane?
Of course the big question is where is it with respect to
the unit circle?
OK.
Take k equals 0 --
0 frequency the DC term, the constant start -- what is that
number equal when k is 0?
1.
So that's normal that k equals 0.
We're right at the position 1.
That's just telling us that a constant
initial value stays unchanged.
The true growth factor and the van Neumann amplification
factor both 1.
But now let k be non-0.
Where is this complex number?
And now I'm going to impose the condition that r is
between 0 and 1, because I know from Courant that is
going to be required.
What's the point here?
OK.
Then you see here is a 1 minus r.
So there's 1, 1 minus r will be a little bit in here, and
then what happens when I add on this number?
So the 1 minus r I've done.
It put me here.
It was real.
But this number is not real.
This just number is r times complex number,
but what do I know?
I do something critical here.
What do I know about this number?
I know its absolute value is 1.
So that it gets multiplied by an r.
Look, look, look, look.
I start with a 1 minus r.
So I go to 1 back to r, and now I add in this part, which
is somewhere on a circle, a radius r that's everything in
this problem.
It's a circle of radius r around the point 1 minus r.
I'm inside the circle, and you could say well you could have
seen that clearly by just using the triangle inequality.
The magnitude of this can't be bigger than the magnitude of
that, which is what?
r.
Plus the magnitude of that, which is what?
1 minus r.
So the magnitude by the triangle inequality can't be
more than r plus 1 minus r which is 1.
But now wait a minute, where did I
use the Courant condition?
The way I said that there sounded like it
would always work.
This has magnitude r.
What's going on here?
Suppose r is bigger than 1, what's going wrong?
In fact, what's the picture if r is bigger than 1?
If r is bigger than 1, then my 1 minus
r is out here somewhere.
So this is the unstable case, 1 minus r, and then a circle
of radio r, bad new right?
Every frequency is unstable in fact.
When I'm reaching too far, and that's just telling me I'm not
using the information I need.
I don't have a chance of keeping
controlling any frequency.
OK.
I've got just enough time to show how bad this second
method it is.
So nobody wanted his name associated
with that second method.
Why is that?
So what's the van Neumann quantity
for the second method?
OK.
Just space here to write.
I want to know the number for that second method.
So now the method is ujn plus 1 is ujn from the time
difference plus r over 2, uj plus 1n minus uj minus 1n.
This is going to be the bad method.
It looks good because the center difference is more
accurate than a one sided difference, but it's also
unstable here.
Can you look at that and see what the van Neumann number g
is going to be?
Think of an exponential going in.
What exponential comes out?
So g is a 1 from this.
Now this is our guys that's shifted over, so that would be
in r over 2 e to the ik delta x just as it was before.
And this one will be a minus r over 2 times e to the -- so
it's back one step, so there will be and e to the
minus ik delta x.
So what's with this now?
g is 1 plus this quantity, r over 2 times e to the ikx
minus d the minus ikx.
What am I seeing there?
I'm seeing 1 plus r, and everybody recognizes this
minus this over 2 as i sin, ir sin k delta x.
That's the amplification factor, and is it
smaller than 1?
No way.
Let me raise it up here.
So I don't care whether r is less than 1 or not, I'm lost.
This is a pure imaginary number.
This is a real number.
So I have the sum of squares to get the magnitude, it will
be the square root of 1 plus r sin k delta x squared.
So if I draw the bad picture then, what's the bad picture
is, now it goes to 1 and then it goes way up
the imaginary axis.
Well maybe not way up but up.
Sorry I shouldn't have made it quite as bad as it was.
If I reduce r, I don't go so far up, but no hope.
So do you see why that one is bad, because the amplification
factor van Neumann tells us what Courant did not tell us.
that the amplification factor here has
magnitude bigger than 1.
It's outside the circle, and every exponential is growing
and there's no hope.
OK.
So the second last lecture on this section 52 of the notes
will be about Lax-Friedrichs and Lax-Wendroff order of
accuracy, stability and actual behavior in practice.
Ok.
See you Wednesday.
Thanks.