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 integral 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 x_0.
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, two characteristics
are hitting, here 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 shock, by the way, that's a big deal--
up 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
and the characteristics speed from the left.
And notice, it's in between.
That's crucial.
That tells us that this figure 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 the shock is at-- capital
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 u_L 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 u_L,
and here I have-- the derivative of that is zero,
the derivative of x_R is zero.
Here I have a minus, so it's minus u_R.
dx/dt times that number 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
u_R minus u_L, 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 Rankine, Hugoniot 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 over t.
This solution will be a simple-- a simple ratio x over 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 over t expression.
So a simple 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 now, how do we get the finite 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 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 Burgers' 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 a while
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 four separate integrals?
The total variation will be this difference plus this difference
plus this plus this.
I'll have four-- in this example--
have four 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, where,
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 monotone 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 than 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, oscillations 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 PDE part.
Now comes the experience of quite a lot of people working
hard to identify a way to get second-order accuracy where
it's smooth 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 begin discrete methods with a new word,
finite volumes.
And with a new interpretation of capital U.
So the discrete value, capital U_(j, n), which,
up to now we've been thinking of as that's an approximation
to the true solution, at j delta x, n 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 delta 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.
I'll 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 of flux at one end
of the interval minus a flux at the other end.
And notice that we are getting for the flux, half--
so a staggered grid is presenting itself naturally.
And so that f_(j+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, u_j,
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 write a-- 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 f_(j+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 2 U_(j+1) plus U_j, now, plus--
or rather minus actually-- comes absolute value of a over 2
times u_(j+1) minus u_j.
If we're patient we can track down the two cases, a positive
and a negative, and recognize that this form switches
direction correctly on 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 u_(j+1).
f at u_(j+1) minus f at u_j.
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+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 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 at 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 U_(j+1) 1
to U_(j+1) minus U_j.
All at time n.
That ratio tells us-- that ratio is close to 1,
if-- 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 people's 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.