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 you see that I thought maybe we've got enough
good examples in theory where we need some practice.
We need some experiments with some of the finite difference
methods that we've spoken about.
So I created this homework in which you would like choose
one of these that we've done.
Well, actually, number two, the Schrodinger equation,
which has that imaginary i in there we haven't spoken about
yet, but you could do the kind of Fourier analysis that we've
done for the heat equation.
But, of course, it's going to be different
because of that i.
You'll get an e to the minus i k squared t, I guess.
And that means that's a number of absolute value 1 instead of
a rapidly decaying coefficient.
OK, so my idea is to choose one of those groups, or number
five, if you want.
And say, by next Friday, just give me some
experiment with a code.
In this case, you could use -- well, we have four methods and
three of them are already in code, and you're welcome to
use those or your own.
And I raised questions last time -- and I'll put this on
the website with a little more detail.
So just to give you an idea of what's coming.
So this is the problem that we're in the middle of
speaking about, near the end, the mixture of
convection and diffusion.
And this is Schrodinger, which is new.
This is the advection equation, or the one-way wave
equation, which we've started with, and this is what's --
this'll be the subject of today's lecture and next time.
And the best reference I can give, as well as the textbook,
the applied math book, has significant discussion of this
conservation law.
But also it's very well presented in the 16.920
lecture notes that are on OpenCourseWare,
lectures 11 and 12.
So I'm going to use those as a guide in this lecture.
Anyway, I thought I'd just put that up as, you know, time to
see what various difference methods will do, And again,
it'll be on the website.
OK.
Now to say another word about this one, this convection
diffusion equation.
So last time I wrote down a difference equation very much
like that except what I wrote down last time was upwinded so
that was just j.
Now it's j minus 1.
And the choice there is not clear.
This gives us higher accuracy, of course, for that x
difference because it's centered at j.
And we still have this centered, of course.
The second difference is centered.
So I wrote down here what the coefficients are of the three
old values at j plus 1, at j, and at j minus 1.
But there is one important point here.
The price for that extra accuracy may possibly be a
negative coefficient.
This could be negative, or it could be positive.
It depends which is more important: the small r, the
convection part, or the big R.
And let me follow up.
So you remember the price.
If we see a negative coefficient, then we expect
some oscillations to appear.
And actually, they were supposed to appear in Figure
5.12, but it hasn't been done yet.
And if you would choose this one and do Figure 5.12, I'll
be incredibly grateful.
So Figure 5.12 was to show the evolution from a u naught for
the convection diffusion.
OK, so can I just -- do you see that in comparing R with r
over 2, those are
dimensionless quantities, right?
Can I take another minute to try?
I'm a total amateur at dimensional analysis.
I did always think, however, that Einstein
could have proved --
well, not maybe proved, but could have conjectured that e
was mc squared, the most important equation in
relativity, just by checking the dimensions.
I don't know if that would --
I mean, right?
That would like -- might have occurred to him to realize
that mc squared has the dimensions, if we wrote those
out, of energy.
I don't know if that put an idea into the system which
produced, of course, the fact that -- and, of course, there
would be a dimensionless constant, which amazingly
happens to be 1.
OK, so dimensional analysis could like change the world.
This won't quite change the world, but a
little touch of something.
OK, so can I just --
you remember that this Peclet number, and I'll write his
name down again.
P-E-C-L-E-T. The Peclet number.
In the differential equation, the Peclet number was a ratio,
essentially a ratio, of c to d.
But to get the --
I think we needed --
was that what we needed?
We needed a length to make it dimensionless.
c is the coefficient for --
c over d sort of tells us roughly the importance of
convection and diffusion, but there's a
length scale to be included.
Now the cell Peclet number, which I'll just call P, I'm
just going to take to be r over R because those are
dimensionless.
And you see that it's going to play a part here.
Let me see, I forgot.
I think I want r over R. Maybe I want a 2 just because I
would like that to --
I'd like it to come out simply there.
OK, so I may put in a 2.
So that are would be 2R p.
Where did that 2 come from, by the way?
Oh, I guess it came from there, right?
It came from that minus 1/2.
So probably I want that 2 there just so that I have a
very simple decision of whether this coefficient is
positive or negative, yeah.
So let me put in a 2.
OK, so r over 2 is RP then.
This quantity is R times P, so this is R minus R times P, or
it's R, which we know is OK, times 1 minus P. In other
words, P smaller than 1, that coefficient --
all the three coefficients are positive.
It's certainly stable.
Not only stable, it can't oscillate.
We can't have -- we're taking a combination of the three old
values and oscillations can't show up.
OK, so just to see, what is this?
And you remember that r is c delta t over delta x, and
capital R is d delta t over delta x squared.
So just to complete this, the delta t's cancel.
One of the delta x cancels.
The other one comes up.
So it's c over d, and then we have a delta x and a 2.
So that's why it's called this cell Peclet number because
it's the cell size, or in my convention, half the cell size
that gives the length scale and decides whether P is -- so
I'm just going to write down the conclusion that the figure
should show.
P smaller than 1.
That means that this coefficient is also positive.
So that's positive coefficients, no oscillation.
And I think you'll see from the output that P bigger than
1 produces oscillation we don't like.
In other words, we really don't want P larger than 1.
And of course you could say, well, look at P. It's got a
delta x in it.
That's not going to be a big number.
And the method will work fine as delta x and delta t go to
zero, but we're going to do a finite calculation.
And so we're seeing for the first time, like, because we
have two terms that are not dimensionally identical, c and
d, that in that competition between convection and
diffusion, the length scale is crucial, and the
mesh size is crucial.
So I guess the answer is we can resolve it by taking delta
x small enough, and hopefully, that's not a problem.
But, of course, you know that there are many physical cases
in which the diffusion coefficient or the viscosity
coefficient is small, is quite small.
And so getting that below 1 might be hard because this d
could be a very small number.
OK, so maybe can I also mention on the -- as something
that we haven't done here, was to move the second derivative,
this guy, to time n plus 1, make it implicit.
In the diffusion part, one often does that.
The diffusion is the term, is the one that's forcing us to
take this coefficient R below a half, delta t of order delta
x squared, and the way to avoid that is to move the
viscosity term, at least half of it, up to the new time.
OK, so that's a third method.
So we've done upwind last time, centered, explicit this
time, but then that's a third one that has to be looked at.
OK, so that's where we stand on linear problems.
I'm sort of ready to jump into nonlinear ones, recognizing
that -- well, there's always more to say about linear ones.
What are some of the things that we have
not properly discussed?
Boundary conditions, above all.
Boundary conditions, and we could do something with those,
and the notes will say more about boundary conditions.
So there are several types of boundary conditions.
I guess I'm going to say something about them now, but
this'll be a board that kind of gets covered up.
So boundary conditions: So what are the different types?
Well, at a real physical boundary, we might have u
equals 0, say, in heat flow.
We might hold the boundary.
It might be the edge of a freezer or something.
We might hold the temperature at zero or
at some given value.
We might prescribe u at a boundary, or we might
prescribe the derivative.
And just in case we're in several dimensions, I'd better
write it as the outgoing derivative.
It could be prescribed, say, 0.
So what are those two types of boundary conditions called?
One's called absorbing, an absorbing boundary, and the
other is an insulated boundary where no heat goes through.
To maintain this, heat probably does go through.
If our body is warm and we're holding it at zero on its
boundary, then heat is going to travel out.
And in fact, as time goes to infinity, we're going to
approach cold -- a uniformly cold.
Here, heat's not allowed to travel out, so what will
happen as time goes on is it will approach a constant.
The temperature will approach a constant because it can't
leave and it diffuses.
And now I wanted to mention a third type of boundary, which
is a purely computational boundary.
We talked last time about Maxwell's equations in free
space, and the big application of Maxwell's equations, or one
of the big ones historically, has been to -- like the
exterior of an airplane, to radar, or other
electromagnetic signals.
So what do we do if we have an exterior problem?
See, here, I was always thinking
about an interior problem.
You know, we're inside 0 to L or something.
But now, what do we do if a region is the outside, it's
free space all the way?
Well, we have to create a
computational box of some sort.
So we create a box in which we compute, but the boundaries of
this box, the x-boundaries of it, are purely artificial.
They're just because we can't compute to infinity.
So waves travel in this box, and they travel to the -- so
you see what I mean?
We're solving it in all this space, except that's too much,
so that's meant to be a very large box, not a small one, as
large as we can afford.
But waves are still going to hit the boundary of that box.
OK, what do we choose as boundary condition there?
We don't want it to reflect back because
it's not a real boundary.
So we have to create --
I'll put ABC for absorbing boundary condition.
Absorbing.
That A is for an absorbing boundary condition.
And the creation of these boundary conditions is an
important topic within finite differences.
For Maxwell's equations, a good idea was produced and
keeps being developed, called perfectly matched layer.
The idea in that perfectly matched layer -- let me just
write PML -- and I haven't seen this so much in other
application areas.
The idea of this perfectly matched layer is to create an
artificial dielectric, I guess, a dielectric medium it
would be in the electromagnetic case.
An artificial thin layer, which has just the right
properties to stop the wave and swallow it.
OK, I'll put on the website references to that.
It's a bright idea, but not totally simple, and when we
change difference methods, then that perfectly matched
layer has to go with it.
OK, that's the end of my thoughts about linear
problems for now.
We're coming back, of course, in a few weeks to solving
large systems, large linear systems, the kind that we
would meet if we use an implicit method, or the kind
that we meet if we're solving a steady state problem.
But for now, let me go to that model equation.
And you see right away the difference between that
equation and the advection equation that
had a constant velocity.
Here we have of a velocity that -- so c
is minus u, I guess.
c is minus -- think of c -- this u as being
comparable to minus c.
So, in other words, if u is positive, I'm now going to
have waves going to the right.
So let me just put this down. u greater than 0 will mean
waves going to increasing x to the right.
And the idea will be -- so we have to discuss that problem.
And as always, we would like to solve it analytically.
Find whatever formulas we can, understand what's happening in
terms of these characteristic lines because when we have one
equation and one space variable, we really can catch
the essence of the solution by understanding the
characteristic lines.
And then how do we do it numerically?
OK, before class, I wrote just a few equivalent forms of it.
This is an equivalent form, and this is, so to speak, the
right form.
This is the form where we see a time derivative of u, as
always, and we see a space derivative of u squared over
2, which, of course, produces that.
But somehow this is the quantity, it's called the
flux, that's physically meaningful, f of u.
So I've chosen a simple flux function, u squared over 2.
So that's the differential equation.
But you'll see that the differential equation can
produce two solutions at the same point,
and we have to choose.
In other words, the differential --
the solution can become discontinuous.
A perfectly smooth starting function, after awhile -- this
is the essence of the problem now.
I might have a smooth starting function, at least continuous.
The solution is constant along characteristic lines.
But what happens if two characteristic lines run into
each other?
I've drawn that possibility here.
And it'll take a little time to see exactly what's
happening, but the
characteristics we can see here.
So the starting values were 1 up to that point.
So when the starting value is 1 and u is constant along
characteristic lines, then u is 1, and those lines, it's
just like having c equal to minus 1, I guess, in the
one-way wave equation.
We have a one-way wave going to the right now along
characteristic lines.
So u is 1 in this part of the x, t plane, right?
We have to come back to the formulas to
confirm what I just said.
And then over here on the far right, the starting value is
0, and the characteristics in that case, the c associated
with that is 0, so it's standing.
It's not moving.
The characteristics are just straight
lines in the time direction.
Nothing's happening, right?
When you think of a characteristic that's going
straight up, it's just carrying the value of u
without moving left or right, so u stays 0.
So u will be 1 here, u will be 0 here, and now something's
happening in here, and we have to figure out what it is.
But maybe I can anticipate the main point before any algebra.
The characteristic lines in this region have a slope
somewhere between the slope 1 and the slope 0.
You see, we're starting -- let me graph u0 of x.
So let me graph. u0 of x is --
the key points are 0 and 1, 0 and 1.
So u0 of x, I start with a constant initial function, a
constant profile to the left of 0, and a constant profile
to the right of 1, and linear between.
OK, so it's that starting values that I wrote here: 1, 1
minus x, and 0 are these three pieces.
And then here I'm in the x, t plane where here
I'm in the x, u plane.
So I drew characteristics, and we still have to track.
We have to say more what characteristics are and
understand that these have slopes in between, but you're
not surprised because the values of u are in
between 1 and 0.
So they go that way and they meet here.
So up to time t equal 1, characteristics tell
everything.
Up to t equal 1, I know that u will be 1 along these
characteristics, and when I find these characteristics, u
is whatever it starts at.
And u is 0 up all of these characteristics, so up to that
time, algebra --
the differential equation is going to be totally OK.
But after that time, the differential equation is in
difficulty.
Because at that moment t -- so what's happening as t
increases from 0 to 1?
This is the picture at t equals 0.
The picture at t equal to 1 --
t equal to half, let's say.
What would be the picture at t equal to half?
Well, this t equal to half, this state 1, has penetrated
this far to halfway.
Then in here is this converging characteristic that
come from a linear profile.
And over here, it's all zeros.
So that zero is still there, but you see what's happening.
This guy is -- the profile is steepening.
This is at t equal to half.
Steeper.
We still have the algebra to do, but I think the first
thing is to see roughly what's happening.
Then at t equal to 1, this is totally steep.
So this is the situation at t equal to 1.
And I could've given you that as the initial condition, but
I, in the 16.920 notes, thought, OK, let's have a
little peace for awhile, see what those characteristics are
doing when they're converging, but they're not crossing.
The trouble is, after t equal to 1, they are crossing.
And now, what's the solution in this crosshatched region?
Do I take the value of 0 that's coming up along these
characteristics?
Do I take the value 1 that's coming along these
characteristics?
What's happening in here?
So there's a crosshatched region there, which is a big
question mark because my rule of following characteristics
has led me to two answers.
And we want a physically relevant answer.
We're sort of expecting, and we'll see, that this step
function is a step function.
I mean, the answer will be a step function.
It will be 1 and then drop to 0 along some path in the x, t
plane, which is not a characteristic.
The characteristics are going at 45
degrees and at 90 degrees.
And actually, so this is the region I'm interested in.
So can I blow up that region?
So if I blow up that region, that's where the shocks -- the
shocks are meeting everywhere in this region.
This was the point.
What was that point?
I guess that was the point x equals 0 and
t equal to 1 there.
OK, and the question is -- and I know that I have u equals 0
over here, and I know -- oops, u equal 1 over here, and I
know that I have u equal 0 here, but the question is
where and -- well, what happens in between?
You can think of at least two possible scenarios, and one of
those is that there's a shock line that goes, let's say,
somewhere in the middle.
And so u stays 1 all the way to that line and
drops to 0 after it.
In other words, this wall of water travels to the right.
The shock, the discontinuity, moves along, not necessarily
at the speed that this one.
This would be moving along the speed of 1 with it pushing.
Not stationary either, which is what this right side would
like, but somewhere between.
So that's a shock line.
I'll just, for short, call it a shock.
Then there's another possibility.
And that will happen in this case, actually.
That will happen in this case.
But there is another possibility that will happen
in other cases, and that would be a fan of characteristics.
Let me see that when it comes, but I'm just mentioning it to
say that this business of the shock line, which is what
we'll pursue now, which because it's what happens with
this problem, is not the only --
is not always the way to match, the way to deal with a
problem of lack of information coming from the
characteristics.
OK, so I guess our question is, what is that shock line?
What's the path?
What's the slope, and what controls it?
And the answer will be that what controls --
we run into problems with the differential equations.
And the way to see what continues to happen is to go
to the integrated form, the integral form of the equation,
where, of course, by taking an integral, I make everything
smoother, and also I preserve the physical law, the
conservation law.
So this is really the conservation law.
What is this?
Well, let's first of all see where it comes from.
It just comes from integrating this equation
with respect to x.
So I have a time derivative that I keep out here, but my
space derivative is going to disappear because I'm
integrating with respect to x.
So I can call d by dt or partial d by dt the integral
of u with respect to x, and I'm integrating from somewhere
-- from some point to some other point.
Instead of a and b, I'll call them --
I really should --
I'm sorry, I should say x on the left and x on the right.
Some point x at the left and some point x at the right.
And now I've integrated this one.
I've called -- this is f of u so I'm now thinking of any
flux function.
And the integral, of course, gives me f of u at the upper
end minus f of u at the lower end.
OK, so that's the conservation law.
Let me write that word down.
This is the conservation law.
What is it saying?
It's saying that the change in the quantity in the integral
is given by -- if there's any change, any time derivative,
the rate of change in the interval --
it's the conservation of this integral meaning that any
change in that integral happens because there is flow
going out the right-hand side of the integral or there's
flow coming in the left side of the interval.
So that's why we have a minus sign between them.
So that's a statement of conservation, and --
oh, there's one more example that I'll mention now.
Let me give this particular equation a name, just because
you see it.
You often see it written -- called Burger's equation.
But actually, the full name is the inviscid, no viscosity,
Burger's equation because it's a 0.
So just one word about Burger's equation.
So Burger allowed a uxx term here.
So he had one of our -- nonlinear convection and a
diffusion term, and his convection was this simple
expression uux.
It's the kind of thing you see in the Navier-Stokes
equations, right?
The nonlinearity in Navier-Stokes is sort of in
this term of this general sort, but this is a major
simplification.
Such a simplification that Burger could solve, even
though the equation was nonlinear.
Well, two people at the same time discovered the change of
variable, a simple little trick.
Hopf and Cole: I'll put their names down because they had
nice, short last names.
They both discovered --
in the case of a diffusion term, they both discovered a
little -- a change of variables that made the
equation linear.
Then they could solve the equation.
That's in the textbook in the section on conservation laws.
And then the important thing that Burger could do was let
the diffusion coefficient go to zero.
So then in the limit, he got the inviscid Burger's
equation, the one we're looking at.
And he could look at the solution and let
the limit go to zero.
And that's another way to figure out
what's the right solution.
That's the viscosity method, a fundamental method, for
finding out when we run into a problem, shocks
run into each other.
We don't know what to choose.
One way to find out in this 1D problem and also in much
tougher nonlinear multidimensional problems is
put a little viscosity in, and let it be small.
Let it even get smaller.
Then with a little viscosity, officially this term -- in
some sense, this term, because it's a second derivative, is
officially, mathematically speaking, sort of dominates
this one as far as --
I mean, this will produce a smooth solution.
This term, as we know doesn't smooth the solution.
It just advects it, just carries it along.
But this term will smooth it.
So this is the dominant term, the second derivative.
No surprise.
So when we introduce that, that smooths everything.
We have no problem to say there is a solution, at least
for simpler models, and then we can see what happens when
that term approaches zero.
OK, so that's another method.
And, of course, that's kind of what we're doing with finite
differences.
Finite differences, if we put in a little viscosity, a
little second-space difference, and take delta t
small enough for stability, then we can solve it.
And if we take it very small, we should get something close
to the inviscid limit.
OK, I was just giving Burger's name to
that particular example.
So this isn't the only example, but it's
the simplest example.
OK, another example of conservation laws that
everybody likes, I think, is the traffic flow, traffic flow
in one direction, say, on the Mass.
Pike.
So what's the equation of traffic flow?
Let me put that on a board here, and then I have to come
back to this.
Well, I'll put it here.
So Example 2 -- well, I don't know if you can call
that Example 1.
We didn't give it a physical meaning.
Example 2 will be traffic flow in 1D, and what's the unknown?
Well, it's the thing I've been calling u, but for traffic
flow, the natural unknown is the density of the traffic,
the density of the cars.
And everybody uses rho for density.
So can I use rho rather than u?
And then in the equation will be a car -- a velocity that
I'll call v. That's not a new unknown.
In fact, the velocity experimentally depends -- oh,
there's a relation.
The density and velocity are, by experiments, if you run an
experiment on the Mass.
Pike in heavy traffic -- that could be a homework problem,
too -- to find the relation between v and rho.
So I'm going to keep this rho was the unknown and imagine
that we have some experimental velocity v of rho.
And what would we expect?
I guess we expect that as rho increases -- as the density
increases, what will happen?
As the density increases, the traffic is going to slow down,
ultimately come to a stop, where if the density is very
small, so the density 0 --
I mean, we just have a single driver, the police are not
out, so he's going vmax or 2vmax.
But as the road gets crowded, the velocity drops and maybe
drops linearly.
That would be a possible and not that terrible relation
between v and rho.
So what would be the flux function in this
traffic flow problem?
The conservation law is like conservation of cars, right?
If we look at a unit, a piece of the turnpike, then the
conservation law says that the total number of cars inside,
which is this, changes by cars leaving and by cars entering.
And so the flux function I think will be -- so we have to
know like sort of how many cars are leaving, and I think
the flux function would be probably v times rho.
I think the flux function of rho now is my unknown.
I think it would be v of rho times rho, the velocity times
the density.
That tells me how many cars are going out.
OK, so I guess what I want to say is that we will meet --
that we meet in this traffic flow example exactly these
possibilities that -- we meet them mathematically just the
way we meet them in actual driving -- of a shock.
So this shock corresponds to -- so what that's?
Like stop-go --
this corresponds to I suppose like meeting a red light?
Coming up to a red light, traffic you're
behind normally stops.
It has to slow way down and stop.
Then the light turns green.
The first car takes off.
The cars begin to spread out.
That's the solution that I don't see with this starting
function, a sort of fan of cars.
So maybe I could write those words down: shock or a fan.
The shock is what happens when characteristics come together,
cars come together.
The fan is when the characteristics fan out.
Actually, I could illustrate a fan here right away.
Suppose I reverse 0 and 1 just to see what --
OK, so here's my -- this is my x again.
So here I'm going to have u0 to be 0 now, u0 to be --
well, I can even -- it's called a Riemann problem when
there are just two values here.
u0 is 1.
So the so-called Riemann problem is the cleanest
example of the whole conservation law theory when
there are just two starting values 0 and 1 or 1 and 0.
And it makes a big difference.
1 and 0, we saw when we reached time 1, we had a 1 and
a 0, and after that, a shock formed.
Here 0 and 1, the characteristics here
are going as always.
So u is 1 here.
The characteristics here, this is staying at 0.
There's no reason to change.
But here we have now a region, just like the other one, but
with a major difference.
That now, it's 0 on the left and 1 on the right, where
before, it was 1 on the left and 0 on the right.
And the question is, what happens in here?
Well, you might think a shock, a step up from 0
to 1 at some point.
And I would have to agree that you could do that in such a
way that you maintained the conservation law.
But that's not what happens.
This is going to be a fan in here.
We're going to have -- what we'll have is
characteristics in here.
This will --
let me draw the profile.
So the profile is 0 there to 0 at some time.
Then over here, where this characteristic is met, it's 1,
and this is what happens: Instead of getting worse, less
smooth, the solution gets better, more smooth, when the
direction is reversed.
So we see that our problem had to be nonlinear.
A linear problem couldn't do these two opposite things in
the same thing.
So here we would have a fan and the characteristics go
there, and it's reflected there.
OK, so there's a picture of the phenomena that we're
looking at and we would like to see numerically, too.
So next time I want to find -- so I can say now what I have
to do next time.
One is to locate the shock, the shock line, when there is
one, and two, to decide -- how do we decide between shock or
fan when from the point of view of pure conservation, we
could choose.
But nature makes a choice.
And one way to decide, make that decision, is Burger's way
of putting in a little diffusion and letting it
disappear, but we need a direct rule.
So we need a direct entropy condition.
So this is going to be a shock speed -- a formula for shock
speed s. s equals shock speed we'll find
if there is a shock.
And number two is going to be -- this decision is going to
be based on a quantity called entropy, which appears only in
nonlinear problems, really.
OK, so that's a beginning picture, which these 16.920
notes are an excellent reference for.
The book also and Monday's lecture will pick up there.
So we can see the two situations here, and the
question is what do we do at this point.
OK, see you Monday.
Have a good weekend.
And any email questions about the homework, just you could
send to me or -- and/or to -- maybe and to Mr. Cho, who will
be looking at the homeworks that come in.
Good, thanks.