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 a 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 m c 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
m c 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 r 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 a 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 R*P 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 the 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 a while --
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 xt-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 a 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 u_0 of x.
So let me graph. u_0 of x is -- the key points are 0 and 1,
0 and 1.
So u_0 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 xt-plane
where here I'm in the xu-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, the 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, at t equal to half, this state 1,
has penetrated this far, 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 and 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 xt-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 this 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 interval
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,
Burgers' equation because it's a 0.
So just one word about Burger's equation.
So Burger allowed a u_xx term here.
So he had one of our -- nonlinear convection
and a diffusion term, and his convection was this simple
expression u*u_x.
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 as 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 v_max or 2*v_max.
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 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, the 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 u_0 to be 0 now, u_0 to be -- well,
I can even -- it's called a Riemann problem when there are
just two values here.
u_0 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.