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: All right, so I'll start with this review board,
which shows the conservation law.
Scalar, one unknown, one dimension, x, one space
variable, and the flux function is f of u.
That's the differential form.
We have seen that at some points that differential form
produces a discontinuous solution, a discontinuity.
In that case we are better off with the
integrated intgeral form.
So I just integrate that one from a point x left to a point
x right from an a to b and that gives me--
I have then the integral.
So if u was density for example, the integral of
density will be the mass.
That's the quantity that I'm seeing here, so this is the
rate of change of mass in the interval, comes from flux
going out minus flux coming in.
So that's the conservation law in its pure statement.
And we learned that solution for this-- while we're looking
at the differential problem-- we can follow these
characteristics from the initial time zero.
So at time zero we look at a point x and the value at that
point just travels--
is sort of safely carried up the characteristic line, which
is that equals x0.
And the problem arose when characteristics either
collided into each other or separated.
Those are the two difficulties and they have different
resolutions.
So this is what we met in the traffic example when
a light goes red.
When a light goes red cars are coming up to the light, in
front of the light the density is low.
The cars have to brake quickly.
So something happens-- that's a question mark there.
What do you do when characteristics collide?
When the value here cannot be both the initial value that's
coming up that characteristic and the initial value that's
coming up there because they're different?
So we need a rule for what to do there.
And we also need a rule for what to do in
case of an empty space.
When no characteristics--
here are two characteristics are hitting, here are none are
touching there.
Characteristics are separating in this case and the question
is, what solution do you fill in the middle when you haven't
got a characteristic?
It's not coming from the start.
OK, what solution to choose?
The answers will be different.
And let me take the first one, the collision question, first.
So what will happen is a shock will form.
There'll be a shock line so that these characteristics
carry information into the shocked, by the way,
that's a big deal--
until they hit the shock and these characteristics carry
their information.
And across the shock we have different values, u on the
left and u on the right.
So that's one possible form of the solution because that will
solve the equation in this region where the
characteristics from the right are following the formula.
It'll solve it from the left, so everything comes down to
the jump condition at the shock.
What's the shock speed?
Where does this shock go?
It's to find that shock in the xt-plane that I have to use
this integral form.
I must use the integral form because I've got a
discontinuity and the differential form just doesn't
make sense.
OK, I wrote down the key idea here or the conclusion, but
not the reasoning.
I'm looking for this jump condition at the shock, which
is going to locate the shock.
And the key quantity that I have to find is the shock
speed, which will tell me where that shock is going, how
it's moving up in time.
So time is going that way, x is going
this way in this picture.
So this has a speed, s of t.
And here's the conclusion.
There's the jump condition.
So what does it mean?
This notation, this bracket notation is the jump.
The shock equation is beautiful.
It says that s, the shock speed, times the jump in u is
the jump in f of u.
And let me take this example that is sort of the Burgers'
equation example, so this is the example that comes from
Burgers' equation.
That's the example where the flux is this nice parabolic
function u squared over 2.
Then the change in flux, the jump in the flux function is u
right squared minus u left squared over 2 and I divide by
the jump in u, which is u right minus u left and that
cancels beautifully and leaves me with the 1/2 of u right
plus u left.
It tells me that the shock speed is exactly halfway
between, it's exactly the average in this Burger's
example of the characteristic speed from the right than the
characteristics speed from the left.
And notice, it's in between.
That's crucial.
That tells us that this thing here is good.
It tells us that characteristics are entering
the shock from both sides.
Characteristic lines are entering the shock.
The shock is somehow in this sort of nonlinear case, it's
somehow kind of cutting out bits of the solution, reducing
the energy, reducing the entropy, reducing the total
variation as we'll see it in a minute.
So that's the answer here and I guess I have to justify,
where did this come from?
Well, it came from this integral form.
Suppose I integrate from a little bit to the left of the
shock to a little bit to the right of the shock.
Of course, that's where I want to integrate in order to
capture the shock.
OK, so if I just take that equation and write it over
here and I'll imagine-- because I'm not going to be x
left and x right are going to be pretty near the shock--
u up to the shock will pretty much have a
constant value, u left.
So approximately that is the shock position.
Let me say that shock is at captial X is the
position of the shock.
So I have x minus x left times u on the left of the shock.
I'm doing that integral.
Oh, I need d by dt of all this.
So d by dt.
It's d by dt of the integral of u dx.
Let me just put here what I'm computing.
I'm computing that from x left to x right.
Well, I'm going up to x, up to capital X where the shock is
and then onward.
So here I went up to x using the value u left.
And now will be the onward part, which will be x right
minus x, the half of the interval that's on the right
of the shock times u right.
And then copying this term is exactly our jump in f of u.
Plus the jump in f of u is zero.
Those words, total variation are a reminder to myself for
what comes next.
OK, so is that OK?
That will be approximately zero because u wasn't exactly
constant on the left and right.
It didn't stay at the value ul, but we're squeezing the
interval down where it's essentially constant.
OK, so now just look and see what we've got.
dx dt is the shock speed, so that's s.
Multiplying ul and here I have the derivative of that is
zero, the derivative of xr is zero.
Here I have a minus, so it's minus ur.
dx dt times that number is e plus the jump in f of u, left
to right, is 0.
And that's our equation.
If I put this on the opposite side it becomes u
right minus u left.
Shall I do that?
So I'll keep this on this side, got a good sign, but
I'll move this onto the other side so it will be the shock
speed, d capital X dt times ur minus ul, which
is the jump in u.
So that's great.
We now have a way to deal with shocks that collide.
So that was a derivation of this jump condition, which
first had a name, two authors [? Rankin Whogonio ?]
in gas dynamics, but then it was seen as the right
condition in all conservation laws.
OK and now I'm ready to comment on the second issue.
What do you do when characteristics separate?
Well, let me say what you don't do.
You could, without wrecking the equation, you could try to
put a shock in there and give this the value, u on the left
even though a characteristic isn't telling us what's going
on in this space or this one.
You could try to put a shock in there.
You could maybe satisfy the jump condition if you put the
right shock in, but the characteristics would be
coming out and that's wrong.
And it's called the entropy condition, which says that
this can't happen.
That the derivative of entropy has a certain sign and it
rules this out.
It rules out this possibility of characteristics emerging
from the shock and we need a totally different solution.
And it turns out to be a fan in there.
So we put in a fan of characteristics coming from
here and that gives a solution that actually
will depend on x/t.
This solution will be a simple ratio, x/t enters that.
In other words, that gives us a smooth solution and I guess
thinking about the traffic example, that's
exactly what happens.
So this is the shock part, now for the fan part, what happens
when we have a bunch of cars sitting there, ready to go,
the light goes green, the first cars start moving, cars
behind them see the car in front moving, they start
moving and so the initial profile would be, in this
case, for cars, would be high density and then zero density.
So this would be maximum density sitting at the light
and zero density in front of the light and then the light
goes green now and what happens?
Well, as we all know those first cars pull ahead and the
profile is these cars haven't yet started to move.
These have started to move and there's a profile there and
then zero density.
So this represents the first car, the lead car.
In front of that car is zero density and that's the fan.
That's the x/t expression.
So a sample linear expression, which will flatten out as time
goes on, these cars will be going faster and faster and
then eventually up to full speed.
OK, so that's the fan.
Maybe now is my chance to speak about this concept of
the variation in a function.
The idea will be you now, how do we find the difference
method to do this automatically?
Well, OK, one way to do it is put in a little viscosity.
This rule, the rule that picks out this shock in the one case
of colliding characteristics and this fan in the case of
separating characteristics, one way to get to those two
solution would be to put in a little bit of viscosity and
then let epsilon go to zero.
And the result would be these two solutions.
And in a way that's what finite
differences is going to do.
Finite difference or finite volume is going to be a
numerical viscosity and the epsilon in the Burger's
equation would become delta t or a delta x in the numerical
method in the finite difference equation.
And then what's the concept?
So we have to find finite differences.
And we want those--
the point is finite differences, sort
of takes more care.
We can't just create something that's consistent in the
smooth case.
Up till now we've just had no hesitation.
We replaced derivatives by differences without too many
worries, really.
And that was fine as long as the solution is smooth.
But now if we're going to replace derivatives by
differences we have to do it in such a way that this
entropy condition, these solutions emerge
automatically.
And really, in a way that our finite difference or our
finite volume methods should somehow conserve mass the same
way the true equation does.
So those are two different points.
First point is this question of total variation.
What's the total variation in a function?
Well, let me draw a function.
The total variation is-- well, let me see.
Total variation in a function u is going to be, by
definition, the integral of the absolute value of du/dx.
OK now we have to understand this concept.
It's beautiful that the analysis has led to a quantity
as simple as that, as a guide.
So when we create finite difference methods we will
ask, are they total variation diminishing?
So TVD will be total variation diminishing methods.
Methods where this quantity does not increase in time.
Let me just understand this quantity first.
The absolute value is crucial.
If the absolute value was not there then I would just have
the integral of the derivative, so it would be u
at that point minus u at that.
But here I'm taking absolute values, so for awhile the
derivative is positive up to there.
So the integral gives me this value minus this.
Then in this region the derivative is negative, but
I'm taking absolute value.
Do you see that that integral will split
into 4 separate integrals?
The total variation will be this difference plus this
difference plus this plus this.
I'll have 4--
in this example-- have 4 positive
contributions to the variation.
And that's the key quantity there.
And what a shock will do-- can I maybe mention in this figure
what a shock does, can do, would be to somehow cut out
part of the variation.
When a shock occurs it might happen that these things are
going into the shock, these are going in from that side,
and the shock kind of jumps from this value, u left, to
this value, u right following the shock
condition when it emerges.
And this additional variation is wiped out.
You see that the variation for the solution that contains a
shock is still-- this part is still there, but this part is
all downhill, so it's just this.
It's that minus that.
And the middle two parts have been cut out.
OK, so that's the total variation in a function.
Let me just take the most important example that's not
total variation diminishing, which is Lax-Wendroff.
The one even in the linear case.
The one experiment we've done with that profile, if you
remember we had a wall of water and the exact solution
brought the wall to the left with velocity c, with speed c,
but Lax-Wendroff, which did quite well staying near the
discontinuity near the jump had a
little oscillation behind.
That was essentially our one objection to Lax-Wendroff was
that oscillation showing up.
And so that's telling us that the variation--see, the
variation in the true solution is just that step.
The variation in the Lax-Wendroff answer is this
bigger step plus this one plus this one plus
this, larger than 1.
Roughly the match is that monotone methods.
So here's the general picture.
We can get monotone or I'll say, all positive coefficients
and those will be total variation diminishing.
Great.
Why don't I draw--
recall a monotome example was Lax-Friedrichs or upwind.
Upwind will be our model, so let me draw the upwind one.
The upwind one didn't increase the variation, it stayed at 1.
We had no oscillations, but it was low accuracy.
So first order.
We can't do better then first order accuracy, the order of
accuracy I'm measuring in the smooth part of the solution
because that's where I can take taylor series and match
terms and see which is the first one that doesn't match.
And Lax-Friedrichs and upwind was first order.
Then our problem is that if I allow a negative coefficient
then variation increases, ocsillations will come in, but
I can get the accuracy up to second order or higher.
And Lax-Wendroff is the example, the key example there
that I'll continue with.
OK, that's the idea of total variation.
It just gives us a quantity to work with when we see
oscillation.
And of course, I'm interested in oscillations, so this is a
linear case with a discontinuity.
I don't think it's quite fair to call it a shock.
At least, the shocks that we meet in conservation laws in
the nonlinear problems, they're created--
the ones we saw here are created for nonlinear reasons.
Colliding characteristics, whereas the linear case has
characteristics all parallel.
So here we're in the lecture, we're really into the second
lecture into the scientific computing part of our problem.
So that's the analysis part, the pd part.
Now comes the experience of quite a lot of people working
hard to identify a way to get second order accuracy where it
smooths and monotonicity or something like it where the
shocks are.
So that's essentially the goal.
Our method is going to turn out to be nonlinear.
It's going to have some--
it'll be called a flux limiter.
The key quantity will be the flux.
And for a smooth function when u right and u left are very
close, the flux is small, we don't have a problem.
But at a shock where there's a significant jump in u and
especially in f of u, the flux limiter is going to come in,
into the discrete method.
Let me began discrete methods with a new
word, finite volumes.
And with a new interpretation of capital U.
So the discreet value, capital U jn, which up to now we've
been thinking of was that's an approximation to the true
solution, j delta xn delta t.
Now because we're dealing with
discontinuities we better integrate.
We integrate over a delta x interval.
We take the average over a little delta x interval
instead of hitting just the point.
So the delta x interval from j minus 1/2 to j plus 1/2, we
take the average of u at the time, n and delat t and that's
what finite volume--
and let's think of as capital U, the discrete approximation.
The appropriateness of that is in the integral form of the
conservation law.
So the integral form of a conservation law leads us to
that natural discretization.
d by dt of u, of the average.
You're keeping your eye on the integral form of the
conservation law.
Just come back to it again.
The derivative of the mass, the integral of u, so in the
discrete case the difference of capital U plus the
integrated term, a flux at one end of the integral minus a
flux at the other end.
And notice that we are getting for the flux, 1/2.
So a staggered grid is presenting itself naturally.
And so that f j plus 1/2, the fluxes, which is the critical
thing is on the j plus 1/2 grid.
All right, now this could be exactly defined as the flux at
that grid point.
Capital U squared over 2 in the Burger equation at that
grid point.
That would give us a totally straightforward, first order
numerical method.
First order, I see first order here in time and it might be
better than that in space because in
some way that's centered.
But if we want higher accuracy, if we want to
include Lax-Wendroff for example, that includes several
values of u, we go to a more general flux function.
The whole job here is create a smart flux function.
So a flux function depends not only on u in the center, uj--
at time n, this is all at time n.
It can also depend on values to the left of it and values
to the right.
It may use those values, it may not, but we allow that
possibility then that the flux--
so this is the numerical method that you aim for.
Is to choose a flux function and of course, this should be
consistent with the original problem and that simply says
that if all the u's were the same that the numerical flux
function would give you the flux function
in the given equation.
OK, so this is our general pattern.
And now the question is really, maybe the first step
is I would like to write upwind and Lax-Wendroff.
I'd like to find their flux functions.
I guess I'm saying that when we look at an equation in this
format, we can take methods that we have already used and
see, OK, even in the linear case, the
one-way wave equation--
what are we picking for the flux function?
OK, now the one distinction is--
the one-way wave equation we had a constant and upwind for
example, went in the direction decided by that speed, c.
If the speed reversed sign so that the waves were going the
other way, the wind is coming the other way, the upwind
direction is the other way so our method would change.
Let me just write down the flux function that turns out
to come for upwind and then for Lax-Wendroff.
OK, so for upwind the flux turns out to be f, turns out
to be 1/2 of--
because it's at this, this is fj plus 1/2, really.
It's at the center point.
The flux function will be, I'll put in the a.
a delta t--
sorry, a over 2uj plus 1 plus uj.
Now plus, or rather minus actually, comes absolute value
of a/2 times uj plus 1 minus uj.
If we're patient we can track down the two cases, a positive
and a negative, and recognize that this form switches
direction correctly.
The difference method.
So it keeps us upwind in other words.
If the wind changes direction our flux changes appropriately
to follow that direction.
So that's an OK flux function.
It's a monotone one and it's first order.
Now let me put in Lax-Wendroff.
What would Lax-Wendroff be?
Would be the same, the same upwind flux, if I call that up
for upwind and then a correction, a Lax-Wendroff
term that gets the second order accuracy.
OK, and that turns out to be, it turns out that I multiply
this term by that ratio, r.
By the absolute value of r.
So it turns out that it's the Lax-Wendroff thing, to get
that extra accuracy.
You might remember, had an extra delta t over delta x.
And so it's this same thing times a over 2.
In the nonlinear case it'll be f of uj plus 1.
f at uj j plus 1 minus f at uj.
Next time I'll come back to this point and pin down this
more carefully.
All I want to do in completing today is introduce this idea
of a flux limiter.
So what people say when they look at Lax-Wendroff is that
this term that's been added on is fine in the case of a
smooth solution, but in that case of a jump like our
picture here, that term is increasing variation.
It's an anti-diffusion term.
Instead of smoothing things out the way f upwind does,
that extra term has the effect of increasing the jump.
It's anti-diffusive and therefore when it's a
significant term the flux limiter cuts back on it.
The result will be a numerical method, which is not linear
anymore even if we apply it to a linear problem.
The result won't be strictly linear and it will produce a
profile that doesn't have that oscillation.
So it has what we want.
It gives us good answers when the solution is smooth, it
stays near the shock the way Lax-Wendroff does, but it
doesn't oscillate.
OK, so what's the flux limiter?
The flux limiter is, I multiply in here by some
function, the limiter function, phi j plus 1 and I
multiply in here by phi j.
In other words, there's an extra factor, this flux
limiting factor.
So the question is, how to construct that?
What formula do we use for the flux limiter?
So the flux limiter, after a lot of experiments is chosen
to depend on--
so we want it to be 1.
So what do we want?
We want it to be 1 or close to 1 and smooth part and we want
it to be greater than 1 at a discontinuity.
OK, so our question is, how do you identify that
discontinuity?
Well, the convenient way has turned out to be look at the
ratio-- look at differences.
Look at the ratio of-- so what comes in here will be the
ratio r say j plus 1/2.
That'll be the ratio of u, say j plus 2 minus uj plus 1 to uj
plus 1 minus uj.
All that time n.
That ratio tells us, that ratio is close to 1.
So the ratio r is close to 1 for smooth and our function,
phi of r then, we want our function phi when the ratio is
1 to give the value 1.
So if I try to draw a graph of-- let me draw a graph of
peoples flux functions--
of the phi of r.
So here's a graph of phi of r.
First on the graph I'll put upwind and Lax-Wendroff
without anything.
So Lax-Wendroff would take--
phi Lax-Wendroff is identically 1.
This is the ratio r and this is phi.
So Lax-Wendroff before fixing had no flux limiter function
and it just had 1.
And actually, if I can keep the same formula, phi upwind
would just be zero.
Then it'll turn out that--
let me draw these ratios.
There is a region here and here where the flux limiter is
doing the right thing.
So in that region it's acting the right way.
It's limiting the anti-diffusive flux.
It's controlling the variation, it's controlling
the oscillation.
So various people have proposed--
let's see, no.
I haven't got this right.
I think it's- because I know that part is wrong for
Lax-Wendroff.
So maybe it's there.
No, I've forgotten--
sorry that I'm not, I've lost this because I know that one
possibility is this as the flux limiter.
One possibility would be this one and then there's another
possibility that is OK that goes through here.
Let me stop for today with this idea of a flux limiter
and you see that it's a task to fix the method to do the
right thing at smooth parts and also the right thing at
discontinuities.
OK, I'll pick up on that next time.
Thanks.