Lecture 18: Finite Elements in 1D (part 2)

Flash and JavaScript are required for this feature.

Download the video from iTunes U or the Internet Archive.

Instructor: Prof. Gilbert Strang

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

PROFESSOR STRANG: So. I got the preparation for finite elements. Again we're in one dimension, because that's where you can see first and most clearly how the system works. So the system was, really, to begin with the weak form that I introduced last time. The Galerkin idea that I introduced just at the very end of last time, and that idea is that instead of the continuous differential equations where-- Galerkin's idea is how do you make it discrete, and he'll choose some trial functions. Those are functions. And some test function. And I didn't get to say, but I'll say now, very often they're the same. So very often the phis are the same as the V's. And probably they will be in all my examples today. And then what's new today is, what choices would you make? You have pretty wide choice, but there are some natural ones to start with. And so that's where we are today. And then also the other part of today is how do we get from all that preparation to the equations that we actually solve. The KU=F, where does the K come from, where does the F come from? The F is going to come, of course, somehow from this right-hand side and the K is going to come from the left side. OK, so that's the overall direction if you want to know how finite elements work, that's the key line there.

And then here I've recalled what the step was, here's our differential equation. And there's its weak form. This is the weak form. And let me put this in. Here I've reproduced what we reached last time, the weak form. That was the strong form with its boundary conditions; now we get the weak form with its boundary conditions. And the weak form involves u, so it's the boundary conditions on u, the fixed ones, that get used in the weak form. The free ones don't appear in the weak form, but by the magic of integration by parts, the free ones will sort of come out in a natural way. But we have to build in the essential, Dirichlet boundary conditions like that fixed right-hand end. OK, and because I think of v as a little movement away from u, but my u's are all fixed at the end, therefore v has to be fixed, too. OK, so this is important. But we haven't made it down to Earth yet. We were doing a lot of ideas last time, which are the center of things, but now let's get down to Earth. And down to Earth really means, what functions do we choose? And then how do we get the equation. So that's the job today, the principal job and this is of course what you would actually code up to run a finite element simulation. You would make a decision on these functions.

And let's start with piecewise linear functions. So a typical phi is maybe centered at a node, so it goes up and down. So if that's node two, let's number these nodes: one, two, three, four, and five. Well, that's five there but I'm going to have, let's see. What do I have here? At a typical node, two, my function is zero except in the intervals that touch node two. So here's phi_2. phi_2(x). That's its graph. It comes along, it's piecewise linear up, it's piecewise linear back down, and it's zero again. So that would be phi_2. Then phi_1, there'll be a phi_1. Like so, exactly similar. The whole point is keep the system simple. That's what the finite element idea is. Use simple functions for the phis. And I'm taking them also to be the test function v. Keep those simple. Now, what about at the boundaries? OK, well we've said our functions have to be, I'm doing free-fixed. So fixed comes to zero there and all my functions will do that. But this end is free, so watch this. This is another, what I'll call a half-hat. If I call those hat functions is that OK? It makes a nice short word. So these are hat functions. And then a fuller description would be piecewise linear, but hat functions is clear.

OK, now notice I'm sticking in here what I maybe could call a half-hat. Because my V's and my phis, my functions are not constrained at that left-hand end. That's free. So there's a phi_1, and there's a phi_0. And a phi_3, and a phi_4. So altogether, I'm going to have five. And I've started the numbering at zero, I guess it just happens. So let's accept that. So I'm going to have five functions. And they will also be my five test functions. So let me first think of them as the trial functions. Phi. So what was the point about trial functions? The point about is, my approximate, my finite element solution U(x) is going to be a combination of those. Some U_0, times the phi_0 function, up to U_4 times the phi_4 function. And these are my unknowns. U_0, U_1, those numbers. I have five unknowns. And I'm using hat functions. So I should, really, why don't I draw a graph of U(x). So I don't need the words piecewise linear. You see, I have a kind of space of functions. My functions are all, my approximations are combinations of these fixed basis functions. You could call these phis, trial functions, basis functions, you're always choosing, in applied math, a bunch of functions whose combinations are going to be, is what you're going to work with.

And here I'm making them hat functions. OK, so what does a combination of those guys look like? I'll even erase the word hat functions; we won't forget that. So what would a combination of-- Here's the same interval zero, one, with the same mesh, one, two, three, four guys inside. What would-- It just helps the visualization to see. If I combine these guys, so those were, despite the way it might look, that's one separate function. Two, three, four and five. And now suppose I take a combination. What kind of a function do I get? How would you describe this, a combination like this? Of those five guys? What will it look like? Between every node it will be a straight line, right? Because all these guys are straight, between those. So a typical function will start at what? This height will be U_0, because, well let me draw a few things. OK, so I'll put U_0 a little higher because I want to end up down at zero. So that might be the height U_0, and U_1 might be pretty close. U_2 might be coming down a bit. Coming down a bit, coming down a bit. And ending at zero. So those were meant to be corners. That wasn't a very good bit. There, OK. It wasn't meant to be, yeah. So that height is U_0. Now, notice why? Why is that? Why is the coefficient of phi_0 exactly the height, the displacement, at zero?

Because all the other phis are zero. All the other phis are zero at this point. Only this guy's coming in. So we'll only see in this, at this point, see, even phi_1 has dropped to zero. So we're only going to see phi_0 times U_0. So, I should have said, we take all these heights to be one. Height is one. Course it doesn't matter because we're multiplying by u's, but so let's settle. They all have heights one. And this is sort of a key part of the finite element idea. You see, Galerkin, when he created just two or three functions phi, he tried to follow the exact solution. He had an idea in his mind what the solution to the problem will be. And he wanted to take two or three functions that would get him close to it. Here we choose the functions, they're in the finite element library before we even know what the equation is, or the boundary conditions. These functions are like, the hat function choice. And they have this beautiful sort of a connect to nodes. In a way where-- In fact, I got involved in finite elements in the first place just to understand what's the difference between finite elements and finite differences.

Because the finite elements, as you'll see, are associated with nodes. phi_1 is the only guy that's not zero at node one. And therefore, what will this height be? What's that height? So that maybe comes down a little. What's this height? U_1. It's the coefficient of phi_1, because phi_1 is sitting there. And all the other phis are zero there, so this height is really U_1. And U_2, and U_3, and U_4, and U_5 is zero. So that was U_1, U_2, and so on. What I'm saying is, and we'll see it happen, is that this KU=F equation that we finally get to is going to look very like a finite difference equation. But it's coming from this different direction. And this way allows many more possibilities, gets things sort of more naturally right. It takes less-- With the finite difference equation, we had to go in there and decide what it should be. With finite elements, our decision is just the phis. Once we've decided the phis, Galerkin tells us what the equation is. And we'll get to the KU=F, what it is. But here I'm getting you to see what our approximate solution can look like. And this beautiful fact that the coefficients in here have a physical meaning. They're actually the displacements at the nodes for these simple functions.

OK. You got a picture of what the trial functions are, some people would think about the functions as these guys. Other people might think of this picture, the combinations. So those are the individual basis functions. That's a typical combination of the typical one. OK, so we're looking for an equation for these U's. Five equations, of course. Because we've got five U's. So that's my final step here. What are the five equations for the five U's. And those are the equations that I'm going to call KU=F. OK. So here's now a critical moment. Where do the equations come from? Well, the equations come from the weak form. So I take the weak form. And in for u, for u here, I guess it's just there. I put this. I put capital U. So can I just copy, this is the weak form for before we've made it discrete. Before we've chosen n phis. OK, now let's choose the n phis, so now this'll be the weak form. Again, the weak form with Galerkin. After the decision. So it'll be the integral, from zero to one, of this c(x), whatever it is. Times the U(x)-- sorry, times the dU/dx, right? dU/dx, so what is dU/dx? Oh, you have to pay attention. This was a true solution. Little u. But now this is where I'm working. I'm working with capital U. So instead of the d little u dx, it's d capital U dx. Maybe I'll put it in there. d capital U dx.

And then I'll put down here what it is. What is it? It's U_0*phi_0, can I use prime just to-- Or d phi_0/dx, whatever? Can I use prime for derivative? Just save a little space. So this is the derivative of my guy. U_1*phi_1'(x), up to whatever it was. U_4*phi_4'(x). That's what that term is. And that multiplies dv/dx. Where v is, here v is any test function. Any test function, only required to have v(1)=0. But now, I'm going discrete. So instead of any test function, I'll use these five functions. So I've got five functions. The phis are the same as the V's, then. V_1, V_2, V_3 and V_4. Same guys. So now I'll put in dV, can I say dV_i/dx? And on the right hand side I have the integral from zero to one of f(x)*V_i(x)dx, i is zero, one, two, three, or four. I'm testing against five V's. So I have this equation for five different V's. So i equals zero, one, two, three, four, gives me my five equations. Here are my five unknowns. This is my five by five system.

Let me just step back a minute so you see what happened there. So what do you have to do? You chose the basis functions, the phis and the V's. Then you just plug into the weak form, you plug in dU/dx, is coming from there. So this is dU/dx. You have to plug dV/dx, you have to do the integrals. You have to do the integrals, that's something we didn't have in finite differences. Finite elements involves doing the integrals, left side and right-hand side. OK, so, and here we have five different integrals to do. We have f(x) times each V. This will be, this number will be F_i. So my F vector is going to be an F_0, F_1, down to F_4. The five guys that I get from these five integrals. Alright. And the K matrix is sitting here somewhere. That's the last thing, that's the final thing, is to see what is the K matrix. Which is coming. This is somehow K's times U's are sitting here. F's are sitting over there. So it may be good to see the F's first. So do you see this now? We made the choices, then what's our job? Our next job is to do all the integrations. Integrate my function against V.

Let's make the natural first example, let f be one. First example, let f be one. All right. So if f is one, I'm going to find the system KU=F. So if I know f(x) is one, then I have everything I need to find these numbers. OK, actually we can probably do those by, you can probably tell me what they are. If f(x) is one, now. So this is example one. f(x) is a constant. One we have solved before. What's the integral of V_0? Right, that's what I have to do. What's the integral of this, f(x) being one, I'm just asking, I just have V_0(x)? The integral of V_0(x), and let me again draw V_0, which is phi_0. It's a half hat and then it goes along at zero. What's the integral of that function? One, yeah. How do I think about that integral? It's the area of the triangle. It's the area. That's what an integral is, it's the area. So the area is, I've got delta x there. Right, delta x as the base. One as the height, and you see the formula for the area of a triangle, it's got a half in there somewhere, right? A half. OK, can I factor out the delta x? Because the delta x is going to come in. I think there's a half there. And then what about F_1? What's the integral of this times V_1(x), the next V?

It's the area under this dashed function. Which is now the basis 2 delta x, so I get a one. Is that right? I get a one, one, one, one. OK, so that was obviously not too tough, right? That was straightforward and notice something. Even here, in fact, we see it here. I don't know if you remember about that half. Do you remember something about when we did finite differences and we had a free boundary? And we lost an order of accuracy if we didn't do it right? Do you remember that? At the fixed boundary we were fine, but with finite differences at a free boundary, where I was using the matrix T. With one, minus one at the top row, I lost an order of accuracy. Unless I made some change on the right-hand side. Look what's happening. The finite element method is making the change for me on the right hand side. So the finite element method is going to automatically keep the second order accuracy. Keep the second order accuracy. So that's a key point. That these piecewise linear functions are associated with second order accuracy. Later we'll move up to parabolas, to cubics; that will move up the order of accuracy in a nice way. Where with finite differences we would have had to create new finite difference formulas. Our minus one, two, minus one formula, that was good for second order accuracy. Then we would have to figure out, in the quiz had partly started it, what if there's a c(x) in there, what do you do? Finite differences, there's more thinking involved. Finite elements is like just press the button. Well, there's a little more to it than that of course, because it's taking a full lecture. But in the end it's more systematic, you could say.

So that's the F. Now are you ready for the K? So this is the key part, OK? So you have to can get this thing to simplify. So what am I looking for here? This whole left-hand side should be K times U. So I'm looking to see what multiplies-- I'm looking to make sense out of this. What's the first equation? Right, so the first equation or the zeroth equation, I guess. The zeroth equation, the one that'll run along and have this right-hand side, the zeroth equation is the equation when i is zero. It's the equation that comes from testing our weak form for V_0. For that particular form. Maybe I'll just start over on this board. Then I can write a formula, but I'd rather you see how it comes. So I'm looking at equation zero. So take i=0. So I have my left side is my integral, of c(x). Times this combination that I wrote, U_0*phi_0' plus... U_4*phi_4', times dV_0/dx*dx equal, and on the right side is where I got the F_0, which I already figured out to be delta x times a half.

It's the left side that I'm worrying about. OK, you see what's happening here? This is some matrix. Its zeroth row is what we're finding -- multiplying U_0, U_1, U_2, U_3, and U_4 -- equaling the F vector. I'm supposed to be getting the first row of the matrix, the top row of the matrix, from the top V. OK, so let's just do these. We've got integrals to do again. Alright, what is dV_0/dx? Do you see? Let me just see? What number is going in here? What number is going in there? Yeah, if we see that we're golden. What number is going in there? That's the thing that multiplies U_0 in the first row that means I should use V_0, so this is the point, this is K_00. And what's its formula? You realize I'm starting the count at zero because all these counts-- So what is K_00? It's an integral. Of what? c(x), good. Times this guy, because it's multiplying, times this guy, V_0', dx. That's what you have to do. That's what you have to do. c(x) times phi', it's phi_0', that's what would sit there. And maybe, well, let's figure that one, shall we?

I have to know c(x), right, that's part of the problem. What would you like me to choose for c(x)? One. Thank you. I'll choose one. Let this be one. Or I could make it capital C and you would see a capital C appearing everywhere, but let's make it one. So what are we doing now? What's our equation? Our right-hand side is one, our c(x) is one, our equation has reduced to -u''=1. The first equation in the course. So we're back to September the 3rd or whatever it was. But doing it now by finite elements. OK, so let c(x) be one and tell me what this integral is. So c(x) is now, we're taking-- In our problem, we're supposing it's one. Let me just say: Suppose it's function. Then we have lots of integrals to do involving that function. And we might not do them exactly, that would be alright. It's certainly totally OK to do the integrals approximately, because we're doing everything else approximately. So we just have to be sure that we do the integrals with sufficient accuracy so that we don't lose accuracy in the integrals. Of course, with a one we're going to do the integral exactly. But if c(x) was some variable function, I wouldn't have to do it exactly, I would just have to do it with enough accuracy so that I don't lose extra accuracy beyond what I'm losing in the whole Galerkin approximation.

OK, ready for that number. What number comes out of that? phi_0', let's graph phi_0'. And of course it's the same as V_0', so can I put a little graph here? Here is zero to one, and I'm going to graph phi_0'. So what's phi_0'? Oh, it's negative, isn't it? My little graph isn't going to work. I didn't leave enough room. It's negative. So I'll just write it in words. It's the same as V_0', and what is it? Tell me what it is, what's the derivative of that function? It's what? Negative one. Wait a minute. Yeah, it's going to have a certain value, yeah. You can tell me what it is beyond that point real fast. So it's something up to delta-- So what is the slope? What's that slope there, of phi_0? It's not negative one, because remember, what's the base here? That's not the point one, I'm sorry. All these were delta x's. Those were just numbering the nodes. But the actual length scale is the delta x scale. So now tell me what it is. The derivative is negative one over delta x, right? It dropped by one when it went across by delta x. And this is only up to node one. Up to delta x and then zero afterwards. This is a key point. That all our functions are local. Our functions are local. What does that mean? You can tell me what am I going to get when I integrate, for example, when later on I might be integrating phi_1' against V_4'. What's the answer? This is the key point. Later, when I'm looking for the one, four entry, when I'm looking there, I'm going to do an integral of phi_1'. I'll erase for a moment and do this in my head. When I integrate phi_1' against V_4'. Maybe it's the fourth row. And the first guy over, maybe it's this guy I'm doing. Doesn't matter a whole lot. Because the answer is, when I integrate phi_1' against V_4' just, it's nice to get the easy ones. It's zero. Why is it zero? Why is the integral of phi_1' against V_4' zero? Because these phis are local. phi_1' is only non-zero here. V_4' is only non-zero over here. The two don't overlap. Anywhere the one is not zero, the other is zero. So that's a zero there. In fact, our overlaps, I'm just sort of looking ahead here. Our overlaps, a phi overlaps itself, of course. And its right-hand neighbor and its left-hand neighbor. But nobody two or three or more away. I think our K, all our integrals are going to be zero outside, we'll have another tri-diagonal matrix. We're going to have zeroes all here. And we'll only have entries of phi against V when they're either the same or just differ by one. So we'll only have three diagonals.

OK, we were about to find out what that number is. So the slope of this is minus one over delta x, and that's-- I'm sorry, let me go back to zero, zero. OK. This is what we're keeping our fingers crossed for. What's that number? So I have this thing, actually is it just squared? And that's the slope. And then the phi and the V I'm choosing the same, so that's the slope again. I think I'm just getting one over delta x squared for that times that times the one. So what's K_00? One over delta x. Where'd the delta x come from? Because we're only integrating over, it looks like zero to one but they're all zero. We're really only integrating, the only reality was out to node one. Out to delta x. You see that the number there, the number here on the diagonal is one over delta x.

OK, how about doing K_11 for me? So again, now these guys will be the same guys. It's a square. No it's the integral phi_1' against V_1', they're the same. And what is phi_1', which is the same as V_1'? What's the derivative now? It's, ah. What's the slope of this function? It goes up and goes back down, right? I have a plus part, so the slope going up is the one over delta x. And then the slope coming down is minus one over delta x. So this was up to delta x and then to two delta x, and then zero. That's a much more typical thing, the slope goes the function, the hat function goes up to the top of the hat. Back down. The slope up and the slope down are easy. And now the integral's easy. So I'm just squaring, well, when I square it, this squared is the one over delta x squared. This is the same, because the minus will get squared. So what's K_11? What's K_11 now? Have you got K_11 in your head? This is one over delta x squared. And now what is the integral? Two over delta x. Because now we're integrating from zero to two delta x, because that's where my functions are going out from zero to node two. If the function's numbered one. So it's two delta x times this; I think we get a two over delta x on that diagonal.

Would you care to guess the rest of the diagonal? Yes, you tell me what's K_22 and K_33 and K_44? They're all the same. We're just shifting over. So two over delta x goes down there. Alright, one more to do. One more integral to do. This next guy. So can you tell me what do I get now for K_01? K_01, so now this is the case where I'm in row zero, so this should be V_0, because that tells me the row I'm in. But phi_1. What happens when I integrate, just see the picture here. Let me just draw it small. phi_1', so let me draw phi_1. And V_0. OK. But it's the derivatives that I want. It's the slopes that I want, OK? So what do I get from here on out? Zero, because this guy only got to there. That's the half hat, the first guy stopped at delta x. So whatever is happening here is going to be multiplied by zero. So it's just here. One delta x interval for this one. They just overlap in one interval, of course. This guy and its neighbor only overlap in one interval. And what's the deal about the two slopes? They're opposite. One's coming down, one's going up. But the slopes are one over delta x and minus one over delta x. Do you see what's happening here? I'm integrating. Here I have a slope of one over delta x, and here I have a slope of minus one over delta x, so I should multiply those. Minus one over delta x squared, integrate, what goes in K_01? What's that number? It's that times that integrated, but now the integral is only going really out to delta x because basically I'm just stopping there.

So-- But there's a minus now. Because it's not the square. It's this times its neighbor. One's going up, and one's going down. So it's delta x squared, and then the length of the interval, this is a minus one over delta x. Would you care to guess the rest of this matrix? What's the rest of that diagonal, above the main diagonal? It's all the same. That stays the same, because when I do phi_2*V_1, my picture is just like shifted over. But I still have one coming down, and one going up. When I do phi_3*V_2, same thing. And if I do phi_1*V_2, it'll be the same. I'm going to get this minus one over delta x all the way on that diagonal, also. Symmetry. It's going to come out symmetric. Actually, since the course started by speaking about properties of matrices, let me just say K is going to turn out to be symmetric positive definite. And what's more, for this example we recognize K completely. You will say: Why did you take so long to get to this result? K is the one over delta x part times, what's the matrix? It's T. It's T. So it's one, minus one, minus one, two, minus one, minus one, two, minus one, minus one, two, minus one, and I guess we had five of them. Oh, but nobody there, right? That's not there. That fixed-- Why do we not have a minus one? Because we've got no, there isn't a six. That would be column-- We've got one, two, three, four, five columns; there's no U_5, there's no V_5, we've got them all.

And the F, so that-- So KU, this thing multiplies U_0 to U_4, to U_4, and it produces F, which is one over delta, which is what? Oh, delta x is in the numerator, right. Times a half, one, one, one and one. That is the finite element system KU=F. For this simple problem. It's exactly what finite differences did. So you can see why my first introduction to finite elements was with the question: what's the difference? The finite element community at that point, this was like the golden age of finite elements, all this was just beginning to be created. These elements were being used. Especially in civil and structural engineering, that's where a lot of the earliest papers came out of. And then, in a model problem it didn't look anything new. It looked like our original finite difference matrix. But there were some new things. First, there was this new 1/2, that we hadn't particularly noticed with finite differences. We, we could catch onto that. Here's a minor difference. You notice that the delta x is-- Strictly speaking the delta x is up here then, but when I divide by delta x then I'm back to the finite difference. I have the one over delta x squared, it looks like finite differences again. So everything looks the same. But, of course, if c(x) isn't one or if f(x) isn't one, oh yeah. If c(x) isn't one then I've got integrals to do. I would approximate those. And I could then come out with something that would look like a finite differences.

Let me take our other favorite model problem. What would be the F if, yeah, here's a question. What would be the right side if my vector, instead of being one, what's my other favorite choice? Delta. So I take delta at x minus, let me take delta at x minus 1/4, first. Suppose that's my f. Then I've got to change all these guys. And what would they be? What would be the new right-hand side when I have this point load? I have to go back to the integrals, right? I have to go back to these guys. These integrals, with that new f, this is now delta of x minus a 1/4, times each V, times dx. I have to integrate delta of x minus 1/4 against every hat function. And see what it equals? And what will I get? You're going to tell me right away. What are those integrals? That's a point load at node one. Times the V integrated over the whole thing. What do I get? I get a one, yeah. That integral is going to pick out the value at a quarter. Right, that's what the delta function does, the spike is at a quarter. Has area one, so it picks out the V_i at a quarter. V_i at a quarter will be? One for the-- I think we get a [0, 1, 0, 0, 0].

Again a little bit what our finite differences suggested that we should do. Alright, here's one final one for today. Suppose the delta function is not at a node. Suppose it's at 3/8. Or it could be at any point a, but let me just take a typical, a special one where I can do it. Suppose the load is at 3/8. What do I get for the integrals now? So now, it's delta of x minus 3/8. The spike is at this point here. That's where delta is now, spiking at 3/8. Is that right? We had 1/5, tell me what-- oh, I should have had 1/5 before, sorry. Change that on the videotape. All those 1/4s were-- That 1/4 was 1/5, and now what do I want? Three? I wanted to take a nice one that was halfway. I just forgot what halfway was. Where is halfway there? 3/10 now for delta. So that was-- Before I had it for when delta. So previously was delta at x minus 1/5 and now delta at x minus 3/10, what's the F now? What's the F now? So the spike is right in the middle between one and two. What's do those integrals come out to be? If I integrate delta function times the different hats, what do I get? What do I get, yeah. I get a zero for this first guy because it didn't touch the half hat. And then what do I get there? Half. And what do I get at the next one? Half again. And then the other guys it doesn't touch. You see, it automatically does it. Does those smart things. It automatically makes the smart choice. And if the spike was at a, at any point a, then at that typical point a wherever it is, like there, spike could be there. Then I would have, what would I have if the spike was there?

I'd have a little bit of phi_3, and a big bit of phi_4. And the two parts would add to one. It would take the right proportion. It would be the proportion, by however much this spike was near there, it would give that extra weight to phi_4. OK. So there is the finite element method. It produced something that you might say, oh, we knew that. But you've got to see that it deals automatically with c(x), it deals automatically with f(x), it deals automatically with the free boundary. You see, the solution there is going to take sort of a balance, a pretty close balance, of this half hat with this one, and the solution will actually be a pretty close to free. It'll be pretty close to having the right zero slope there. OK, good.