Recitation 8

Flash and JavaScript are required for this feature.

Download the video from iTunes U or the Internet Archive.

Instructor: Prof. Gilbert Strang

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 ocw.mit.edu.

PROFESSOR STRANG: Okay, shall we start on this review session? And I'm figuring that the main topic would be finite elements, because the last review session focused on the truss problems. One note about codes. So probably you've looked at the CSE page, where some codes that we had earlier are put there, section by section. Now people are contributing a truss code, so I don't see any reason why that shouldn't eventually go onto the CSE page. If you have other codes, like finite element codes, for specific problems, think about, if you get them with some comments, if they're appropriate to go on the CSE page, that's certainly a possibility. And the homework solutions, which are now on the 18.085 page, will migrate to the CSE page, for their permanent status. So that's just to say, this is a page that's growing, and any help is very, very welcome.

So, I'm ready for questions about finite elements, maybe specifically, or other things, too. Thanks. Yeah, lots to say about finite elements. Yeah?

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: So this is finite elements for heat transfer. Okay, which course was that?

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: Okay. Right, free or fixed boundary conditions. Okay, I see. So the question is, how do we handle boundary conditions which are not something equal zero, but something equals a constant, or something is given. Right. Good point. So, somehow when we're given numbers, they're going to-- Input, source terms, they're going to show up on the right-hand side. Right? Eventually, they'll be part of the right side vector F, in some appropriate way. So that we'll be solving KU=F, where U will represent the unknown displacements, whatever they are in the problem.

So let's see, what do I say about getting things correctly onto the right hand side? My general thought would be, so we create a matrix that I might call K_0. That's a matrix where we have not knocked out any rows or columns. It's probably, almost surely singular. It's got all possible unknowns, even those that are not actually unknown, associated with the columns of that matrix. So we start with that. Start with that. Okay. Now, suppose-- Let's take things in order then. Suppose U_7 is set to zero. So what will that do to our matrix K? How do we-- Having created this matrix K, let's suppose we've got it, how do we include in it the boundary condition U_7=0? Okay, well that column of the matrix, the column number seven, I guess, or the column corresponding to that unknown, is going to go. Is that right? Yes. That's what we've done, we've struck out a column. And we should also, I guess, strike out a row. We would really like to end up with symmetric, and normally positive definite, matrices. So this would be a column row-- I'm speaking roughly, because I didn't have advanced preparation on this question. So column/row removed.

Okay, now suppose U_7 has some value A. How do we account for that? So what do we do? U_7 isn't still unknown, so we want to get that column out of our matrix. But something has to happen to account for the A. So what is that? So my picture is that in my original K_0, which had-- My original K_0, and then all the U's equal all the F's. So this had a column for everybody, including U_7. So there was the other U's, including U_7, and these F's maybe came from body forces, stuff inside the interval. Now here comes perhaps a boundary condition. Oh, I wrote it in the middle anyway.

So what happens when I set U_7 to A? So I think that becomes, I strike out that unknown, strike out that column, strike out that row, I think. But what should happen? How does the A come in? Well that A multiplied that column, right? That A multiplied-- This equation expressed equilibrium. Those equations were all valid. You could say we haven't imposed the boundary conditions yet. Now when we do-- So that A multiplies that column. So what happens now, what do you think should I do? It seems to me that I should bring A times that column over to the other side, with a minus, of course, because I put it on the other side of the equation. So I would think I do this: A times column seven. I think. I think that would seem right.

Now you're going to ask about -- you unfortunately did ask about -- a free condition. Say u'-- So what am I going to do? I'd better think here. What's cooking now? So I could prescribe u', let me say one, equal zero, let's just think about that. What happens with that? Hmm. So let's see. Well, so u' was not one of our U's, so this isn't just simply, OK assign a value, zero or A, to a U. u' being zero, well that came in as a natural boundary condition, right? If I remember my own lectures. So that u'=0 came in as a natural boundary condition. And what did that mean? That meant that I don't have to impose it. The finite element method should deal with that. This is really our w. Our dual unknown, w. It's better to think of it as w. Maybe cu' is what it is, but it's good to think of it as w.

Okay, so I think of that as something we don't have to impose. In other words, the way I'm thinking now, I don't have to think about that. I've got trial functions whose combinations will get me close to the correct answer. And the equations will naturally do that.

Okay, now you're asking this serious question? What about w=A? What about that? Or B. Let me make it something different. Okay. Let me say right off, without having prepared a good answer to this question, I'm not going to do it as well as I should. So maybe what my real answer should be, let me think through the best answer and do it another time. It's going to contribute some term, some new term. When we integrated by parts to get to the weak form, this guy just showed up. And up to now, we've just assume that B was zero, and we didn't worry. But somewhere in that integration by parts, a term is going to appear that I can't just throw away. Yeah. So let me think about that more properly. I mean that's a very good question, and the lectures didn't allow for that. Okay, if I can put that one off and give you a good answer. Because it deserves it.

Okay, other questions. Other discussion. Yeah, you have one?

AUDIENCE: I have maybe a silly question.

PROFESSOR STRANG: You have a silly question? Good. Okay.

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: Compare finite elements with finite differences. Okay, yeah. That's not a silly question. So finite differences is certainly associated with values at a mesh. Here would be one way to think of it. So the question is, finite elements, finite element method versus finite difference method. And of course, there are other methods, too.

Okay, I could think about finite differences. One way of seeing them together is, I could think about finite differences this way. My test functions, you remember, which I always integrate against. Instead of taking the test functions to be the trial functions, which is giving us our standard finite element K, I could take the test functions to be delta functions. What would happen if I take the-- So, again just think about that.

Suppose my solution U is a combination of the phis, as always. As before. But my test functions, V, are delta functions at mesh points. Because up to now, we've always taking the V's to be the same as the phis. So the V_i, the test functions, and I'll want n of them. I'll have n mesh points. My n mesh points, those will give me n delta functions, n test functions. And what happens then, when I look at the integral? The finite element deal says, take the weak form. Well, OK. And you'll see why, I'm going to take the weak form before I take the integration by parts. Why is that? Because when I integrated by parts, do you remember what the effect was? The effect was to take a derivative off of U and throw it onto the test function V. But now that I'm choosing delta functions for the V's, it's too risky. I don't want a derivative of a delta function. I better keep it all over here.

So I'll have -cU'', times V. This is what I had originally. I just took the equation, you remember how we got these equations? I just took the strong form, multiplied by V, integrated, and now I usually integrated by parts, but I'm backing off of that. Because my V is too singular to do it. So again, what would I do? My U would be a combination, as always, of these. And what will I get for these integrals? What will I get for those integrals? Well, of course, if that's a delta function at a mesh point, when I integrate something against a delta function, I get-- What do I get when I integrate against a delta function? If V_i is the delta function at a particular mesh point i, as always I get the value of that at i. At i. So I'll get -cU'' will be whatever it is. It will be the sum of U_i's, the unknown constant, times the phi_i double prime.

And what will I get on the right-hand side? When I integrate f times V_i? I'll get f at that mesh point, f_i. In other words, I think what I've got there is essentially a difference equation. I'm back pretty much to finite differences. I'm taking the value f, not integrated now, but just its value at i, because that's what I get when I do a delta function. And I'm taking this. Notice, a lot of things are different here. And there's a name for this. This is called collocation. I guess what I'm saying is collocation, this idea, is somewhere in between. I'm interpolating -- collocation -- between the finite element, the Galerkin principle, and finite differences, where I don't even think about trial functions. In between is something that really looks like-- [SOFT MUSIC] [LAUGHTER] I don't know if that wonderful music comes through on the OpenCourseWare tape, but anyway. This is MIT. We get a brass band here.

Okay, that would be a very finite-difference-looking set of equations to approximate the original differential equation. The point I want to make, this gives me a chance to make a point here. Could I use hat functions? Could I use hat functions in this as my trial functions, if I was not going to integrate by parts? When are hat functions allowed? Because hat functions have a jump in the slope. Right? Hat functions have a jump in the slope. Then my answer, you can see it coming, is no. Hat functions, they don't have good second derivatives.

So I wanted to make that point, also, for beam bending. And let me make it. So your question is leading me to another topic that I wanted to mention. That hat functions, phi' is OK. Phi double prime, second derivatives of hat functions, are basically, let me just say it in a word, not OK. So this is OK, this is for second order equations. Well, second order equations have a second derivative. Looking bad. But when I integrate by parts-- After integration by parts. And that's exactly how we've used them. We've used them for our second order equation. We integrated by parts. We had terms like cu'-- c phi primes times V'. First derivatives, no problem. But not OK for collocation, this method that I've just spoken about. Because we're not integrating by parts. Our test function V, the delta, was not smooth enough to move a derivative over there. And, now my additional not OK, for bending.

So this is just bringing us into today's lecture, but maybe it completes a point. So in bending we did integrate by parts. So what did we get? Our bending weak form, our fourth order equation weak form. Say u fourth derivative equal f(x). When I integrated by parts twice, you remember, I had u fourth times a v, times a test function, dx. And my idea, so really we're seeing it here for the first time, was to do two integrations by parts, get two derivatives off of u, put them onto v, and maybe there's a c(x) in here, too. This was our integral to be done. That's the basic integral that leads to all the entries of K. But the trouble is-- So that's in the bending problem. In the fourth order problem, our integrals are going to look like that. And we have to use functions that are good enough.

So if I use hat functions, they're not good enough. Why's that? If you suppose u and v are the same, then I'm integrating the second derivative squared. And what is the second derivative, what's its square, and what's the integral? This is really worth recognizing. So suppose phi is the same as v, and it's the hat function. And if I want to integrate some constant, or some nice function times phi'' times v'', what do I get? What's the second derivative? What's the first derivative of my hat function? It's up, right? It's the first derivative. So let me do phi'*v', is, let me draw its graph. The first derivative is zero up to there, then the first derivative jumps up to something high, up to the middle. Then it drops to the negative. And then it comes back to zero. Right? That's my first derivative. I'm OK integrating that. I'm OK squaring it and integrating. But the second derivative, if I graph phi'', which will be the same as v'', what's the second derivative of that function? It's zero, and then what is it at this point? The derivative now? Delta, right? Spikes up. Then nothing until this point. At that point what happens? Spikes down twice as far. Really way down. And then finally there's a spike up.

Okay. I can integrate spikes times nice functions. What I can't integrate is a spike times a spike. So the answer for this integral is what? It's infinite. The integral is no good. If I have a spike times a spike, you see-- Delta functions, I'm OK with a nice function, f(x) times a delta, is OK. But not OK for the integral of delta of x times delta of x. That's not OK. That one is infinite. That one is not in our space of functions that we can integrate.

OK. So this is the difficulty you meet with higher order problems, is that your functions, some of the old functions that worked fine, are no longer OK. I need two derivatives. Or I need to be able to integrate the second derivative squared. It can have a jump, but it can't have a spike. So what functions have we chosen? How would we deal with bending problems? And I guess that'll probably come, maybe quickly, into Friday's lecture. How would I do finite elements for bending? I'm not allowed to use hats, so what will I use? Do you remember that construction, Monday? Those cubics. Remember those cubics, those C^1 cubics would be OK.

Why's that? Because they have a continuous derivative. They're one degree smoother than hat functions. These were C^0 hat functions. And they were still C^0 when I put bubbles in there. I had to get up to that extra degree of smoothness, this continuous slope. Then I can take the second derivative. The second derivative may have a jump. Jumps are OK. A jump in the second derivative there and a jump there, I can multiply two jump functions and integrate. So those C^1 cubics that we constructed, you remember at every node I had two unknowns, the height and the slope. So there were two cubics for each mesh point. But the key idea was that they had this extra continuity.

And let me just say, looking ahead, these splines that are also coming Friday have one more degree of freedom. I'm sorry, one more, there's C^2. So splines, cubic splines -- I'll just put it here because it fits -- C^2 for cubic splines, that's coming Monday. That's functions that have two continuous derivatives, and the jump is in the third derivative. So those could be used, well, they could be used for sixth order problems. Of course we're never going to see a sixth order problem. I hope. But they could be used there. Or they could be used in collocation for fourth order-- Well, could they be used for fourth order problems? I'm not sure.

Anyway, the more smoothness, as I go up in the derivatives, I have to get more smoothness in all the finite element constructions. And it gets a lot harder. Actually, I can do it in 1-D, but it gets seriously hard in two dimensions to figure out function surfaces that-- This is what CAD is all about, actually. Is figuring out, how do you get smooth surfaces? The whole CAD/CAM world of NURBS and piecewise, rational things. It's a very complicated construction, to get the shape of some engineering piece approximated, including its curved boundary, by smooth piecewise polynomials, or piecewise something. That whole world of CAD is not trivial, to do that.

Well, OK, that's my speech. You got me started with your question. Did I answer your question in some way?

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: It went beyond your question. Yes it did. That's right. These are things that I may not have time to develop in Friday's lecture, because I really have to get on to two dimensions, and gradients and divergences, all that good stuff has to start Monday. So we'll see something about splines Friday.

But the whole world of different elements is sort of interesting. Actually, can I say one more thing about the elements that we've already constructed. I want to say one more point that occurs to me. And then I really will stop and answer questions. This extra point occurred to me. You remember our example where we had hat functions, and then we also had bubble functions. Right? So where did the bubble go? These were only C^0, hats plus bubbles, plus parabola bubbles, quadratic bubbles. Parabola bubble functions. So how did those go? Between two mesh points-- Sorry, I'd better leave myself more space for that bubble. OK, so one of these guys goes up, down. One of these guys, another hat goes up and down. And then between two I threw in a bubble function, like so. Right? So those were all in my list of trial functions and test functions. Two kinds, hats and bubbles. Now, that's all fine. And we can do these integrals.

But one point I sort of went past. So let's call this phi left, phi bubble, and phi right. OK. So my trial functions, my approximate solution, U, is going to be some amount of that left hat, some amount of the bubble, some amount of the right hat. And more functions on other intervals. Now, my point is just a simple one. When we only had hats, then U and the coefficient of this hat was exactly the value of U at that mesh point. Why was that? Because all the other hats were zero there. So when I only had hats, these U's, U left and U right, were exactly the values of my approximation at the nodes. That's why it looked so much like finite differences.

What I want to ask you is, so if I take the halfway point as the sort of mesh point for the bubble, what is the value of U of my approximation, at that halfway point? I just want to plug in halfway point, evaluate everything at the halfway point. So I want to find U at this halfway point. Right there. Do you see what's happening? U bubble, what is U bubble, or phi bubble at the halfway point? One. The bubble goes up to one at the halfway point. So the coefficient that multiplies it multiplies one. And you think, ah. You might think, all right, that's the value of my approximation at the halfway point. But that's wrong, correct?

What is the value of my approximation at the halfway point, here? Well, this guy is not zero at the halfway point. This one is not zero. This should have been R, sorry. These hats are not zero there. My basis is not quite as beautiful. It's not quite as beautiful, because my trial functions are not one trial function per mesh point, including halfway points, as before. No problem. I can do all the integrals, I can get the values, I can solve KU=F to find the weights. And then I just have to remember that my basis is not as perfect as before, because at that mesh point, three functions are contributing. Probably I'm making this sound like a big deal. It isn't a big deal, but you have to remember that that's the case, that that coefficient of the bubble is not your solution at that point. You have to remember that the left hat and the right hat are also contributing at that point.

Or I could rethink my functions, and make them-- Now that's the cool thing. So that's the idea. Just so you see that I could do that. Here is my left point, here's my halfway point, and here's my right point. So I'm going to have the same combinations, which are the piecewise parabolas that are continuous at the nodes. Let me show it to you. I could change from hats and bubbles to, I'll keep the bubble, because that's one at the halfway point. That's good. But I have to fix those hat functions. If I want the hat functions to be focused-- If I want the phi left function, I call that a hat function, but it's going to change shape, it won't be a conical hat. If I want a phi function, I want it to be one at that point, its own home point, but what else do I want? I want it to be zero there. I don't want it to contribute there, if I want this perfection. Or there.

So you see my hat, instead of coming down like that and being 1/2 at that point, what shall I do? I replace that hat by some combination of the hat and the bubble to make a parabola that does that. And of course, this other way will be the same. So that's not a hat anymore. What name should I give it? Anybody suggest a name? It's sort of a hat. Sombrero. Sombrero. [LAUGHTER] OK, so I now have sombrero functions, right? OK. Good. And there'll be a sombrero function, instead of the hat, it was replaced by a better hat. Now what about the one here? Just tell me, how do I draw that one? What's that function? I'll do it with dashed lines, it's zero, it should be zero at its nodes-- at the nodes that are not its nodes. So you see each one is going to come down and back up.

Really, I'm just saying you have the freedom to do this. By using combinations of those functions, I get exactly the same as combinations of these functions. The only difference is that with these new sombrero functions-- So now I have a U left sombrero times a phi left sombrero, and a U bubble times its own bubble, I didn't change. And a U right sombrero times its phi right sombrero. And the only point is that, at that halfway point, what does this equal at the halfway point? It equals that. Right. I've just done a little change of variables. I've changed the function, so that would change its coefficients, but it's still the same space of test functions. I didn't change that. My space of test functions, these are all still piecewise parabolas. They always were piecewise parabolas, but some of them were straight parabolas. Now they're more interesting parabolas.

OK, that may be a point you didn't think about, anyway. But it shows you that in constructing finite elements, it's really what combinations you get, more than what individual hats or not, or sombreros you might use. In other words, you could use either of those. If you use the hat functions, then you could build on easy integrals. Those integrals would be a little trickier, not seriously. And then the solution would be like a finite difference thing. Every U would link to a mesh point.

OK. I'll leave that and hope you guys have some questions about homework. Thanks.

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: Oh, a question about this?

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: Yeah, when I say better accuracy, when I put the hats and the bubbles together, the improvement would be noticed in between. Yeah.

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: These guys, you mean? Well yeah. I didn't change a thing. The answer would still be the same. My capital U that would come out would be no different. It's just I've chosen a different basis, you could say, to represent it. But there wouldn't be any difference, and it would still be-- So what would be the accuracy now? This is like a good review. How close will U be to the true solution? Well, I've got up to quadratic, second degree, so I would expect h cubed error. Because that's the term I'm missing. And that's the term that the cubics would get. Right. Yes?

AUDIENCE: [UNINTELLIGIBLE] PROFESSOR STRANG: The question is, why did I pick-- So I had my hats. I guess essentially, I wanted to throw in some additional functions, the bubbles, that would up the accuracy without serious increase in work. And the natural construction was that one. Yeah. Then I got piecewise parabolas. You look for a construction of finite elements that kind of comes out neatly, and this one did. There was another question.

AUDIENCE: [UNINTELLIGIBLE]

PROFESSOR STRANG: Oh, these parabolas would produce-- That's a good question. Yeah, these parabolas would produce a delta function right there. You see, they have a jump in slope. Even if they're parabolas, they don't have smooth slope. So this could not be used for a bending problem.

Or let me say it a little differently, because people commit crimes. I call them variational crimes. People have tried. You know, it's a heck of a lot of work to construct smoother trial functions. So people say, well, just ignore that delta function. Put it in anyway, integrate elsewhere, and just skip the delta function. So that would be called a nonconforming function. There's a very famous way to create nonconforming function. And that's very hot stuff, actually. Professor Peraire and others in Aero are world experts on that. Those are called discontinuous Galerkin. Discontinuous Galerkin. This is fun. You see what this world is about.

So what's discontinuous Galerkin? Those will be ones where there's a jump in the function. Here I don't need much space to show you discontinuous Galerkin. So there's a mesh point, there's a mesh point. Discontinuous Galerkin might do this. Oh no, not there. That's continuous, sorry. Discontinuous Galerkin would use half hats inside the interval. It would use that one. And then another one it would use would be one that jumped this way. In other words, it sticks within an interval, it uses two functions instead of one. There are two unknowns. Do you see what's happening now? My space of trial functions is going to have jumps. Up to now, all my good hat-- my ordinary hat functions were all continuous. They were all C zero. Now I'm dropping back to C minus 1. Somehow my functions have jumps. So this is called DG, discontinuous Galerkin.

And how could it be justified? I mean, right? These are functions-- Now I've got a function that's not allowed in my second order problem. Right? If I looked at my list here, what was OK and what was not OK, if I look at interior half hats, if I look at functions that have a jump, jump functions, shall I say, jump functions. Those would be C^(-1), jump functions. And phi' is not OK. So I won't tell you, and I wouldn't be able to do it well, how do they rescue that? However Professor Peraire and others, others, Professor Darmofal, a whole group in Aero, computing with functions that have jumps. And the way they manage to do it is, they have to add some penalty term that penalizes the jumps, that controls the jumps. If they just took the standard set up, as we have, and tried their functions, it would fail completely. So they adjust the energy to penalize jumps. Anyway, somebody might sometime end up working on some project with them, and DG would be very possible.

So you didn't ask me about any homework problems yet. Last chance. Anything to ask about? You got a lecture about crazy stuff today. Sorry about that, in a way, but I hope it was all right. So homework is now posted, with a MATLAB question using the elements that we've talked about and some questions from the text. And don't forget, then, my original comment that if we end up with some good codes for CSE that could be used by other people, I'm very grateful to get them. OK. Let's end today a little early, and I'll see you Friday for splines.