Lecture 20: Element Matrices; 4th Order Bending Equations

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 education resources for free. To make a donation, or to view additional material from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

PROFESSOR STRANG: Okay. Hi. So, our goal, certainly, to reach next week is partial differential equations. Laplace's equation. Additional topics still to see clearly in 1-D. So one of those topics, these will both come today, one of those topics is the idea, still in the finite element world, of element matrices. So you remember, we saw those, that each bar in the truss could give a piece of A transpose A that could be stamped in, or assembled, to use the right word -- I think assembled is maybe used more than stamped in, but both okay -- into K. For graphs, an edge in the graph gave us a little [1, -1; -1, 1] matrix that could be stamped in. And now we want to see, how does that process work for finite elements. Because that's how finite element matrices are really put together, out of these element matrices. Okay, so that's step one, that's half today's lecture.

Step two, problem two now, coming from the next section, is fourth order equations. Up to now, all our differential equations have been second order. Are there fourth order equations that are important? Yes, there are. For beam bending. So I'll describe that application, which leads to fourth order equations. Do they fit our A transpose C A framework? You bet. You know they will. Each additional application in the framework kind of gets us comfortable, familiar with that framework, what it can do. The A transpose C A.

So let me start with element matrices. I'll get the homework, those numbers will get posted on the website later today. I just thought I'd put them down, and I have to figure out what would be a suitable MATLAB question. So my idea for this homework, as it really was for the last homework, was that you get some, not a large number of ordinary questions, paper and pencil questions, and one MATLAB question. When I wrote "MATLAB last", people interpreted that as "last MATLAB". But those words are not commutative. Right? I just meant -- and it doesn't really matter, it was a dumb thing to say -- I meant you to put the MATLAB question at the end after the regular questions. Sorry. So that MATLAB is due-- Maybe a bunch of those turned in on Monday didn't include the MATLAB. No penalty. I'm talking now about the MATLAB for the trusses. How many have still got a MATLAB for the truss still to turn in? A number. Oh, not too many. Okay. Anyway, just when you can. In my envelope is good. Okay. So it'll be a similar thing for this week, and because it's short I said okay, let's get that cleaned up by Monday. And then we're ready for vector calculus, partial differential equations, two and three dimensions, the next big step.

If I want to illustrate element matrices, the best place I could do it would be to go back to our piecewise linear elements in 1-D, and see how an element now is just one of those little intervals. Little pieces of the whole structure. If it's a bar, or whatever it is, I've cut it into pieces. Here's one piece, here's another piece. With finite differences, if I gave you unequally spaced meshes, if my little h is different from my big H, we'd have to think again. I mean there would be some three-point-- Instead of minus one, two, minus one for the second difference, it would be a little lopsided, of course, when the mesh is unequal. We would have to think that through for finite differences; for finite elements, the system thinks for us. So I just want to show how that happens.

So I'll take this. What's the element matrix for our standard equation for that element? Okay. And then it would fit into the big matrix K. I'm going to come up with a matrix K that I could do in any other way, but this is the way it's really done. So I look in that interval. I'm focusing on that interval. So here I've drawn the basis functions. But now I'm not going to operate so much with basis functions. Let me just think about that interval again. That's the same interval. So those were the basis functions that went up to height one. One started at height one and went down, one started at height zero and went up. But it's their combination, of course. Let me number this node number zero and this node number one.

So we have some height, U_0, here. This is coming from U_0 times that hat. And some other height, U_1, here, which is coming from U_1 times that hat. But inside this interval, my U in this interval is U_0 times that first hat function, the phi_0, the coming down hat. Plus U_1 times the phi_1 function, the going up hat. It's a linear function, so it's going to be there. That's my graph-- No. Yes. This is U(x). U(x). Okay. Right? I'm focusing inside one element. One interval. I've numbered them zero to one, but of course that distance is H.

Okay, so there's my function. Right? Everybody agrees that this combination, it starts out right, it ends up right, at height U_1, and it's linear. So that's got to be right. So really, what's the contribution from that element? I look at the quantity, that-- And I'm always taking the phi functions to equal the V functions. So I can write this as c(x) times dU/dx squared, dx. Over the whole thing, over the whole interval, that would be the U transpose K U. This, everybody remembers the K part, of course, comes from the left side of the problem. It comes from integrations. And this d-- U is now a combination of all the-- With weights U_0, U_1, U_2, U_3. If I plug in for capital U and do all the integrals, I'll see what K is.

Now what's the point? The point is, when I did those integrals, the question is, how am I going to do the integrals? Do I do them-- The way I did before, was I watched what phi_0, one phi times another phi and I integrated. That was successful. But the new way is, integrate it an interval at a time. Do the integrals. So this would be K_global. Now I'm going to go just from zero to H. If I go from just zero to H, that's going to give me the K_element piece. The little piece that comes from this little interval. And on that little interval, this is my formula for U. I'm just hoping you'll sort of see this as a reasonable idea, and then when we do the integral you'll see it. It clicks.

Okay, how does it click? So what's my plan? My plan is, here's my function. There is its picture. I'm just going to plug it into here. Oh, it's going to be simple, isn't it? What's the slope of that function in that element? On that interval? What's the dU/dx? So instead of doing every full integral for each separate phi, I'm doing every element integral for both phis. See, these two phis are both coming in to that element. Okay, so what's dU/dx.

I didn't realize how neat this that's going to be. So what's dU/dx in this element? I'll use a different board here. So it's an integral, then, from zero to H of whatever my c(x) might be. And I'll make a comment on that. But my focus here is, what's dU/dx? What's the derivative? And then I'm going to square it. For this function, for that picture, the slope is obviously (U_1-U_0)/H. Right? The slope is (U_1-U_0)/H. And I'm squaring it. dx.

Okay. Let me take c(x)=1. Just to see clearly what's going on. So c(x) is just going to be one. Okay, so I'm claiming that this is my U transpose. My little piece. Why is it only a little piece? Because it only involves two of the U's. It's going to be a little two by two element matrix, that comes from this element. And then it's going to to be put into the big K, the global K, in its proper place. Okay, well. It's trivial, right? This is a constant, this is a constant, the integral is just H times U_1-U_0, squared, over H squared. Because that's getting squared, the interval was length H. I think I just have a 1/H. So this is my U transpose K element U. And now I want to pick out, what's that matrix? What's that little two by two matrix that only touches these two U's? What's the matrix-- Well, since it only touches these two U's, you can tell me. I want a [U_0, U_1].

Sorry, let me make a little more space. And you can tell me what matrix goes in there. I want this all to match up. [U_0, U_1]. Now here is the two by two element matrix. [U_0, U_1]. What's the two by two matrix that will correctly produce this answer? I'm just shooting for that answer. What do I have here? This is U_1 squared minus 2U_0*U_1, plus a U_0 squared. Right? And I have to remember the 1/H, so it's automatically going to come out right. So there's a 1/H, shall I remember that first off. 1/H is part of my element matrix. And then, what are the numbers that go inside that matrix?

We had practice with this. You remember when we talked about positive definite matrices, way back in chapter one? The point was that we could look at eigenvalues, or pivots, or something. But the core idea was energy. The core idea was that energy, that quadratic, and that's what we're looking at again. That's the energy. That's the energy, right there, and this is the energy, this is the energy in the finite element subspace. All I'm saying is, what matrix, what two by two matrix, goes with U_1 squared minus 2U_0*U_1 plus U_0 squared. Just tell me what to put in that matrix.

What do I put in here? One, right. Because it's multiplying U_0*U_0. What do I put there? Minus one. Good. Because I have a minus two, it comes in, minus one comes in twice, and a one goes there. So there, with the 1/H included, is K_e. K_element, for that element. For the big H element. You'll say big deal, because we've seen this thing before. Notice what is nice. First of all, notice how nice that is. It's particularly nice, of course, because I took c(x) to be one.

So let me make a comment. Suppose c(x) was not one. What would I do? Suppose c(x), suppose I have some variable stiffness in the material. Suppose the material could be changing width, so its stiffness would change. So in other words, I'd have a variable, c(x), that I should do the integral. Probably finite element systems aren't set up to actually do the exact integral. What would they do? They would take, for that simple integration, they would probably just take c(x) at the midpoint. So there's a numerical integration here in the creation of these element matrices. And numerical immigration is just, take a suitable combination of the values at a few points. You do know Simpson's rule? Simpson's rule, that's a pretty high level rule. My suggestion there was just a midpoint rule. Just take c(x) at the middle of the interval. Then it would factor out.

So I should really put a c here. A c should really be coming in there. And you expected that, right, from the A transpose C A. We always saw a C in the middle. It really should be there. When I took c to be 1, I didn't see it.

So what am I doing? I'm approximating c(x) by c at halfway. Approximately. I would replace this unpleasant, possibly varying function by c at a point and use that value. So numerical integration is one part of the picture that-- We won't go into all the different rules. There's a rectangle rule, there's a trapezoid rule, very good. There's the Simpson's rule that's better. As I get higher order elements, like those cubic elements I spoke about, the numerical integration has to keep up. If I was integrating cubic stuff, I wouldn't use such a cheap rule. I would go up to Simpson's rule, or Gauss's rule, or somebody's rule.

Anyway, that's the c part. Here's the part that stamps into the matrix. Notice, by the way, when I stamp it in, tell me how it's going to look stamped in. And then I've completed it. So here's my big K. I wish I had a little more room for it. Okay, here's my big K. So that interval that I drew there, the H interval, will stamp in here, some two by two. Right? Now where will the similar thing coming from this guy. So, maybe it has to be numbered minus one. Sorry about that. That's another little interval. I'll do the same thing on that interval. I'll get a little element matrix, two by two guy, for K for that element.

And where will it fit in to the big K? Does it fit in up here, let me just ask you. Does it fit in up there? Yes, no? I'm assembling, stamping in the small two by twos into the full n by n. And if I draw the picture you'll see it. So when I do the two by two for that big H interval, it goes there, let's say, then I just want to say where does the two by two go for the interval to the left? Does it go there? Nope. How does it go? It overlaps. Right? It overlaps.

Why does it overlap? Because the phi_0, this guy, is acting on the right, and also acting on the left. The U_0 is active, is partly controlling the slope this way, and also that way. The U_0 is in common the two intervals. Anytime any unknowns, any mesh points that are in common to two elements, we're going to have an overlap when we assemble. And so it'll just sit, it'll sit right there. And so there is the diagonal guy. And maybe you could tell me what number will go there. What number would actually go there? And then you'll see the whole point of assembly, stamping in.

Well, what number goes there from this? One times the c/H. So in here would be the c -- can I call it c_right, or c_H, c on the big H interval -- divided by the H. Because that one, we had to get it right. And then what will go in that very same spot? So add it in coming from the small h interval. The same thing, it'll be this one, like shifted up, moved over. So it'll be coming from that one, so it'll be a c on the little h interval, divided by little h. That would be the diagonal entry of K. That's what we would see right there.

Over here, should I write a typical row of K? Typical row of K, when I do that, is going to have this c_h/h plus c_H/H. That's like the two, right? That's like the two. And what goes here? What will the entry be that sits there when I assemble? Just this guy. Right? Just this guy, times that. The entry here will be the minus c_H/H. That'll be the entry. This is the diagonal one, this is the one to the right, and what's the one to the left? What's the one here? Well, you know what it is. It's going to come from the minus and it'll be the minus c little h over little h.

Look at that. That just shows you how it works. And again, you can look at that, page 299 to 300 in the book. You see that if the c's are the same, if the h's are the same, then we're looking again at our minus one, two, minus one. Times whatever c over h, to keep it dimensionally right. Do you see that? It's just simple. Simple idea. The point is that each interval can be done separately. It's a simple idea in 1-D. It's a key idea in 2-D, where we have triangles, we have tetrahedra, tets. We'll see this in two dimensions, later in this chapter, when we're doing Laplace's equation. It's just fun to see it work. You'll have different triangles, say call them triangles.

So that phi-- Do you want me to look ahead? Just ten seconds, to triangles? So imagine we have triangles here, so we have piecewise, we have little pyramids. Instead of hat functions, they grow to pyramids. So there's a pyramid guy whose height is one there, and drops to zero in all these places. And that's our phi. Our trial function, test function, will be a pyramid function, then. And I can do integrals that way, or I can take the integral over a typical triangle. So a typical triangle is involved with three, now I've three functions, in the linear case, controlling inside that triangle. So what will be the size of the element matrix? Can you sort of see how the system is going to work? And then we'll make it work in 2-D. Every node, every mesh point, corresponds, has a pyramid function, has a U that goes with it. Those U's are the unknowns. And how many of those unknowns are operating inside that triangle? Three. So what will be the size of the element matrix, the non-zero part of the element matrix? Three by three. What else could it be? So we'll see what it looks like. And we'll have integrals over triangles.

So that's good. Okay, thanks. Exactly halfway through the hour is exactly that first topic of element matrices. Done. Okay. Let me take two deep breaths, and move to fourth order equations.

Fourth order equations. For the bending of a beam. So I'd better draw a beam. This is a 1-D problem, still. This is a 1-D problem, still. To keep it 1-D, this better be a thin beam. So this is a thin beam. And the loads, what's the difference? What's here? It looks like a bar, pretty much, right? But the difference is, the load is acting this way. The load is acting that way on the beam. Maybe two loads. Maybe a uniform load. Maybe the weight of the beam. But it's transverse. It's in the perpendicular, it's transverse to the beam. It's this way. So the beam bends. Let me do a fixed-free. So this'll be a fixed-free beam. Fixed-free, the word for fixed-free would be cantilever beam.

Okay. So what happens if I impose those loads. Well, the beam bends. So the displacement is now downwards, is now not the direction of the rod. The displacement I'm interested in is perpendicular to the beam. Downwards.

Okay, so we can start on our framework. So this is displacement. u(x). I'll stay with the same letters. So you know I'm going to have an A, that will take me to whatever this is going to get called. What's the quantity there? It's going to be, let's see. What happens? So this is just geometric now. Let me put in the easy part, C. That'll be sort of the bending stiffness. Right? This'll be the bending stiffness. Because the beam is going to bend. And over here I'll get a suitable w. What I get here-- So when the beam bends... Its curvature. Curvature of the beam is what's produced. It's not stretching of the beam; it's curving of the beam. So this quantity, e, will be the curvature that comes from the displacement. If I displace these beams, suppose I put a load here, it's going to bend that down. The bar will curve. So the curvature e. Now then, the question is, what is this A; what is the curvature?

Well, do you remember? You're on the ball if you remember the formula for curvature. It's a horrible formula, actually. But that's only-- We're going to make it better. Do you remember the curvature, it was in calculus. Yes, you all remember this. Suppose I have a graph. I know its slope, that's become easy now, right? Calculus. But the curvature of it, what derivative did it involve? Second derivative. And was it the second derivative exactly? No, unfortunately there's some term which is horrible. One plus the first derivative squared, all square root. But I'm just going to take u double prime. So this is an approximation.

So what is A now? What's my matrix A, that gets me from u-- or my operator A, that gets me from u to e? What's the e=Au equation? A is just second derivative. That's something new. A is second derivative. And why do I do that? Because I assume small curvature, small displacements. I assume that u' squared is very small compared to one. So it's just slightly bent. A beam that goes way down here, I'd have to go nonlinear. But if I want to keep things linear, I approximate. This will be much smaller than this, so one is fine, that term goes.

And now the next step will be easy. What's the bending moment? This will be called the bending moment. And let me use the letter w again. I should really use capital M for bending moment. That will be the stiffness times the curvature. So that's the force, the way the spring had a restoring force by Hooke's law. This is the equivalent of Hooke's law. But this restoring force is not pulling back, it's bending back. It's torquing back.

And then, of course you know, there'll be an equilibrium equation to balance the load. So this load is the f(x). And you know what that'll be. Because you know it'll involve A transpose. And the transpose of that, do you want to just make a guess? I shouldn't use transpose, of course, but I'll use it again. It'll be the same; it'll be second derivative. I mentioned we had a minus sign with the first derivative, but now we're going to have two minus signs, so it's going to come out symmetric, second derivative, and the equation here will be w''=f(x).

Our framework is working. We have plus boundary conditions on u. And here we have plus boundary conditions on w. So those parts we have not yet mentioned. And of course, that depends on my picture. While we're at it, why don't we figure out, what do you think is the boundary condition here? And how many boundary conditions am I going to look for altogether? Four altogether. Because I have a fourth order equation. There will be four arbitrary constants until I plug in boundary conditions. So I'm looking for four boundary conditions, two at this end, two at that end now. And what will be the two at the fixed end? At the fixed end, obviously, it's built in. Built in. Slightly different words sometimes for the beam problem. Here I'll have u=0 and u'=0. Those apply with A; those go with A. Those are the essential conditions, the Dirichlet conditions, the ones I must impose all the way. And now at this free end, what do you think?

Well, it's great. It's w=0, and w'=0. It's just beautiful, the way it all works. So that's a completely fixed-free. Then why don't I draw in, just while we're talking about boundary conditions, an alternative. So here's my beam. And now you see, it's under a load, f(x), transverse load. Now that would be different boundary conditions. Anybody know the name of a beam that's set up like that? Simply supported. You don't need beam theory, and I don't know beam theory, to tell the truth, to do these problems. So that's a simply supported beam.

And what are the boundary conditions that go with that? Well this is u=0. No displacement, it sits there on that support. And what else is happening at that support? There's no bending moment. Nobody's here. Right? So it's w=0. And at this end, too. Also at this end, u at one is sitting there, and w(1) is zero. Yeah. That would be the boundary conditions, four of them, for a fourth order equation that we'll just write down in a minute, for simply supported.

And we could have a mix. This could be simply supported here, free here. I think. Or maybe, could it be, or maybe that's too risky. Would that be a singular case, simply supported? Huh. So as always with boundary conditions, some are unstable. Some are not going to determine all four constants. Just the way free-free didn't work, right? Free-free for a rod didn't determine anything, it left a whole rigid motion. Maybe u=0, w=0 at one end, and free at the other end, is that--? It sounds risky to me. But we can see.

Okay. So, do you get the general picture of the beam? So what's the equation? What's A transpose C A, when I put it all together? I'll use that space to put in A transpose C A. Continuous, we're talking here. Right now we've got differential equations. So what's the differential equation, A transpose C A? So it's the second derivative. I'm just going backwards around the framework, as always. The second derivative of this, and this is c(x) times e(x), and e(x) is second derivative of u, equal f(x).

Good. In ten minutes, we've written down the framework, some possible boundary conditions, and the combined A transpose C A equation. I mean, we're ready to go. We've got the pattern to think about this. So let's see, what should we do first? I would say the first thing to do is, let c be one and solve some problems. Let c be one, and consider-- So if c is one, it's a fourth derivative equation.

Should we take uniform load? Yeah. How does a beam bend under its own weight? So it's just one, or whatever constant. So it's constant load, it's just its own weight, it's going to sag a little in the middle. What's the solution to that equation? And what shall I take as boundary conditions?

Let me do the simply supported one. Because that would be kind of nice. So it's simply supported, it's sagging under its own weight, with u(0)=0, u''(0)=0, because that's the w. And u(1)=0, u''(1)=0. Whatever. I don't know that I'll have the patience to go through and plug in all four boundary conditions to determine all four constants. Just get me to that point. Get me to a solution u(x), the general solution here that's got four constants in it is what? Okay.

Okay, think again. What are we looking at? We're looking at a linear differential equation. Linear problem. I'm asking for the general solution to a linear equation. What's the general set up? General set up is, particular solution plus nullspace solution. Right? You see an equation like that, you're looking for the general solution, tell me one particular solution and then tell me all the solutions when it has zero on the right, and we've got everybody. So that was true for matrices, it was true for Ax=b, it's just as true for differential equations. So what's one particular solution? What's one function whose fourth derivative is one? Yes? What am I looking for here? 1/4 of x to the-- No, what? 1/24, is it? 1/24 of x to the fourth? x to the fourth over 24. Because four derivatives-- So we're thinking we're in the polynomial world here. Just as we were with u''. With the bar it was x squared over two, the particular solution, now we're up to x fourth over 24. So that's the fourth derivative is one, good. So I'm seeing a fourth degree bending there.

And now what about the null space solutions, the homogeneous solutions. This accounts for the one, now what are the possibilities if it was a zero? You're going to tell me the whole bunch, right? A plus Bx plus Cx squared plus Dx cubed. Because all of those have fourth derivatives equal zero. So that's the general solution. Okay. So whatever the boundary conditions are, they determine A, B, C, D. We're not that far away from Monday's lecture on fitting cubics. Actually we're really close to it. When we use finite elements, we're going to use exactly those cubics. I'll get to that point.

Let me take the other model problem, that everybody knows what's coming. What's the other right-hand side that this course lives and dies on-- lives on. Delta function. Right. Delta at some point. So that's a point load then. I can make it-- Whatever the boundary conditions are.

Right. Good point. This boring stuff will just repeat, right? That's the null space solution. But now, what is a particular solution? The particular solution has become interesting. The particular solution here was straightforward, simple, a good one to do first. What's a particular solution to a point load? So instead of having distributed load here, I'm putting a heavy weight here at this point, a. And it's heavy weight, I'm multiplying the delta function by one, I could multiply by some L for load or something, but let's just keep it simple.

What's a solution to that? Fourth derivative equals delta. So that means I've now got to integrate, one way to get the answer here would be to integrate four times. Right? If I integrate delta four times, then I've got something whose fourth derivative will match delta. So do you remember the integrals of delta? Okay, so I integrate. First integral is: step. Second integral is: ramp. Third integral is: quadratic, right? This was linear, boop boop, linear pieces. The next integral is going to bring me up to quadratic ramp. And the next, fourth one, is going to bring me up to cubic ramp. Cubic. So it's going to be cubic ramp, is what I get there. So one particular solution would be a function that's zero, and then at the point a, it suddenly goes up cubically. So it's zero here, and it's x cubed over six there, I think. If I integrate three times, I'll be up to x cubed over six. Is that right? Yes, cubic.

So that's an interesting function. Of course, these parts will tilt the function, will change it. So our solution won't look like this, because I've only got one particular function, and I'll need these to satisfy the boundary conditions. So there is one particular solution.

The general solution-- Yes, good for us to think out the general solution. What does that picture look like when I add in this stuff? Very, very important. I'm sorry I don't have more space for this highly important picture. Okay, so here's my point a. Keep your eye on that point a. Okay, so to the left of it, I've got some curve, whatever, dut dut dut dut dut, the beam. And to the right of it, I've got some other curve, whatever it is, dut dut dut dut dut. And what's cooking at point a? What's the jump condition at point a? That's the critical question. What changes at point a? You remember, this is the corresponding thing, the analog of our ramp. So what changed for the ramp at point a? What jumped? The slope. Okay, now the question is: what's going to jump here? What jumped there? Here, did the function jump? Certainly not. Did the slope jump? Certainly not. Did the second derivative jump? No, no. The second derivative was zero and then zero. What jumped? Third derivative. The third derivative is allowed to jump. And of course. A jump in the third derivative produces a delta in the fourth. Right? It just works.

So this is a cubic of some sort, coming from that jump. This is another, a different cubic, coming from this sort, from this junk, and this. So it's cubic in each piece. Why is it cubic in each piece? Because, what's the equation in the middle of that piece? What's our differential equation if I look here? Here's my equation. What is it in the middle of that piece? u fourth equal zero. The delta function is zero there. And u fourth is zero here. So of course this is a cubic spline. We're meeting that neat word, cubic spline. Those turn out to be very, very handy functions for other things, too. So we see them here as a solution to fourth order equations with point loads are cubic splines. Because the big key point is that there's a jump here in u triple prime, the third derivative. A jump in the third derivative, that's what we saw here. And we'll see it if we have any cubic meeting any cubic.

Let me just say, a jump in the third derivative, your eye probably won't notice it. I mean, it's a discontinuity, somehow. We don't have the same polynomial from here to here. But that discontinuity in u triple prime, it's pretty darn smooth still. The slope is continuous, so your eye doesn't see a ramp. And even more, the curving is continuous, the curvature is continuous. e and w are good. It's just a jump in the third derivative. Okay, so I want to speak about splines, and more about this, and about finite elements for beam problems on Friday. And then that will take care of 1-D and we'll move into 2-D. Okay.