# Session 11: Numerical Methods of PDEs: Finite Volume Methods 2

Flash and JavaScript are required for this feature.

Download the video from Internet Archive.

- About this Video
- Playlist
- Transcript
- Lecture Notes
- Download this Video

Description: This session continues discussion of finite volume methods, and works through an example of upwinding using a traffic jam simulation.

Instructor: Qiqi Wang

The recording quality of this video is the best available from the source.

Session 1: Course introduct...

Session 2: Numerical Integr...

Session 3: Numerical Integr...

Session 4: Numerical Integr...

Session 5: Numerical Integr...

Session 6: Numerical Integr...

Session 7: Numerical Method...

Session 8: Numerical Method...

Session 10: Numerical Metho...

Now Playing

Session 11: Numerical Metho...

Session 12: Numerical Metho...

Session 15: Numerical Metho...

Session 16: Numerical Metho...

Session 21: Probabilistic M...

Session 23: Probabilistic M...

Session 24: Probabilistic M...

Session 25: Probabilistic M...

- Download English-US transcript (PDF)

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.

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]

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. 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.

Numerical Methods of PDEs: Finite Volume Methods 2 (PDF)

## Free Downloads

### Video

- Internet Archive (MP4 - 116MB)

### Caption

- English-US (SRT)