The following content is provided under a Creative
Commons license.
Your support will help MIT OpenCourseWare
continue to offer high quality educational resources for free.
To make a donation or to view additional materials
from hundreds of MIT courses, visit MIT OpenCourseWare
at osw.mit.edu.
QIQI WANG: So today, we're going to be continuing finite volume.
We started on finite volume on Monday.
And just a very quick recap, finite volume
doesn't solve all the differential equations,
but it solves a class of differential equations
that we care most, that is conservation laws.
Conservation laws are differential equations that
can be returned in this form.
And today, I'm just going to be more specific.
Let's say, rho is the density of a conservative quantity.
The time derivative of rho plus the divergence of a flux,
which is a function of rho, is equal to 0.
So this is a general conservation law
when there is no [INAUDIBLE], when there is no generation
or it's appearing off any material within the control
volume.
Everything that increases has to come from outside the control
volume.
Anything that decreases has to go out of the control volume.
And this is the differential form.
And in one dimension of cos, the divergence
is just spatial derivative.
It's just a partial x of F of rho.
And today we have going to be solving
one of the equations like this.
And in order to apply finite volume,
we need to make it an integral form.
We have to really integrate this in space.
And instead of a divergence operating space,
we get this-- d dt of an integral over any control
volume.
Let me describe it as omega.
Omega is a control volume.
So we are basically integrating this whole thing in omega.
And this is the first thing we get.
And the second thing we get is that now we
are doing this in impossibly multiple dimensions,
like last week [INAUDIBLE].
Now if you integrate the divergence of something--
you're integrating a divergence of a vector
field-- in a control volume, what do you get?
A listing something like the integral
of a divergence of a vector in a control volume
equal to something [INAUDIBLE].
Something can relate a divergence.
[INAUDIBLE] the divergence 0.
[INAUDIBLE] is so cool.
There are two theorems in method popular that I
can see that are very important, both
relating a spatial derivative to a spatial integral.
And the divergence theorem is the more fundamental one.
I can feel a more fundamental one,
because it basically tells you the spatial integral
of a divergence can be expressed with something that doesn't
involve derivative at all.
It can be represented as the integral over the surface
of the control volume.
The surface normal starts with just the flux vector
is equal to 0.
So divergence theorem says that the integral of this
is equal to this-- divergence theorem.
So what you see is that this is the 2D or 3D version
of the integral form of the conservation law.
What we are saying is that now, again, this
is the total amount of mass, or stuff, within the control
volume.
And today we're going to be doing half.
We see that the total [INAUDIBLE]
half within a control volume.
And this is the flux through anywhere in the surface.
So there is a normal.
And this normal is also at normal of the control volume.
So n points towards the outset.
So this F of rho is the flux towards the outset.
That makes sense, because if F points outward--
if the flux is towards the outside-- then
this term would be positive, because F
is going to be aligned with n.
This term is positive.
And the time derivative will be negative
because things are going out, the stuff should decrease.
If asked if the flux vector points inward, then
it is going opposite to the outward point of normal,
opposite to n.
And this integral would be negative.
So if this thing is negative, then we have to add to 0
so this can have some positive, which means the amount of stuff
is going to increase as a function of time.
Is it clear?
This is now generalizing that in 2D or 3D.
And of course, we can still write that in 1D.
In the one dimensional special case,
what is this partial omega?
Partial omega is the boundary of omega.
What is the boundary of the one dimensional control volume?
The point, or two points-- two points, right?
And the two points, I have R and L. And outward pointing normal
is going to be just a 1 at the right hand side at R,
is going to be minus 1 at L. So I'm just
going to be saying this is F of rho at R minus F of rho at L
is equal to 0.
Or as we said on Monday, the time derivative of total stuff
is equal to the stuff coming in from the left minus the stuff
coming out from the right, if you moved both of this stuff
[INAUDIBLE] side by side.
So this is the kind of differential equations
we are solving, and the scheme is finite volume.
And finite volume basically follows exactly this formula.
Instead of integral over control volume,
we are storing the average of the control volume, which
is simply the integral divided by the size of the control
volume.
It could modify the cell average with the size of the cell,
with the integral.
And in order to think of the DDP of the cell average,
you just divide this also by the cell volume, or cell size.
So you get the evolution equation
that allows you to control using finite volume.
Now this equation is not approximate This is exact.
What we need to approximate in finite volume is this flux.
This flux is a flux evaluated after the cell boundaries
or cell [INAUDIBLE].
And the finite volume is [INAUDIBLE].
You only have the cell average.
So you need to approximate the interface value from the cell
average.
And that is the approximation we have been talking about.
So that's kind of a review.
The first new thing we're going to talk about
is, actually, why do we use finite volume.
We use finite volume because it captures shock waves,
because if there is discontinuity in the solution,
finite difference doesn't work.
I'm going to show you with a-- what?
Cancel.
I'm still testing around with the code you wrote.
But let me go to today's directory,
and I'm going to show you two demos.
I'm going to show you a simulate-- so first,
I'm going to tell you what this is doing.
So I have two functions, one is simulate fv,
and one is simulate fd.
So these two functions are identical,
except for this single letter.
One is fd, one's fv.
But otherwise, I'm just doing the same.
I'm using [INAUDIBLE] 45 to integrate in time.
Conditions on that I'm [INAUDIBLE]
having space, using either finite volume
or finite difference.
I see there is a question on what this ''at'' do.
This "at" is basically take a function I wrote,
which is either du dp fv or du dp fd
and feed that into the ODE 45.
Anytime you use ODE 45 you need this "at" to make it work.
Now the difference between finding fd and finding fv
is that in this we're using the standard finite difference
to solve the differential form of the Burger's equation.
So du dp is equal to minus u times du dx.
This is du dx.
Now you should be familiar with this.
I'm taking u minus a specific version of u--
this is assuming [INAUDIBLE] conditions-- divided by dx.
So this is a [INAUDIBLE] space on a difference of du dx.
So here I am solving the differential form
of the equation of the finite difference.
In finite volume-- let's review what we do in finite volume.
du dx [INAUDIBLE] flux.
Now, this is not using the differential form,
but using the [INAUDIBLE] form because we
are computing the flux.
And here, I'm assuming the opposite direction
is always the left.
Here, I'm assuming like I have an initial condition that
is opposite, and it's easy to illustrate
using a case like that.
So the wind always come from the left, u is always positive.
That means when [INAUDIBLE] it is a minus [INAUDIBLE].
And I'm doing du dt equal to flux at the left minus flux
at the right divided by dx.
So this is finite volume.
So let's see, what is the key difference between the results
of these two schemes.
And in either one, I'm starting with something positive,
so sine x plus 1 over 2 because I'm assuming also
in finite difference and finite volume,
the wind comes from the left.
So that I'm doing [? backwarding ?]
space is a finite difference or finite volume.
So I have a positive solution.
I told them to see how it does.
So I'm going to [INAUDIBLE] on a different sort,
and somebody new will tell me, what
is wrong with a finite difference solution.
Watch carefully.
I'm going to start.
You need to tell me what's going wrong.
So it looks like any Burger's equation.
I mean, we have a shock wave.
Things are going.
So what's wrong with this?
Is anything wrong?
Solution [INAUDIBLE].
Let me run again.
Anything looks wrong?
Many different circles it'd be tracing, but it's not.
AUDIENCE: [INAUDIBLE]
QIQI WANG: You have sharp eyes.
Now, let's look at the finite volume solution.
Exactly the same, but using finite volume.
What's the difference?
This [INAUDIBLE] the shock wave is correct in this case.
While the shock just [INAUDIBLE] it forms,
the finite difference just stands there.
The speed of the shock wave is finite difference,
is wrong in finite difference.
So by using the differential form of the equation
and disregarding the conservative form,
you are not guaranteed to have the correct speed of the shock
wave.
It actually solves the other part correctly.
The only part that it doesn't capture correctly
is how the shock wave is supposed to behave.
And if we use finite volume, because we
are looking at the conservative form of the differential
equation of this form, the other integral form.
This form of [INAUDIBLE] is to have shock waves.
Solve this one, and you take the divergence of the [INAUDIBLE]
that is discontinuous, for it acts derivative of something
that is discontinuous.
That is no longer a correct thing to do.
And therefore, you get wrong behavior
when there is discontinuity like shock waves.
Any questions?
Oh, and by the way, another thing is this.
I should let you look.
In addition to the shock wave, there
is another thing that is not right in the finite difference.
Lets look again at the finite volume.
The Burger's equation's a conservation law,
which means that [INAUDIBLE] under the curve should be
[INAUDIBLE].
The integral of the solution u, of course, the whole domain,
should be the same, because anything that comes out of here
should come back into here.
It's like when you step out of the door, you come back.
You're not disappearing.
But if you look at the finite difference,
at this point, it [INAUDIBLE].
Things just disappear out of nowhere.
The solution is not conserved.
So two things that finite difference is missing.
One is, it doesn't have to behave as a shock wave.
Two, it is the solution is not conservative.
Any question on why we do finite volume?
So this is why we do finite volume.
And the next question is another question,
probably, you have when you are going through the last lecture
and reading through the material is, why the hell
do we do upwind?
Why do we look at the upwind direction?
So the upwind of flux is this.
So we look at a few volumes.
k minus 1, k plus 1.
The finite volume, we store the average value in the volumes.
So we have k minus 1/2 and k plus 1/2
as the volume interfaces or the cell interfaces.
And the upwind flux is basically saying F at k minus 1/2
is either the F of rho at k minus 1 or F of rho at k,
depending on where the wind comes from, depending on if,
should say-- let me use C as the speed of the wind.
So if c is greater than zero, the winds
have positive speed, which means the wind goes
towards the right.
And then you look at the upwind direction,
you should be looking at when you are computing
the upwind value of the plus [INAUDIBLE] minus 2,
you should look towards the left, [INAUDIBLE].
So when k is greater-- so when C is less than 0,
that means the speed of the wing is negative,
opposite direction is towards the right.
The winds are [INAUDIBLE] opposite direction for this,
you should be looking at the right
and use the k value instead.
And here, c is what?
c is this.
c is the c that appears in this.
c is the value that appears in the differential form
of the equation, over here.
So as you can imagine, if we have
[INAUDIBLE] like the [INAUDIBLE] vacuum equation we have,
this c exactly determines whether your solution
moves towards the left or moves towards the right.
So if this c is positive, the solution will go to the right
and winds go towards the right.
If c is negative, it goes the opposite direction.
But in general, in a conservation law,
how do we convert our conservation law,
which is partial F rho partial x equal to 0, into to that form?
But if we convert this into this one, what is c?
Can we always convert the conservative form
into this, like a vacuum scheme looks like, leaving
a vacuum looking like form?
AUDIENCE: Just take the derivative.
QIQI WANG: I just take the derivative, exactly.
If I take c equal to partial-- not partial.
I do too many pd's.
dF d rho, then I get my C. In the Burger's equations,
c is [INAUDIBLE] rho because F is equal to half of rho square,
right?
So when I take derivative, I just get [INAUDIBLE].
As you see if you're looking at all the equations,
you'll see what we did [INAUDIBLE].
And by the way, you may be thinking
of why did I use these here.
Because if you look at a gas dynamic equation, the C you get
is going to be the speed of sound.
OK.
Why do we need do upwind?
I mean, the real mathematical reason
is something I'm going to talk about [INAUDIBLE].
But, like, I think Dr. [? Vikram ?] [? Dag ?] is going
to be giving the lecture right after the midterm.
But he's going to talking about stability.
And this is for stability.
I'm just going to demonstrate it empirically
in the shared folder we had-- not this one--
we had from last lecture.
OK.
So I had a driver here.
And I have a bunch of people's schemes.
Who's scheme is working by the way?
Who knows?
Which one of them should I use?
Who kind of knows more or less your script is working?
OK.
Yours?
OK.
[INAUDIBLE]
AUDIENCE: [INAUDIBLE]
QIQI WANG: OK.
All right.
So I'm going to change this to-- OK.
Thanks.
So I can run it here.
And I get a solution.
I mean, you probably didn't treat the boundary condition
correctly, but that's fine.
So there is something on the boundary
that is not working correctly.
I think that you might have set the [INAUDIBLE] here to be 0.
So you just set the [INAUDIBLE] here
to be 0 instead of [INAUDIBLE] the [INAUDIBLE]
blowing with the [INAUDIBLE].
Because I'm having a [INAUDIBLE] find the word, so [INAUDIBLE],
right?
[INAUDIBLE]
But that's fine.
OK.
So I think it should be relatively easy to modify.
F-- yeah, yeah, yeah.
OK.
So I just need to treat the last one.
So F_k_plus_half(and), right?
The last one is equal to-- OK.
I will just cheat here, F_bad_k(and).
Because I know my solution is positive,
so I just take the left one.
OK.
So here is k_plus_half.
AUDIENCE: That should have taken care of--
QIQI WANG: Oh, that should have taken care of that.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Sorry?
k_plus_half of i minus-- that is end, right?
AUDIENCE: Yeah.
That's [INAUDIBLE].
QIQI WANG: I think you shouldn't be-- F_bar_k, right?
OK.
And else it is equal to F_bar_k(1), right?
OK.
So let's see if that works.
Oh, that works.
OK, nice.
OK.
So we have this scheme over here.
We have an upwinding scheme.
So as we did in finite [INAUDIBLE],
will know that backward in space or forward in space anything
you [INAUDIBLE] have a one-sided [INAUDIBLE],
you get first order accuracy.
But if you take essential difference
or take the average of left and right,
you get secondary accuracy.
So why don't we use that?
I'm going to change your code to do exactly that.
I'm going to instead of doing all this conditional thing,
I'm just going to do this is equal to half of F_bar_k
and F_bar_k plus 1.
This thing's divided by 2.
So I'm not going to be using your conditions.
I'm just going to do k_plus_half equal to actually k_plus_half,
k_plus k_plus over 2.
And instead of doing this condition--
apologize first for going to make the [INAUDIBLE] unstable.
But k_plus_half(and) is going to be F_bar_k(and) plus F_bar_k(1)
divide by 2.
So now if I do this, I always set the interface value
to be the average between the two [INAUDIBLE].
What do I get?
Now, it looks like correct.
Oh.
All right.
So what we see here?
We see here in the beginning, it looks fine.
Let's look more carefully what actually happens.
Let's set dt equal to 0.001 to be more-- let's see
what happens.
OK.
Usually, the evolution- let me [INAUDIBLE]
it, so it doesn't [INAUDIBLE].
[INAUDIBLE] I'm going to do ylim(0,2),
so that doesn't [INAUDIBLE].
OK.
In the beginning, we see it is behaving correctly.
And I am probably getting secondary accuracy, right?
But as soon as the [INAUDIBLE] is formed,
something bad is going to happen.
Now, [INAUDIBLE].
All right.
[INAUDIBLE]
So, yeah.
So we've got these calculations.
And it looks pretty much like in [INAUDIBLE] analysis
you get this [INAUDIBLE] phenomenon, like right or wrong
[INAUDIBLE].
So this is actually why we need to take an upwinding scheme,
especially if you have a nonlinear differential equation
with potentially shock waves.
I mean, if you go to math class, I
can probably talk more about shock waves being
a dissipated phenomenon and a second order scheme,
like essential difference that has no numerical dissipation.
And you need the numerical dissipation [INAUDIBLE].
I'm just going to demonstrate numerically.
This type of thing is going to happen if we don't do upwinding
and you have shock wave.
So for something that doesn't have shock wave
it's probably fine as long as you use an appropriate time
integrator.
OK.
So today-- I mean, this is why upwinding-- we
are going to be doing something like this.
So I tried to stimulate distance.
This is a lot of parts, right?
And if you do a simulation of individual parts,
yeah, you can do this.
But [INAUDIBLE] this [INAUDIBLE] it really
looks like-- I mean, each [INAUDIBLE] more
looks more or less like a molecule.
So you really get a continuum of cars or [INAUDIBLE].
And one of the ideas is to instead of treating each car
as an individual agent and simulate that, we just
look at from a macroscopic point of view
and simulate this like gas dynamics,
where each car is like a molecule.
OK.
And so a macroscopic description of the cars on a highway
is just P as a function of x and t, right?
So x is space.
t is time.
And P is the number of cars per length of the highway.
So P is a converse quantity, because cars are converse,
right?
Cars are not disappearing into [INAUDIBLE].
Yeah.
Cars are not disappearing as airlines sometimes do.
AUDIENCE: [INAUDIBLE]
[INTERPOSING VOICES]
QIQI WANG: OK.
So we are also assuming the velocity of cars
depends on the density.
Well, let's first say so P equal to 0 is an empty highway.
P equal to 1 is this.
OK.
P equal to 1 is like a parking lot.
P equal to 1 is like the maximum density
you can pack onto a highway.
OK.
So this is P. And we are assuming the velocity
of the cars is a functional P.
So again, a velocity equal to 0 is like this.
So this car cannot move when the density is at [INAUDIBLE].
And let's pick the unit of velocity,
so that velocity equal to 1 is something
like 65 miles per hour, where you would have-- or 80--
where you would be driving if there is no other cars.
Or 85.
OK.
So what would be an easy way or intuitive
way you would model u as a function P?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah.
It can be more curved.
But let's make it simple.
Did you want to [INAUDIBLE]?
So would you [INAUDIBLE] code [INAUDIBLE]
like something like this?
Or would you [INAUDIBLE] coding [INAUDIBLE] like U(P)= 1 minus
P?
This captures the right behavior.
So when P is equal to 0, you drive at 85 miles per hour.
When P is equal to 1, even if you want to drive faster,
you're stuck, right?
OK.
Now, here's my first question.
What is the flux?
What's F(P)?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Integral--
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah.
I know what you're thinking.
You are thinking about U being the d f in P, right?
So it has to be the integral.
But that's not true.
It's easier than that.
You are thinking more-- you're too smart.
AUDIENCE: [INAUDIBLE]
QIQI WANG: P times U. That's that correct answer.
Can you say why that is this case?
AUDIENCE: [INAUDIBLE]
QIQI WANG: That's [INAUDIBLE].
Yeah.
Because flux is number of cars, the amount of converse quantity
that flows through a point per unit amount of time, right?
So if you just [INAUDIBLE] the control volume like this,
this is my-- oh, I shouldn't use green, because you can't see.
If you true the control volume like this--
and the flux of this [INAUDIBLE] is the number
of cars that goes through this line per unit amount of time.
So this is cars that have density of rho
goes through the [INAUDIBLE] U. Of course, the flux is
going to be P times U.
The flux is going to be P times U, which is
P minus P squared in this case.
Because U is 1 minus P. P times U is P minus P squared.
Is it clear?
Any questions?
So we derived our differential equation, dP dt
plus d of P-- I'm just going to write F(P), it's easier.
F(P) dx is equal to 0 where F(P) is equal to P minus P squared.
OK.
Now, what is C?
What is the [INAUDIBLE] characteristic C?
Now, this becomes highly [INAUDIBLE].
What is C?
AUDIENCE: 1 minus [INAUDIBLE].
QIQI WANG: It's 1 minus 2P.
It is equal to dF dP, right?
So derivative of P is 1 [INAUDIBLE].
OK.
So we have a relation like this.
We have a relation of F versus P.
P is all the way from 0 to 1.
F is something like this, where the maximum takes 1/4.
OK.
And we have a C that looks like this.
We have C that is equal to 1 when there is no cars
and goes all the way to minus 1 when
there is a parking lot of cars.
This is 1 and minus 1.
OK.
So what is happening?
C, again, is the speed of the local solution.
It is intuitive that when there are no cars,
things move at 85 miles per hour.
So the [INAUDIBLE] should be 1.
What happens here?
When there is a parking lot of cars,
things move backwards at 85 miles an hour?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yes.
It's kind of the-- so this C here
is not-- I mean, if you look into physics
you've got two kind of speeds, like a group speed
and a [INAUDIBLE] speed, right?
So here, the cars move forward always have a positive speed.
But it doesn't mean if you look at cars in the backwards
[INAUDIBLE] point of view, it doesn't
mean the wave that appears in the flow of traffic
is going to move forward.
Sometimes you are going to see waves actually
moving backwards.
I mean, we're waiting to see some examples when it actually
solves that numerically.
OK.
So the task for you is to hold up
a simulation of this thing in a periodic [INAUDIBLE] solution.
So let's say this is our highway.
This is our highway.
It just goes in circles, right?
Assuming you have a bunch of cars driving in circles.
AUDIENCE: I could never get off the [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
QIQI WANG: Huh?
AUDIENCE: NASCAR.
QIQI WANG: NASCAR, yeah, yeah, right.
Yeah.
We have a lot of cars driving in home circles.
And let's just [INAUDIBLE] how these work.
And then at sometime, we are going
to put a red light over here.
And say, OK, the cars stop here.
What is the right way to model that numerically?
You can set the flux [INAUDIBLE] to 0, right?
I mean, at some point, that's a [INAUDIBLE] simulation.
And for example, time equal to 1,
that's [INAUDIBLE] 0 [INAUDIBLE] before the red light.
And at some point, let's set the light [INAUDIBLE].
But let's do something simple in the beginning.
Let's just [INAUDIBLE] the solution
of [INAUDIBLE] with no red lights or anything like that.
I'm going to give you the initial position.
And we'll see how it evolves.
And then we are using final volume.
And you are free to take-- the [INAUDIBLE] case
we had in the last lecture is in the shared folder.
And I'm just going to make another shared folder here.
OK.
I'm just going to take a [INAUDIBLE] and the du dt.
I'm just going to take the [INAUDIBLE] stamp [INAUDIBLE]
to this.
Because I know it's working.
And [INAUDIBLE] start by modifying [INAUDIBLE].
So instead of u_bar, you're going to solving for P_bar.
And what is the first thing you would want changed in this?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, yeah.
Yeah, make it upwind.
Yeah, right.
To make it upwind there are [INAUDIBLE].
Huh?
AUDIENCE: [INAUDIBLE]
QIQI WANG: It's no longer that equation.
The equation is [INAUDIBLE].
But there is only one line we can
change to solve for a different equation and find the volume.
What line is that?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Second line.
The only thing that satisfies a different equation is the flux.
You change the flux, you are now solving a different equation.
OK.
Of course, another thing is what is
upwind direction [INAUDIBLE].
In [INAUDIBLE] equation, we know U is positive plus [INAUDIBLE]
the right.
U is negative [INAUDIBLE].
Now, this is different.
Now, this is different, because we have [INAUDIBLE].
What is left and what is right all
depends on the [INAUDIBLE] instead of the [INAUDIBLE].
So any questions?
OK.
If you are done with the-- If you get the solution going,
try to put a red light somewhere.
Let me kind of illustrate how I should modify the driver.
So I have this here.
That, I didn't call it [INAUDIBLE].
Because it's no longer [INAUDIBLE].
Let me [INAUDIBLE].
trafficDriver.
OK, trafficDriver.
And let me see.
So in this ODE 45, I'm going to set t0=(i-1) times dt and t1=i
times dt.
So this is like the beginning and end of this interval.
And when I'm solving it, I'm going
to be solving it from t0 to t1.
So now the difference is if I look at-- so let me do this.
OK.
So I just want to be able to passing the actual time.
So I'm just going to be passing also time into this.
So I also have a time into this.
And then now if-- OK, so let me just
go through what you need to modify.
So first of all, I'm going to modify this,
which is the different flux.
Right?
Instead of half of u squared, I have u minus u squared.
OK.
k_plus_half is still the same.
I'm going to-- so if I don't modify anything,
I think it all wrong.
But like once a shock wave forms, it'll go unstable.
So let me try it.
I doesn't even wrong.
[INAUDIBLE] dimension matrices being concatenated
are not consistent.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, OK.
I'm [INAUDIBLE] different [INAUDIBLE].
Yes, thank you.
I save this.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, i'm solving-- yeah-- from t0 to t1.
Yeah, thank you.
OK.
So this thing goes.
And we're going to see [INAUDIBLE]-- wait.
I shouldn't initialize it to 2, right?
Because 1 is parking lot.
What is 2?
OK.
So I need to also modify the initial condition.
So what I can do is I can set this, let's say, divide by 4.
OK.
So this is going to be my initial condition [INAUDIBLE].
OK.
I also want to-- my ylim, when I show my solution,
I just want to go from 0 to 1, right?
OK.
So let me run again.
OK.
So this is how my solution looks like.
This part goes really fast.
[INAUDIBLE] unstable.
And I need to change this to upwinding.
But we have the [INAUDIBLE] part of the solution correct.
Every [INAUDIBLE] towards a positive direction when
the [INAUDIBLE] move faster.
That's [INAUDIBLE] car density, cars move slower.
So sometimes, a shock wave is going to form.
Because the faster parts is going
to [INAUDIBLE] to the slower parts and decelerate.
So it becomes a shock wave.
So let's start to modify this to an upwinding scheme.
We computed F_bar.
And we have a c.
c is equal to 1-u_bar times 2.
OK.
And the c at k_plus_half is equal to c plus-- OK.
Let me do like this.
So in the condition, instead of u_bar, I'm going to put c.
OK.
Here, too.
Instead of the u_bar, I'm going to put c.
OK.
So this completes my upwinding.
And let's try it again.
We still have the solution.
But the lower density goes faster.
High density goes slower [INAUDIBLE] high density,
but now it's [INAUDIBLE] stability.
So that's good.
And if you want to get better solution,
I can change this to 256.
And I round this.
And I [INAUDIBLE] the top [INAUDIBLE].
OK.
So now, let's say at this point, so at either the beginning
[INAUDIBLE], I have time equal to 1, we get a red light.
And at time equal to 2, we turn the light to green.
[INAUDIBLE] there.
OK.
So we do that.
[INAUDIBLE] then that [INAUDIBLE]
what we need to do to make that happen.
The key thing is now you have time.
You have t available to you.
So you can basically do a conditional statement on t.
If t is greater than 1 and t is less than 2, set the red light.
[INAUDIBLE]
AUDIENCE: What about [INAUDIBLE]?
QIQI WANG: Oh.
Thank you.
That is wrong.
So it should be F_bar_k.
Yeah, it should be F_bar_k.
Thanks.
Yeah.
I should just rename it to something else,
so that-- I'm just going to rename it as 2t.
All right.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Which one [INAUDIBLE]?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Sorry?
AUDIENCE: [INAUDIBLE]
QIQI WANG: This one?
AUDIENCE: Yeah. [INAUDIBLE]
QIQI WANG: Yeah. [INAUDIBLE] plus 1, right?
So because I'm looking, I [INAUDIBLE] the interface.
I want to look at the interface [INAUDIBLE]
decide which side I should upwind the flux at k_plus_half.
In order to make a decision, I want my [INAUDIBLE] speed
at k_plus_half.
And the right way to [INAUDIBLE] the other [INAUDIBLE]
here k_plus_half is take the average of the [INAUDIBLE]
speed [INAUDIBLE] at k(and) k plus 1.
Just to make the code easier to understand,
I just want to change-- I think I want to replace-- oh.
I think I want to replace my i with k.
[INAUDIBLE]
All right.
Yeah.
So instead of-- it's easier to replace [INAUDIBLE]
k_plus_half than [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
QIQI WANG: Why did I half t?
I wanted to half t, because what I like you to do now
is to set a [INAUDIBLE] when t is [INAUDIBLE].
OK.
So what happens on the red light?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Huh?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, the [INAUDIBLE] goes to 0.
That's one way to say it.
Another way to say it is no cars can pass, right?
If no cars can pass, then the [INAUDIBLE]
pass through the red light is what?
0, right?
So what I can do here is that before I compute k minus half,
I just want to say if k is greater
than 1 t is less than 2, F_k_plus_half at the end--
so k_plus_half-- the last one is basically the right hand
and also the left hand-- is equal to 0.
I think I [INAUDIBLE].
OK.
That's it.
That's what I need to do.
Now, in the trafficDriver, let me set it to-- ah, ah, ah.
Let me copy this trafficDriver for the [INAUDIBLE].
Yeah.
I think we're overwriting [INAUDIBLE].
And it's going to save trafficDriver [INAUDIBLE].
OK.
Now, instead of 1 over dt, I'm going
to simulate it for five time units,
so that I actually see the light.
OK.
We see the same thing here [INAUDIBLE].
At sometimes, the light [INAUDIBLE].
Yeah, OK.
See, the lights is in effect here.
All the [INAUDIBLE] simulates here.
And then we have even [INAUDIBLE]
that actually goes backward.
OK.
Now, there is no cars moving.
All the cars are stopped over there.
So this is a parking lot while there is no cars here.
And at [INAUDIBLE] equal to 2 or [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
QIQI WANG: Did I-- did I-- what's happening?
OK.
If t is greater than 1 and t is less than 2--
AUDIENCE: [INAUDIBLE]
QIQI WANG: And never goes passed 2, are you sure?
AUDIENCE: [INAUDIBLE]
QIQI WANG: OK.
I think I know what's happening.
So I actually have come to a situation
where-- hm, let me think.
AUDIENCE: [INAUDIBLE].
What
QIQI WANG: Yeah.
Let me do this.
I'm actually encountering a situation that cannot be
handled by this [INAUDIBLE].
What is happening is that [INAUDIBLE]-- let me do this.
I think I can change this around by doing this.
I can change this around by not letting the red light
to stop for that long.
So I only let the red light to-- if I only
let the red light go from 0.2, I think
I'm going to get rid of the problem.
The problem is that when the red light happens here,
the flux-- the P at one side shock wave is 1.
At the other side of the shock wave, the P is equal to 0.
So both values of P [INAUDIBLE] equal to 0.
Because the flux is P times 1 minus P or P minus P squared.
So the other P is 1 or zero.
The flux is equal to 0.
That's why no matter if you take upwinding or downwinding,
you'll get 0 flux.
So this is a case where [INAUDIBLE] upwinding
wouldn't work.
It's like a singular case where [INAUDIBLE] upwinding
wouldn't work over here.
So what I'm going to do here is am going
to set something artificial.
And let's say when there is a red light,
I think I'm just going to set the flux equal to 0.01.
So that there is a small leakage of cars over the red light,
so that we actually want to see what
is going to happen [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah.
So if you take [INAUDIBLE], you're
going to learn about the entropy conditions of [INAUDIBLE].
And the entropy condition tells you
that there is a physical solution.
There is a non-physical solution.
And over here [INAUDIBLE].
But then over here it is a physical solution.
While this kind of an [INAUDIBLE]
is going to be a non-physical solution.
And that cannot be dealt with by the [INAUDIBLE] scheme.
And you nee a more complex flux scheme
to deal with situations like that.
AUDIENCE: [INAUDIBLE]
QIQI WANG: There is no downwinding.
And what you need to do is this.
What you need do is you need add something
called numerical dispositing into this.
So you need to do F_k_plus_half is
equal to F_k_plus_half plus 0.5 times absolute value of c.
So instead of taking the c's, you
want to take the absolute value of c's and divide by 2.
So this is like a very similar concept to upwinding.
But it operates differently.
You need to take the absolute value of c,
but multiply that with u_bar and plus_half.
So it's minus [INAUDIBLE] u_bar minus 1 0, minus 1.
So you need to kind of-- so instead of doing upwinding,
you need to use a numerical dispositive
or a numerical diffusion to basically achieve
the impact of upwinding, but do this a little bit differently.
But that's kind of beyond the scope of this class.
So, yeah, If you want to learn that,
you're welcome to take the [INAUDIBLE].
We are going to be learning all about [INAUDIBLE] speed.
So let's see if this goes [INAUDIBLE].
OK.
So we have the [INAUDIBLE].
OK.
Now, the cars are starting to move.
[INAUDIBLE] the light goes green, the cars start to move.
To do it properly, you have to change this back
to the central average.
So you do upwinding or add numerical dispositive,
but don't do both.
And doing both, so the scheme works,
but it's not as accurate.
It's suffers from additional numerical dissipation.
OK.
So I hope you have fun playing with this small example here.
So starting from our next lecture, the next lecture
we're going to be doing a review for the midterm.
But starting from after the midterm,
we're going to be talking about first of all, stability
of solution of positive differential equations,
and then talk about finite [INAUDIBLE] method, which
is the third method we are going to start
in fall for [INAUDIBLE].
All right.
We'll see you next week.