1 00:00:00,000 --> 00:00:01,458 INTRODUCTION: The following content 2 00:00:01,458 --> 00:00:04,420 is provided by MIT OpenCourseWare under a Creative 3 00:00:04,420 --> 00:00:06,090 Commons license. 4 00:00:06,090 --> 00:00:08,230 Additional information about our license 5 00:00:08,230 --> 00:00:10,490 and MIT OpenCourseWare in general 6 00:00:10,490 --> 00:00:11,930 is available at ocw.mit.edu. 7 00:00:16,087 --> 00:00:16,670 PROFESSOR: OK. 8 00:00:16,670 --> 00:00:19,730 Should we go? 9 00:00:19,730 --> 00:00:21,550 Great. 10 00:00:21,550 --> 00:00:23,580 Well I hope you enjoyed Patriots' Day, 11 00:00:23,580 --> 00:00:28,070 and now we'll go for the final weeks of the course. 12 00:00:28,070 --> 00:00:31,500 And I wanted start with just a thought 13 00:00:31,500 --> 00:00:36,600 about scientific computing, and the people who work in it, 14 00:00:36,600 --> 00:00:39,680 and the problems. 15 00:00:39,680 --> 00:00:42,200 So just sort of an orientation. 16 00:00:42,200 --> 00:00:53,170 And maybe my point is that in part of the subject, 17 00:00:53,170 --> 00:00:56,350 the situation started like, just give me the matrix, 18 00:00:56,350 --> 00:00:57,790 tell me the equation. 19 00:00:57,790 --> 00:01:01,770 I don't want to know the physics behind it, 20 00:01:01,770 --> 00:01:05,120 just the form of the equation, and I'll 21 00:01:05,120 --> 00:01:08,000 discuss how it could be solved. 22 00:01:08,000 --> 00:01:11,870 And to contrast that with what you may often 23 00:01:11,870 --> 00:01:18,500 be meeting in your theses, even in projects, 24 00:01:18,500 --> 00:01:25,450 where the details, the physical consequences, the actual detail 25 00:01:25,450 --> 00:01:28,730 behavior of this solution is what you're looking for. 26 00:01:28,730 --> 00:01:31,070 So somehow there are those two parts of the world, 27 00:01:31,070 --> 00:01:36,510 and people who know, understand the physical problem 28 00:01:36,510 --> 00:01:38,510 are on the second side. 29 00:01:38,510 --> 00:01:44,540 And people who maybe think just about matrices, 30 00:01:44,540 --> 00:01:50,660 like that matrix for example -- that saddle point matrix, 31 00:01:50,660 --> 00:01:52,460 are more on the first side. 32 00:01:52,460 --> 00:01:55,880 Like, just give me the matrix and I'll figure out 33 00:01:55,880 --> 00:02:03,330 how to solve the linear system. 34 00:02:03,330 --> 00:02:06,170 If it's a large linear system, I'll 35 00:02:06,170 --> 00:02:10,130 invent preconditioners, study stability, 36 00:02:10,130 --> 00:02:11,350 whatever you have to do. 37 00:02:15,280 --> 00:02:18,980 But today's lecture is kind of half and half, 38 00:02:18,980 --> 00:02:23,900 because that is the framework that we're in. 39 00:02:23,900 --> 00:02:32,300 But the problem comes from a real problem in fluids. 40 00:02:32,300 --> 00:02:38,490 And it's called the Stokes problem, not Navier-Stokes. 41 00:02:38,490 --> 00:02:44,000 So Stokes alone is dealing with this question when there is 42 00:02:44,000 --> 00:02:48,930 a diffusion term -- and I've normalized the coefficient 43 00:02:48,930 --> 00:02:54,250 there, viscosity term, whatever we call it, to be 1 -- 44 00:02:54,250 --> 00:03:00,200 and there is the constraint of incompressibility, 45 00:03:00,200 --> 00:03:03,430 constant pressure, divergence p equals 0. 46 00:03:03,430 --> 00:03:07,050 I'm sorry not constant, divergence of p 47 00:03:07,050 --> 00:03:09,140 equals 0 is a constraint. 48 00:03:15,050 --> 00:03:18,810 First, let's just, so, make some comment 49 00:03:18,810 --> 00:03:22,680 on the physical side of this equation. 50 00:03:22,680 --> 00:03:27,460 What we're missing, that Navier-Stokes would put in, 51 00:03:27,460 --> 00:03:30,930 is a convection term. 52 00:03:30,930 --> 00:03:36,320 A term of v dot grad v that takes account of the fact 53 00:03:36,320 --> 00:03:42,330 that there's a flow here, things are moving. 54 00:03:42,330 --> 00:03:45,170 We're not in a stationary situation. 55 00:03:45,170 --> 00:03:48,920 And that term has been dropped. 56 00:03:48,920 --> 00:03:52,040 So that's telling us then we are in the case 57 00:03:52,040 --> 00:03:59,190 of very, very slow flow, where that term could be neglected. 58 00:03:59,190 --> 00:04:03,140 And certainly neglecting that term makes our problem easier. 59 00:04:03,140 --> 00:04:09,390 It leaves us with this symmetric matrix. 60 00:04:09,390 --> 00:04:12,500 It doesn't make it totally easy as we'll see, 61 00:04:12,500 --> 00:04:18,380 but it's a lot easier than when the convection term -- 62 00:04:18,380 --> 00:04:24,290 which by the way is not symmetric and pushes us outside 63 00:04:24,290 --> 00:04:25,370 this framework. 64 00:04:25,370 --> 00:04:28,780 If that's present, then -- of course in computational fluid 65 00:04:28,780 --> 00:04:31,830 dynamics, that's the central problem -- 66 00:04:31,830 --> 00:04:33,370 you have to deal with it. 67 00:04:33,370 --> 00:04:36,730 So we're taking a model problem that 68 00:04:36,730 --> 00:04:38,780 doesn't have the convection, doesn't have 69 00:04:38,780 --> 00:04:42,640 that unsymmetric complexity. 70 00:04:42,640 --> 00:04:47,450 But actually, it's become more and more important. 71 00:04:47,450 --> 00:04:49,570 It's more than a model problem now, 72 00:04:49,570 --> 00:04:53,250 because there are a lot of medical application, say -- 73 00:04:53,250 --> 00:05:03,230 probably the flow of your blood pushed by the heart, 74 00:05:03,230 --> 00:05:04,860 that's moving along at a decent rate. 75 00:05:04,860 --> 00:05:09,270 But flow into tiny capillaries is slow, 76 00:05:09,270 --> 00:05:12,960 and Stokes is a reasonable step. 77 00:05:12,960 --> 00:05:14,990 And then actually, Stokes is also 78 00:05:14,990 --> 00:05:18,960 going to appear in Navier-Stokes calculations, 79 00:05:18,960 --> 00:05:25,070 because often when you have two difficult processes 80 00:05:25,070 --> 00:05:27,710 you do a kind of splitting method, 81 00:05:27,710 --> 00:05:36,120 and you drop one term for half a step. 82 00:05:36,120 --> 00:05:38,390 So you take for example, a Stokes step, 83 00:05:38,390 --> 00:05:44,140 and then you do what the fluid forces you to do. 84 00:05:44,140 --> 00:05:47,121 You move that process with the flow. 85 00:05:47,121 --> 00:05:47,620 OK. 86 00:05:47,620 --> 00:05:53,050 So we're flowing slowly here. 87 00:05:53,050 --> 00:06:00,240 And the nice thing is that this fits our saddle point 88 00:06:00,240 --> 00:06:04,240 framework, that's why I'm calling this block matrix S. 89 00:06:04,240 --> 00:06:10,960 And you see that really, it came from a minimization. 90 00:06:10,960 --> 00:06:13,420 The minimization was the problem that 91 00:06:13,420 --> 00:06:18,340 leads to the Laplacian, that stands for the Laplace 92 00:06:18,340 --> 00:06:24,650 operator, with a minus sign to make it positive definite. 93 00:06:24,650 --> 00:06:28,170 And then this constraint gets imposed 94 00:06:28,170 --> 00:06:32,200 by a Lagrange multiplier p. 95 00:06:32,200 --> 00:06:34,640 So p is the Lagrange multiplier here. 96 00:06:34,640 --> 00:06:38,780 Oh, wrong. 97 00:06:38,780 --> 00:06:41,370 Divergence v is 0, incompressible. 98 00:06:48,130 --> 00:06:49,560 Now, that's on the tape. 99 00:06:49,560 --> 00:06:50,540 It's now on the tape. 100 00:06:50,540 --> 00:06:51,720 Right. 101 00:06:51,720 --> 00:06:54,500 Divergence v is 0. 102 00:06:54,500 --> 00:06:55,170 OK. 103 00:06:55,170 --> 00:07:02,630 So you remember how our equations looked? 104 00:07:02,630 --> 00:07:06,190 I'm going to write the unknowns here, 105 00:07:06,190 --> 00:07:11,330 v and p, so that you see what everything is, 106 00:07:11,330 --> 00:07:15,560 and it equals f and 0. 107 00:07:15,560 --> 00:07:18,060 OK. 108 00:07:18,060 --> 00:07:21,550 So this is now S, so I chop off here. 109 00:07:21,550 --> 00:07:25,260 So there you see the equation in the saddle point form. 110 00:07:25,260 --> 00:07:32,130 So A is the gradient, and it operates on p. 111 00:07:32,130 --> 00:07:34,160 And so p is a scalar. 112 00:07:34,160 --> 00:07:37,330 I better say -- if I put a little arrow here, 113 00:07:37,330 --> 00:07:40,130 v is a vector, it has to components say, 114 00:07:40,130 --> 00:07:47,340 if we're in a 2D problem; f has two components. 115 00:07:47,340 --> 00:07:52,020 p, the Lagrange multiplier for this constraint, the constraint 116 00:07:52,020 --> 00:07:54,820 is just a scalar constraint, the divergence is 0. 117 00:07:54,820 --> 00:08:03,100 You remember what divergence means: dv_1/dx plus dv_2/dy, 118 00:08:03,100 --> 00:08:07,150 so those are the two components of v, the vector. 119 00:08:07,150 --> 00:08:09,820 And its divergence is 0. 120 00:08:09,820 --> 00:08:14,880 And these are, well with the minus sign, 121 00:08:14,880 --> 00:08:17,750 are adjoints to each other. 122 00:08:17,750 --> 00:08:18,250 OK. 123 00:08:21,150 --> 00:08:26,340 So what's new that makes this problem worth some study? 124 00:08:29,190 --> 00:08:34,060 I guess that up to now, we've taken problems like this 125 00:08:34,060 --> 00:08:37,850 and we've done elimination to get to a single equation. 126 00:08:37,850 --> 00:08:39,940 You remember that. 127 00:08:39,940 --> 00:08:47,890 So that single equation is A transpose C*A -- let's see. 128 00:08:47,890 --> 00:08:56,510 We're multiplying that row by A transpose C and subtracting. 129 00:08:56,510 --> 00:08:58,630 Let me do the elimination under here. 130 00:08:58,630 --> 00:09:02,600 I multiply the first equation by A transpose C. You see, 131 00:09:02,600 --> 00:09:07,190 I'm doing the "just give me the matrix" mode for now. 132 00:09:07,190 --> 00:09:09,990 I multiply by A transpose C and I subtract. 133 00:09:09,990 --> 00:09:13,000 So the first equation is still there, 134 00:09:13,000 --> 00:09:18,720 C inverse v plus A*p is f. 135 00:09:18,720 --> 00:09:22,310 But the second equation, now, when I multiply by A transpose 136 00:09:22,310 --> 00:09:25,870 C, so I just put it up there, and subtract, 137 00:09:25,870 --> 00:09:29,880 that produces a 0 there, it produces minus A transpose 138 00:09:29,880 --> 00:09:37,630 C*A*p equals minus A transpose C*f. 139 00:09:40,210 --> 00:09:45,860 It leads us to this crucial matrix. 140 00:09:45,860 --> 00:09:51,140 And in applications up to now, that was the way to go. 141 00:09:51,140 --> 00:09:52,830 That was totally satisfactory. 142 00:09:52,830 --> 00:09:58,840 We eliminated v, we got an equation for p, we solved it, 143 00:09:58,840 --> 00:10:06,510 and then once we knew p we went back and found v. OK. 144 00:10:06,510 --> 00:10:09,720 So why is this different? 145 00:10:09,720 --> 00:10:13,240 Why don't I want to get to A transpose C*A, 146 00:10:13,240 --> 00:10:16,760 the stiffness matrix, if we were in finite elements. 147 00:10:16,760 --> 00:10:22,230 Why do I want to stay with this two-field system? 148 00:10:22,230 --> 00:10:23,850 Which in finite elements is going 149 00:10:23,850 --> 00:10:25,960 to be a mixed method of finite elements. 150 00:10:25,960 --> 00:10:30,720 I'm going to keep v and p in the problem for a simple reason, 151 00:10:30,720 --> 00:10:36,250 that C in this problem is awkward, 152 00:10:36,250 --> 00:10:40,010 C would be the inverse of the Laplacian. 153 00:10:40,010 --> 00:10:46,370 C is the inverse of the Laplacian d 154 00:10:46,370 --> 00:10:49,720 second by dx square plus d second by dy square. 155 00:10:52,670 --> 00:10:55,630 That's what delta stands for. 156 00:10:55,630 --> 00:11:02,800 And that's an uncomfortable operator 157 00:11:02,800 --> 00:11:04,730 to work with, the inverse here. 158 00:11:04,730 --> 00:11:09,470 In matrix language, when we discretize, 159 00:11:09,470 --> 00:11:11,830 it won't be sparse, right? 160 00:11:11,830 --> 00:11:14,570 The Laplacian is sparse, the Laplacian gives us 161 00:11:14,570 --> 00:11:18,260 the five-point matrix, the matrix with five 162 00:11:18,260 --> 00:11:23,260 non-0's on every row -- four on the diagonal and four minus 163 00:11:23,260 --> 00:11:28,650 1's, that matrix that we know, the K2D matrix. 164 00:11:28,650 --> 00:11:32,200 But this is its inverse, C is now its inverse. 165 00:11:32,200 --> 00:11:35,660 It's the inverse of that matrix K2D that we've studied, 166 00:11:35,660 --> 00:11:39,140 and of course it's full. 167 00:11:39,140 --> 00:11:47,700 So taking this step produces a completely full matrix, and all 168 00:11:47,700 --> 00:11:55,010 of our efficient algorithms for large problems 169 00:11:55,010 --> 00:11:58,460 were for sparse matrices. 170 00:11:58,460 --> 00:12:00,230 Whereas this one is sparse. 171 00:12:00,230 --> 00:12:05,320 I mean this the Laplacian. 172 00:12:05,320 --> 00:12:09,340 So we have the five-point matrix sitting here. 173 00:12:09,340 --> 00:12:12,300 Here we have a discrete gradient, 174 00:12:12,300 --> 00:12:14,730 and a discrete divergence. 175 00:12:14,730 --> 00:12:19,410 They're all derivatives, and you think of the derivatives 176 00:12:19,410 --> 00:12:22,970 as being replaced by differences. 177 00:12:22,970 --> 00:12:27,190 All these matrices are sparse, but if we 178 00:12:27,190 --> 00:12:32,510 take that fatal step of eliminating to get here, 179 00:12:32,510 --> 00:12:34,720 the matrix is full. 180 00:12:34,720 --> 00:12:39,490 Because of this fact that the inverse of a sparse matrix 181 00:12:39,490 --> 00:12:41,520 tends to be a full matrix. 182 00:12:41,520 --> 00:12:44,310 And you don't want to work with that. 183 00:12:44,310 --> 00:12:47,290 So we don't want, that's my message, 184 00:12:47,290 --> 00:12:55,390 don't want A transpose A when this is not sparse. 185 00:12:55,390 --> 00:12:55,890 OK. 186 00:12:59,220 --> 00:13:01,420 So that's the point about this example, which 187 00:13:01,420 --> 00:13:12,990 in the notes on the website is the section called 188 00:13:12,990 --> 00:13:14,890 the saddle point Stokes problem. 189 00:13:14,890 --> 00:13:15,880 OK. 190 00:13:15,880 --> 00:13:16,410 All right. 191 00:13:16,410 --> 00:13:20,590 So now we have both a specific example 192 00:13:20,590 --> 00:13:24,490 of this Stokes equations, three equations 193 00:13:24,490 --> 00:13:29,020 really, because this is two equations for v_1 and v_2, 194 00:13:29,020 --> 00:13:34,640 and this is the third equation, the constraint. 195 00:13:34,640 --> 00:13:36,130 OK. 196 00:13:36,130 --> 00:13:41,350 So we have both a specific example and the general case 197 00:13:41,350 --> 00:13:43,520 of -- 198 00:13:46,849 --> 00:13:49,140 So we're now going to work with that matrix as a whole; 199 00:13:49,140 --> 00:13:51,450 we're not going to reduce it to A transpose C, 200 00:13:51,450 --> 00:13:53,430 although in principle we could. 201 00:13:53,430 --> 00:13:58,740 And the question is one of stability. 202 00:13:58,740 --> 00:14:02,300 So like, stability is a question that 203 00:14:02,300 --> 00:14:05,990 didn't come up for our minimum principles 204 00:14:05,990 --> 00:14:07,980 when there were no constraints. 205 00:14:07,980 --> 00:14:10,790 And maybe I'll just put a little mention why. 206 00:14:10,790 --> 00:14:14,080 So this is the first time we've had to face 207 00:14:14,080 --> 00:14:17,970 the question of stability. 208 00:14:17,970 --> 00:14:22,840 What do I mean that S inverse is under control? 209 00:14:27,340 --> 00:14:33,280 I don't mean that it's easy to compute, 210 00:14:33,280 --> 00:14:35,160 the inverse of that block matrix. 211 00:14:35,160 --> 00:14:38,590 But I have to know that the matrix is invertible. 212 00:14:38,590 --> 00:14:41,890 And stability is always like one slight step 213 00:14:41,890 --> 00:14:43,710 more than just being invertible. 214 00:14:43,710 --> 00:14:46,330 It means it's invertible and you have 215 00:14:46,330 --> 00:14:49,070 some control of the inverse. 216 00:14:49,070 --> 00:14:53,960 You know, it's not just on the verge of disaster. 217 00:14:53,960 --> 00:14:58,137 So we want to know about the inverse of that matrix. 218 00:14:58,137 --> 00:15:00,470 That's essentially it, we want to know about the inverse 219 00:15:00,470 --> 00:15:01,920 of that matrix. 220 00:15:01,920 --> 00:15:04,750 And in the process, we really want 221 00:15:04,750 --> 00:15:08,880 to know about the inverse of this matrix. 222 00:15:08,880 --> 00:15:13,590 So we may not want to compute it, that was my point, 223 00:15:13,590 --> 00:15:17,690 but stability is telling us -- we have here a matrix 224 00:15:17,690 --> 00:15:20,530 that's not positive definite. 225 00:15:20,530 --> 00:15:23,590 That the key issue here. 226 00:15:23,590 --> 00:15:28,240 Our saddle point matrix is not positive definite. 227 00:15:28,240 --> 00:15:32,520 It has positive eigenvalues and negative eigenvalues. 228 00:15:32,520 --> 00:15:35,060 Let me just take some examples. 229 00:15:35,060 --> 00:15:35,920 OK. 230 00:15:35,920 --> 00:15:42,970 So here's an S. So I have a positive thing there. 231 00:15:42,970 --> 00:15:46,080 Nonzero stuff there. 232 00:15:46,080 --> 00:15:50,790 There is the simplest example I could think of. 233 00:15:50,790 --> 00:15:52,720 That matrix is not positive definite. 234 00:15:55,640 --> 00:15:59,310 It's determinant is minus 1, actually. 235 00:15:59,310 --> 00:16:06,720 But its eigenvalues are, I mean, they're safely away from 0. 236 00:16:06,720 --> 00:16:09,220 it's not singular or close to singular either. 237 00:16:13,750 --> 00:16:17,390 So that two by two example, no problem. 238 00:16:17,390 --> 00:16:28,050 This equation -- well, I could study the invertibility of this 239 00:16:28,050 --> 00:16:32,640 particular matrix, that would be solving differential equations. 240 00:16:32,640 --> 00:16:35,360 In other words, I could study the case -- 241 00:16:35,360 --> 00:16:41,480 are p and v under control if I have a bound on f? 242 00:16:41,480 --> 00:16:43,490 Well, they are. 243 00:16:43,490 --> 00:16:48,190 But the issue comes when I go to finite elements, 244 00:16:48,190 --> 00:16:50,890 when I look at a subspace. 245 00:16:50,890 --> 00:16:53,300 So have to try to make that clear. 246 00:16:57,570 --> 00:17:01,190 What I plan to do to solve a continuous problem 247 00:17:01,190 --> 00:17:04,860 is to discretize it. 248 00:17:04,860 --> 00:17:12,230 This is my continuous matrix, block matrix, symbolic 249 00:17:12,230 --> 00:17:14,750 of the differential equations. 250 00:17:14,750 --> 00:17:19,600 And it's OK, analysis would show that. 251 00:17:19,600 --> 00:17:23,720 But the question is: if I approximate v 252 00:17:23,720 --> 00:17:33,830 by some combination of finite elements, let me call it v_h. 253 00:17:33,830 --> 00:17:39,480 If I approximate the pressure by some pressure element p_h, 254 00:17:39,480 --> 00:17:42,660 then I'm sort of projecting the problem down 255 00:17:42,660 --> 00:17:47,080 into these finite-dimensional, finite element subspaces. 256 00:17:47,080 --> 00:17:50,310 And of course f will get projected down to f_h. 257 00:17:50,310 --> 00:17:53,840 And the matrices, these operators, 258 00:17:53,840 --> 00:17:59,840 will turn into finite sparse matrices. 259 00:17:59,840 --> 00:18:03,470 And the question is that projection safe? 260 00:18:03,470 --> 00:18:05,150 So that's really it. 261 00:18:05,150 --> 00:18:07,620 It's the stability. 262 00:18:07,620 --> 00:18:10,890 S inverse is under control, this is OK. 263 00:18:10,890 --> 00:18:16,810 But is the projection of S inverse under control? 264 00:18:16,810 --> 00:18:24,270 Is the projection of S -- so if I project it, 265 00:18:24,270 --> 00:18:28,240 is the inverse of that under control? 266 00:18:28,240 --> 00:18:29,100 That's the question. 267 00:18:29,100 --> 00:18:38,810 And this projection is -- what comes up in applications, 268 00:18:38,810 --> 00:18:45,930 this is a projection of the continuous problem, 269 00:18:45,930 --> 00:18:54,940 the differential equation, onto finite element spaces, 270 00:18:54,940 --> 00:18:56,550 finite element trial function. 271 00:19:00,200 --> 00:19:00,700 OK. 272 00:19:00,700 --> 00:19:04,640 So we've spoken a little about finite elements in this course, 273 00:19:04,640 --> 00:19:08,749 and more in 18.085, and you're probably 274 00:19:08,749 --> 00:19:09,790 have seen them elsewhere. 275 00:19:09,790 --> 00:19:18,060 You know that we can approximate velocities by -- well, 276 00:19:18,060 --> 00:19:23,570 here would be a really important example, 277 00:19:23,570 --> 00:19:27,590 because it's an especially simple finite elements. 278 00:19:27,590 --> 00:19:31,780 If I take my region, suppose it happened to be that rectangle, 279 00:19:31,780 --> 00:19:36,590 I divide that rectangle into small squares, 280 00:19:36,590 --> 00:19:39,740 there's a small square. 281 00:19:39,740 --> 00:19:43,150 I'm using squares rather than triangles. 282 00:19:43,150 --> 00:19:50,090 And in each square I make some low-degree approximation, 283 00:19:50,090 --> 00:19:54,320 low-degree trial function, for the velocities 284 00:19:54,320 --> 00:19:55,240 and the pressures. 285 00:19:55,240 --> 00:20:01,100 And the question is what functions can I use? 286 00:20:01,100 --> 00:20:05,420 For example, one I would love to be able to use would be, 287 00:20:05,420 --> 00:20:09,820 so in each square I would like to be 288 00:20:09,820 --> 00:20:14,480 able to use piecewise constant pressures. 289 00:20:14,480 --> 00:20:18,920 Piecewise constant pressure, so p_h equal constant 290 00:20:18,920 --> 00:20:23,920 in each element. 291 00:20:23,920 --> 00:20:25,860 Element is the square. 292 00:20:28,420 --> 00:20:31,900 So a different constant in each element. 293 00:20:31,900 --> 00:20:36,900 So that would give me my vector of p_h's. 294 00:20:36,900 --> 00:20:42,150 There would be one component for every little square. 295 00:20:42,150 --> 00:20:45,390 And now what about v_h? 296 00:20:45,390 --> 00:20:50,850 What trial function would I like to use for the velocities? 297 00:20:50,850 --> 00:20:56,120 Constant velocities wouldn't be good, 298 00:20:56,120 --> 00:20:58,660 that's further than I'm willing to go. 299 00:20:58,660 --> 00:21:04,420 I'll go piecewise linear in each square. 300 00:21:04,420 --> 00:21:08,750 Now piecewise linear in each triangle meant that it had 301 00:21:08,750 --> 00:21:11,120 the form a plus b*x plus c*y. 302 00:21:13,790 --> 00:21:16,010 Strictly speaking that's what I should 303 00:21:16,010 --> 00:21:19,480 mean by piecewise linear, with a different coefficient a, b, 304 00:21:19,480 --> 00:21:21,380 c in each piece. 305 00:21:21,380 --> 00:21:24,590 But because we're in squares now, 306 00:21:24,590 --> 00:21:31,410 we've got four nodal values, and we want a fourth coefficient. 307 00:21:31,410 --> 00:21:33,820 This was appropriate on triangles, 308 00:21:33,820 --> 00:21:39,360 on a square we really include the bilinear term d*x*y. 309 00:21:39,360 --> 00:21:40,870 OK. 310 00:21:40,870 --> 00:21:49,080 So we have four coefficients which 311 00:21:49,080 --> 00:21:51,730 determine and can be determined from the four 312 00:21:51,730 --> 00:21:53,340 values of the node. 313 00:21:53,340 --> 00:21:57,700 You see that I'm describing the whole space. 314 00:21:57,700 --> 00:22:02,970 Now I cut the region into more squares, 315 00:22:02,970 --> 00:22:12,040 and I've described a way of discretizing the differential 316 00:22:12,040 --> 00:22:12,910 difference equation. 317 00:22:12,910 --> 00:22:16,470 The weak form of continuous problem, 318 00:22:16,470 --> 00:22:22,240 when I plug in these functions, would give me 319 00:22:22,240 --> 00:22:24,320 a matrix problem, a discrete problem, 320 00:22:24,320 --> 00:22:25,510 that will have this form. 321 00:22:25,510 --> 00:22:32,100 And the question is do I have any control? 322 00:22:32,100 --> 00:22:36,020 So it will give me this projected matrix now, 323 00:22:36,020 --> 00:22:39,630 finite-dimensional instead of a continuous problem. 324 00:22:39,630 --> 00:22:42,460 And the question is: have I still got any control? 325 00:22:44,990 --> 00:22:50,870 Well the answer for this particular choice of pressure 326 00:22:50,870 --> 00:22:55,950 space and velocity space is no. 327 00:22:55,950 --> 00:23:04,560 Unfortunately the matrix S that appears instead, 328 00:23:04,560 --> 00:23:12,110 the C inverse h, shall I say, it would 329 00:23:12,110 --> 00:23:16,560 be like our K2D inverse, that was pretty 330 00:23:16,560 --> 00:23:18,300 close of what it would be. 331 00:23:18,300 --> 00:23:25,680 And the A_h and the A transpose h and the 0. 332 00:23:25,680 --> 00:23:33,600 That matrix has a famous null vector 333 00:23:33,600 --> 00:23:35,380 for this specific example. 334 00:23:35,380 --> 00:23:45,720 I'm taking big jumps here, because I just 335 00:23:45,720 --> 00:23:48,330 want to say the main point. 336 00:23:48,330 --> 00:23:52,370 The main point here is that a null vector 337 00:23:52,370 --> 00:23:54,280 appears in this problem. 338 00:23:54,280 --> 00:23:59,120 S fails to be invertible, and that null vector is like -- 339 00:23:59,120 --> 00:24:04,380 finite element people know its name. 340 00:24:04,380 --> 00:24:08,470 It has a couple of names actually, maybe checkerboard. 341 00:24:08,470 --> 00:24:15,710 The null vector is a vector when the pressure is 1, minus 1, 1, 342 00:24:15,710 --> 00:24:22,740 minus 1 in a checkerboard pattern over the whole region, 343 00:24:22,740 --> 00:24:27,830 and a corresponding velocity. 344 00:24:27,830 --> 00:24:36,400 Well, the trouble is A_h times this p is 0. 345 00:24:36,400 --> 00:24:40,190 That's what kills you. 346 00:24:40,190 --> 00:24:47,180 It turns out that A times that particular pressure is 0. 347 00:24:47,180 --> 00:24:50,230 And if you lose control of A that way, 348 00:24:50,230 --> 00:24:52,660 you've lost control A transpose CA. 349 00:24:52,660 --> 00:24:55,240 You've lost everything. 350 00:24:55,240 --> 00:25:04,420 So I guess, coming back to my opening comments. 351 00:25:04,420 --> 00:25:07,700 Here's a case in which we have a physical problem, 352 00:25:07,700 --> 00:25:18,470 but it falls into a pattern, where the mathematics actually 353 00:25:18,470 --> 00:25:21,890 decides, are we able to use the trial spaces that 354 00:25:21,890 --> 00:25:25,210 are most convenient? 355 00:25:25,210 --> 00:25:28,220 I say convenient rather than most accurate of course. 356 00:25:28,220 --> 00:25:30,940 If I wanted to improve the accuracy, 357 00:25:30,940 --> 00:25:33,970 I would go up in the degree, I would 358 00:25:33,970 --> 00:25:37,640 take pressures to be piecewise linear and velocities 359 00:25:37,640 --> 00:25:43,510 to be piecewise quadratic with more points. 360 00:25:43,510 --> 00:25:50,150 And I could check the stability or not, 361 00:25:50,150 --> 00:25:53,870 is the matrix invertible or not, and the answer 362 00:25:53,870 --> 00:25:55,840 might be better there. 363 00:25:55,840 --> 00:26:05,010 In fact the notes and any book on finite elements has a list. 364 00:26:05,010 --> 00:26:08,880 And I don't plan on making that list here. 365 00:26:08,880 --> 00:26:15,550 But it has a list of pressure spaces, pressure elements 366 00:26:15,550 --> 00:26:23,910 and velocity elements, on which the first guy that I've 367 00:26:23,910 --> 00:26:28,270 mentioned would be like Q of 0, degree 0. 368 00:26:28,270 --> 00:26:32,820 and velocities Q_1 degree 1, and the letter 369 00:26:32,820 --> 00:26:36,950 Q being used for quadrilateral. 370 00:26:36,950 --> 00:26:41,370 If I divided it into triangles, I would use the letter P 371 00:26:41,370 --> 00:26:44,300 for polynomial in that case. 372 00:26:44,300 --> 00:26:46,680 And this is a fail. 373 00:26:46,680 --> 00:26:50,840 This one fails. 374 00:26:50,840 --> 00:26:55,910 And then another list would be how about Q_1, Q_1? 375 00:26:55,910 --> 00:26:56,960 I don't know. 376 00:26:56,960 --> 00:27:00,720 I wouldn't bet much on Q_1, Q_1, and of course I 377 00:27:00,720 --> 00:27:03,830 can look at the notes to see what the answer. 378 00:27:03,830 --> 00:27:07,080 Or Q_1, Q_2 that's probably got a better chance. 379 00:27:07,080 --> 00:27:14,060 Or many, many -- you know, people are inventing finite 380 00:27:14,060 --> 00:27:17,790 elements all the time, and they're having to check 381 00:27:17,790 --> 00:27:23,760 the stability or not, the uniform invertibility or not 382 00:27:23,760 --> 00:27:26,810 of S. OK. 383 00:27:26,810 --> 00:27:31,000 I just want to say that this one is so attractive, because it's 384 00:27:31,000 --> 00:27:35,190 so simple that people don't give up on it. 385 00:27:35,190 --> 00:27:43,280 They want to go with this, this is a fail, but -- 386 00:27:43,280 --> 00:27:52,100 so stabilize it. 387 00:27:52,100 --> 00:27:54,900 That one is so attractive. 388 00:27:54,900 --> 00:28:03,910 And you know specifically that it's this checkerboard mode 389 00:28:03,910 --> 00:28:06,520 that's in the null space, so you just 390 00:28:06,520 --> 00:28:11,150 think, how am I going to get that mode out of the null? 391 00:28:11,150 --> 00:28:18,310 And the stabilizing method is you stick a matrix here 392 00:28:18,310 --> 00:28:21,490 for purely numerical reasons, and it'll 393 00:28:21,490 --> 00:28:27,010 be a matrix that exactly has this checkerboard pattern, so 394 00:28:27,010 --> 00:28:31,350 that the checkerboard is not any longer in the null space. 395 00:28:31,350 --> 00:28:39,270 Maybe I call it minus some small number times some E matrix, 396 00:28:39,270 --> 00:28:47,660 it's often written E. And alpha you keep as small as you dare. 397 00:28:47,660 --> 00:28:55,270 But the effect is that it takes this checkerboard pressure out 398 00:28:55,270 --> 00:28:57,900 of the null space. 399 00:29:01,620 --> 00:29:06,650 You see the problem was that without that, the columns of A 400 00:29:06,650 --> 00:29:09,100 were not independent. 401 00:29:09,100 --> 00:29:11,310 A times P could be 0. 402 00:29:11,310 --> 00:29:15,340 And if those columns are not independent, you're lost. 403 00:29:15,340 --> 00:29:18,190 Let me put the 0 back. 404 00:29:18,190 --> 00:29:20,890 Why was that a bad thing? 405 00:29:20,890 --> 00:29:25,620 I want to now sort of connect it to the matrix theory. 406 00:29:25,620 --> 00:29:30,350 If I have this particular pressure and A_h*p is 0, 407 00:29:30,350 --> 00:29:34,820 so that tells me that those columns of A_h are dependent. 408 00:29:37,420 --> 00:29:40,490 Then the columns of S are dependent, right. 409 00:29:40,490 --> 00:29:45,160 If I have something in the null space of A_h, 410 00:29:45,160 --> 00:29:48,610 then I just put -- a pressure, this pressure, 411 00:29:48,610 --> 00:29:52,130 in the null space of A_h -- then I just take the v to be there 412 00:29:52,130 --> 00:29:54,500 0, then these guys don't do anything. 413 00:29:54,500 --> 00:30:00,280 And that A times that bad pressure gives 0, so I have 414 00:30:00,280 --> 00:30:04,870 (v, p), a velocity, pressure, that's in the null space. 415 00:30:04,870 --> 00:30:05,780 OK. 416 00:30:05,780 --> 00:30:12,580 So the way it gets fixed is to stick in some small stabilizing 417 00:30:12,580 --> 00:30:14,080 factor there. 418 00:30:14,080 --> 00:30:15,550 OK. 419 00:30:15,550 --> 00:30:16,180 All right. 420 00:30:16,180 --> 00:30:24,170 Now, that's the specifics of this particular example. 421 00:30:24,170 --> 00:30:27,010 Well, I didn't give all the specifics. 422 00:30:27,010 --> 00:30:32,440 If I were really on the ball I would show in detail 423 00:30:32,440 --> 00:30:39,540 why that checkerboard appears and kills this choice. 424 00:30:39,540 --> 00:30:42,530 I would check each of these other choices, 425 00:30:42,530 --> 00:30:49,990 and beside each one I would put fail or success. 426 00:30:49,990 --> 00:30:57,070 And when it fails, I would try to identify the failure mode, 427 00:30:57,070 --> 00:31:00,940 so that if you wanted you could stabilize and get 428 00:31:00,940 --> 00:31:02,950 rid of that failure mode. 429 00:31:02,950 --> 00:31:03,450 OK. 430 00:31:09,380 --> 00:31:15,470 What I'll do instead is to come to the general question here. 431 00:31:15,470 --> 00:31:19,340 So the general question is to deal with this matrix, 432 00:31:19,340 --> 00:31:20,990 S. Can I write it again? 433 00:31:20,990 --> 00:31:26,020 C inverse, A; A transpose, 0. 434 00:31:26,020 --> 00:31:30,260 I want to look it its -- is it invertible? 435 00:31:30,260 --> 00:31:32,330 Can I get some bound on its inverse? 436 00:31:32,330 --> 00:31:33,790 What do I have to know? 437 00:31:33,790 --> 00:31:38,740 And I'll just emphasize here, what's the real problem? 438 00:31:38,740 --> 00:31:40,440 Because we saw it right there. 439 00:31:40,440 --> 00:31:49,600 The real problem is if the columns of A are dependent. 440 00:31:54,350 --> 00:31:58,480 In other words, i.e., in other words, 441 00:31:58,480 --> 00:32:05,760 A*p equals 0 has a nonzero solution. 442 00:32:05,760 --> 00:32:08,140 Then you're wiped out. 443 00:32:10,740 --> 00:32:11,620 I'll just repeat. 444 00:32:11,620 --> 00:32:15,570 If the columns of A are dependent, 445 00:32:15,570 --> 00:32:18,400 then the columns of S are dependent, 446 00:32:18,400 --> 00:32:21,620 the matrix is singular, and you're finished. 447 00:32:21,620 --> 00:32:26,170 So somehow, we have to get some muscle behind A. 448 00:32:26,170 --> 00:32:29,780 But remember, A is a rectangular matrix. 449 00:32:29,780 --> 00:32:35,290 This C is m by m, this A that comes up 450 00:32:35,290 --> 00:32:40,200 in the constraint is m by n. 451 00:32:40,200 --> 00:32:47,340 It has, we hope, more rows than columns. 452 00:32:51,990 --> 00:32:55,630 But still, the columns could be dependent, 453 00:32:55,630 --> 00:32:57,530 and in this situation they are. 454 00:32:57,530 --> 00:32:58,320 OK. 455 00:32:58,320 --> 00:33:01,910 So the question is, we have to get some muscle behind A. 456 00:33:01,910 --> 00:33:07,630 And this is the famous condition so called inf-sup condition 457 00:33:07,630 --> 00:33:11,850 that you would find in Professor Bathe's finite element book. 458 00:33:11,850 --> 00:33:18,850 But maybe these notes kind of focus on the heart of the idea. 459 00:33:18,850 --> 00:33:27,270 And again, in some way it's conditioned 460 00:33:27,270 --> 00:33:34,700 on A really, and on the pressure spaces and the velocity spaces 461 00:33:34,700 --> 00:33:36,670 that you choose. 462 00:33:36,670 --> 00:33:37,370 OK. 463 00:33:37,370 --> 00:33:42,570 So these, you could really say are A_h, p_h, and v_h, 464 00:33:42,570 --> 00:33:46,010 but I'll skip writing h's on everything now. 465 00:33:46,010 --> 00:33:49,140 So I'm down in the discrete case, 466 00:33:49,140 --> 00:33:56,040 but I have a sequence of discrete problem 467 00:33:56,040 --> 00:33:58,490 which are growing as h gets smaller, 468 00:33:58,490 --> 00:34:05,100 and the question is do I have a fixed beta that control this? 469 00:34:05,100 --> 00:34:11,020 And it turns out that this is the requirement that -- 470 00:34:11,020 --> 00:34:15,380 this inf-sup condition says that for every pressure, 471 00:34:15,380 --> 00:34:23,050 there is a velocity -- every discrete pressure now, 472 00:34:23,050 --> 00:34:27,160 there is a discrete velocity such that this holds. 473 00:34:27,160 --> 00:34:31,840 So it's sort of forcing A to stay away from 0. 474 00:34:31,840 --> 00:34:35,820 This beta is a fixed number, and this sort of 475 00:34:35,820 --> 00:34:41,640 is the right normalization that goes with it. 476 00:34:41,640 --> 00:34:44,640 So you see why it's v transpose A*p? 477 00:34:44,640 --> 00:34:52,570 Because A -- you know, if I look at the quadratic form, 478 00:34:52,570 --> 00:34:56,520 v transpose p transpose -- you know, 479 00:34:56,520 --> 00:35:00,280 the quadratic form with the matrix is always multiply 480 00:35:00,280 --> 00:35:03,570 by a vector here, multiply by its transpose here, 481 00:35:03,570 --> 00:35:06,800 where is the A term? 482 00:35:06,800 --> 00:35:10,980 A will multiply p, and it'll be multiplied by v transpose. 483 00:35:10,980 --> 00:35:15,270 In the quadratic part, in the energy, 484 00:35:15,270 --> 00:35:19,290 is this v transpose A*p. 485 00:35:19,290 --> 00:35:23,150 And that's what we can't let go to 0. 486 00:35:23,150 --> 00:35:29,800 Because we are not -- if that part fails us, 487 00:35:29,800 --> 00:35:33,070 if we lose these and we have a 0 there, 488 00:35:33,070 --> 00:35:38,850 this part on the submatrix can't save us. 489 00:35:38,850 --> 00:35:40,530 Right. 490 00:35:40,530 --> 00:35:42,500 So before I say anything more about that, 491 00:35:42,500 --> 00:35:45,740 let me just emphasize what stability 492 00:35:45,740 --> 00:35:52,880 was in the symmetric case. 493 00:35:52,880 --> 00:35:58,290 So suppose I have a positive definite symmetric matrix. 494 00:36:01,630 --> 00:36:10,590 It's got some eigenvalues, lambda_1 to lambda_n. 495 00:36:10,590 --> 00:36:12,800 OK. 496 00:36:12,800 --> 00:36:14,825 What you see on this board is going 497 00:36:14,825 --> 00:36:19,920 to explain why stability is automatic 498 00:36:19,920 --> 00:36:23,390 if I'm dealing with positive definite things. 499 00:36:23,390 --> 00:36:25,320 So suppose I project. 500 00:36:25,320 --> 00:36:33,540 I take this matrix -- let me just project in the most direct 501 00:36:33,540 --> 00:36:34,060 way. 502 00:36:34,060 --> 00:36:38,600 I'll throw away the last row and column. 503 00:36:38,600 --> 00:36:41,250 Suppose I throw away the last row and column. 504 00:36:41,250 --> 00:36:50,340 That part is now gone, so these are the eigenvalues of K. 505 00:36:50,340 --> 00:36:52,920 And of course, being positive definite, 506 00:36:52,920 --> 00:36:54,760 they're bigger than 0. 507 00:36:54,760 --> 00:36:57,030 OK. 508 00:36:57,030 --> 00:36:59,860 So that was the whole matrix, but now what I'm going to do 509 00:36:59,860 --> 00:37:03,730 is I'm going to project on the first n minus 1 510 00:37:03,730 --> 00:37:05,150 rows and columns. 511 00:37:05,150 --> 00:37:10,680 So I'm going to look at the submatrix K sub n minus 1. 512 00:37:10,680 --> 00:37:11,870 And it has some eigenvalues. 513 00:37:14,650 --> 00:37:16,720 Are they the same as those? 514 00:37:16,720 --> 00:37:18,370 No way. 515 00:37:18,370 --> 00:37:24,290 The eigenvalues of this submatrix will be -- 516 00:37:24,290 --> 00:37:26,100 and you see, I'm not in this situation, 517 00:37:26,100 --> 00:37:28,660 this was not positive definite. 518 00:37:28,660 --> 00:37:30,910 So I'm looking here at the positive definite case. 519 00:37:30,910 --> 00:37:38,230 This has some eigenvalues, say mu_1 up to mu_(n-1). 520 00:37:38,230 --> 00:37:45,800 These will be the eigenvalues of the submatrix, K_(n-1). 521 00:37:45,800 --> 00:37:50,510 And linear algebra is just like the right place 522 00:37:50,510 --> 00:37:56,630 to look for some information on the connection 523 00:37:56,630 --> 00:37:59,650 between the mus and lambdas. 524 00:37:59,650 --> 00:38:00,850 They're certainly not equal. 525 00:38:03,550 --> 00:38:07,070 Our first question is if I take a positive definite matrix 526 00:38:07,070 --> 00:38:09,720 with positive eigenvalues, I throw away 527 00:38:09,720 --> 00:38:14,370 the last row and column. 528 00:38:14,370 --> 00:38:16,590 Do I still have a positive definite? 529 00:38:16,590 --> 00:38:20,100 Is that submatrix positive definite? 530 00:38:20,100 --> 00:38:22,110 And the answer is yes. 531 00:38:22,110 --> 00:38:26,240 The eigenvalues are all positive. 532 00:38:26,240 --> 00:38:31,880 And more than that, how about that smallest eigenvalue? 533 00:38:31,880 --> 00:38:35,490 Because the eigenvalue is measuring like how close 534 00:38:35,490 --> 00:38:39,080 are we to disaster, to singular. 535 00:38:39,080 --> 00:38:42,740 Can I say anything about that smallest eigenvalue? 536 00:38:42,740 --> 00:38:45,210 And the answer is yes. 537 00:38:45,210 --> 00:38:52,920 The smallest eigenvalue of K_(n-1) is never below 538 00:38:52,920 --> 00:38:57,270 the smallest eigenvalue of the full matrix. 539 00:38:57,270 --> 00:39:01,030 So if the full matrix is OK, all the projections, 540 00:39:01,030 --> 00:39:07,910 all the subspace, all the submatrices are even more OK. 541 00:39:07,910 --> 00:39:16,650 And up at the other end, at the top end, 542 00:39:16,650 --> 00:39:19,870 the top eigenvalue in the submatrix 543 00:39:19,870 --> 00:39:23,740 never passes the top eigenvalue in the big matrix. 544 00:39:23,740 --> 00:39:29,070 So you see, at both ends I have what I want, I have control. 545 00:39:29,070 --> 00:39:33,710 If I have some range of eigenvalues for the big matrix, 546 00:39:33,710 --> 00:39:39,170 then all its projections are inside that range. 547 00:39:39,170 --> 00:39:40,790 And therefore the condition number, 548 00:39:40,790 --> 00:39:44,360 which would be the ratio of this to this, 549 00:39:44,360 --> 00:39:49,860 is smaller than the condition number of the whole matrix. 550 00:39:49,860 --> 00:39:57,515 So the condition number of K_(n-1) is less or equal 551 00:39:57,515 --> 00:39:59,810 the condition number of K_n. 552 00:39:59,810 --> 00:40:03,000 OK. 553 00:40:03,000 --> 00:40:08,180 And that's exactly why working with positive definite matrices 554 00:40:08,180 --> 00:40:11,730 is so comfortable. 555 00:40:11,730 --> 00:40:15,850 The stability is just automatic. 556 00:40:15,850 --> 00:40:19,190 But here we're not working with positive definite matrices. 557 00:40:19,190 --> 00:40:21,800 The stability is not automatic, and we 558 00:40:21,800 --> 00:40:25,160 need this inf-sup condition. 559 00:40:25,160 --> 00:40:27,170 And I don't know whether to try. 560 00:40:27,170 --> 00:40:30,700 Maybe I'm -- I'm a little hesitant to go through 561 00:40:30,700 --> 00:40:34,620 the steps that lead from there to there, 562 00:40:34,620 --> 00:40:41,170 because this lecture is already asking you to put a lot 563 00:40:41,170 --> 00:40:43,510 of things together without details. 564 00:40:43,510 --> 00:40:46,310 And the notes do the best job I could 565 00:40:46,310 --> 00:40:52,350 do in making those steps simple and clear. 566 00:40:52,350 --> 00:40:55,480 And this is what you actually check. 567 00:40:55,480 --> 00:41:02,860 You check this condition for any proposed finite element, 568 00:41:02,860 --> 00:41:06,390 and any different problem. 569 00:41:06,390 --> 00:41:08,440 It wouldn't have to be the Stokes problem. 570 00:41:08,440 --> 00:41:13,230 Other problems would give us a constraint matrix 571 00:41:13,230 --> 00:41:22,300 A, a sort of stiffness submatrix C inverse, 572 00:41:22,300 --> 00:41:28,631 and we could check the condition again. 573 00:41:28,631 --> 00:41:29,130 OK. 574 00:41:29,130 --> 00:41:40,150 So maybe I'm going to leave that statement for the notes. 575 00:41:40,150 --> 00:41:41,450 OK. 576 00:41:41,450 --> 00:41:45,570 So that's this inf-sup condition. 577 00:41:45,570 --> 00:41:53,940 And of course, you see right away that if A had dependent 578 00:41:53,940 --> 00:41:59,260 columns, which is the big risk, if A had dependent columns, 579 00:41:59,260 --> 00:42:04,490 that would mean that A*p is 0 for some p. 580 00:42:04,490 --> 00:42:07,970 And then there would be no v. No v could possibly 581 00:42:07,970 --> 00:42:10,310 be found that would save it. 582 00:42:10,310 --> 00:42:17,260 If A had dependent columns, that's bogeyman 583 00:42:17,260 --> 00:42:20,610 that we're avoiding above all. 584 00:42:20,610 --> 00:42:23,380 If A had dependent columns, then there would be a solution 585 00:42:23,380 --> 00:42:27,730 to A*p equals 0, this left side would be 0, 586 00:42:27,730 --> 00:42:36,530 the right side wouldn't, and the inf-sup condition would fail. 587 00:42:36,530 --> 00:42:38,940 But the inf-sup condition tells us, well, 588 00:42:38,940 --> 00:42:44,060 what about if there are pressures, that -- 589 00:42:44,060 --> 00:42:50,350 where A*p get smaller and smaller as the mesh size 590 00:42:50,350 --> 00:42:54,410 changes, then the question is can I find a v and so on. 591 00:42:54,410 --> 00:42:58,820 So that's the full condition, but not 592 00:42:58,820 --> 00:43:05,590 to lose sight of the fact that the central question, 593 00:43:05,590 --> 00:43:10,720 the key question is the independence or dependence 594 00:43:10,720 --> 00:43:12,850 of the columns of A. OK. 595 00:43:12,850 --> 00:43:13,960 Good. 596 00:43:13,960 --> 00:43:20,720 So let's see, maybe that's as much as I 597 00:43:20,720 --> 00:43:23,300 want to say about the Stokes example, 598 00:43:23,300 --> 00:43:29,470 and what it leads to in the finite element world. 599 00:43:29,470 --> 00:43:35,520 I guess I'm going to stay for another lecture or two 600 00:43:35,520 --> 00:43:41,180 with these minimizations. 601 00:43:41,180 --> 00:43:44,190 So maybe I'll just anticipate slightly 602 00:43:44,190 --> 00:43:46,740 the start of the next lecture. 603 00:43:50,140 --> 00:43:56,930 Because I'm quite pleased about getting this next topic 604 00:43:56,930 --> 00:43:59,210 into the series of lectures. 605 00:44:03,810 --> 00:44:04,770 Let me say what it is. 606 00:44:07,480 --> 00:44:16,810 So it's going to be a minimization of an A*u minus b 607 00:44:16,810 --> 00:44:19,520 squared. 608 00:44:19,520 --> 00:44:20,170 OK. 609 00:44:27,750 --> 00:44:33,860 Just as this lecture has seen, the possibility 610 00:44:33,860 --> 00:44:39,230 comes up that the columns of A are not independent, 611 00:44:39,230 --> 00:44:41,200 that this problem is in trouble. 612 00:44:41,200 --> 00:44:45,000 So I want to regularize this problem. 613 00:44:45,000 --> 00:44:48,460 How do I save a problem in which -- 614 00:44:48,460 --> 00:45:01,230 so this is a major application -- in which A is, say -- 615 00:45:01,230 --> 00:45:05,620 maybe A transpose A is the easiest thing to say -- 616 00:45:05,620 --> 00:45:11,810 A transpose A is nearly singular, 617 00:45:11,810 --> 00:45:14,270 or just plain singular. 618 00:45:14,270 --> 00:45:17,550 I'll say or singular. 619 00:45:20,680 --> 00:45:24,080 How do we rescue the problem? 620 00:45:24,080 --> 00:45:27,090 Because numerically and theoretically 621 00:45:27,090 --> 00:45:30,220 it's a disaster as it stands. 622 00:45:30,220 --> 00:45:31,650 OK. 623 00:45:31,650 --> 00:45:32,660 It's ill posed. 624 00:45:35,420 --> 00:45:40,360 So let me introduce that word, ill-posed problem. 625 00:45:40,360 --> 00:45:43,670 And let me mention, since I'm just 626 00:45:43,670 --> 00:45:48,320 using these last few minutes as a kind of getting 627 00:45:48,320 --> 00:45:51,470 the big picture, where do ill-posed problems come from, 628 00:45:51,470 --> 00:45:53,380 where does these things come from? 629 00:45:53,380 --> 00:45:58,480 They come typically from what are called inverse problems. 630 00:46:02,270 --> 00:46:08,697 So the application is to inverse problems. 631 00:46:08,697 --> 00:46:09,780 What's an inverse problem? 632 00:46:14,220 --> 00:46:18,100 An inverse problem is, for a differential equation, 633 00:46:18,100 --> 00:46:22,830 you know some inputs and outputs, 634 00:46:22,830 --> 00:46:26,190 you know some solutions, but you don't know the equation. 635 00:46:26,190 --> 00:46:28,520 You don't know the coefficient. 636 00:46:34,720 --> 00:46:37,280 So you're trying to determine what 637 00:46:37,280 --> 00:46:41,410 are the material properties by looking 638 00:46:41,410 --> 00:46:49,550 to see what the outputs are for experimental inputs. 639 00:46:52,280 --> 00:46:56,650 And very frequently, you don't have enough inputs to really 640 00:46:56,650 --> 00:47:01,170 determine what these -- so it's inverse problem to find 641 00:47:01,170 --> 00:47:16,090 the coefficients of the equation given some information on its 642 00:47:16,090 --> 00:47:19,330 solutions, some measurements. 643 00:47:19,330 --> 00:47:23,740 And of course they're probably noisy measurements. 644 00:47:23,740 --> 00:47:31,260 I mean this is an enormous source of applications, growing 645 00:47:31,260 --> 00:47:37,030 in importance, and how do we deal with it. 646 00:47:37,030 --> 00:47:40,890 Well, you regularize it. 647 00:47:40,890 --> 00:47:48,880 So this is typical of this type of problem. 648 00:47:48,880 --> 00:47:53,790 You add on some multiples of something that's safe, 649 00:47:53,790 --> 00:47:56,110 like for example u square. 650 00:48:01,630 --> 00:48:12,720 So the result is, here, that this A transpose A that 651 00:48:12,720 --> 00:48:18,610 comes from this problem, there's an alpha times the identity 652 00:48:18,610 --> 00:48:19,850 coming from this problem. 653 00:48:19,850 --> 00:48:22,760 I mean it's an ordinary least squares problem, 654 00:48:22,760 --> 00:48:29,150 and it leads us to the matrix A transpose A, 655 00:48:29,150 --> 00:48:33,060 comes from the first part. 656 00:48:33,060 --> 00:48:36,890 But I'm just talking about the standard normal equations. 657 00:48:36,890 --> 00:48:39,410 It will produce an A transpose A from here. 658 00:48:39,410 --> 00:48:43,570 But from here it will produce an alpha times the identity. 659 00:48:43,570 --> 00:48:51,340 And of course, just the way this minus E worked up 660 00:48:51,340 --> 00:48:55,280 there, here I'm doing it even more directly, 661 00:48:55,280 --> 00:49:01,350 I'm pushing the problem away from 0. 662 00:49:01,350 --> 00:49:05,320 This alpha, this parameter and the choice of that parameter, 663 00:49:05,320 --> 00:49:11,370 is the million dollar question in this whole subject. 664 00:49:11,370 --> 00:49:14,520 But its purpose is to regularize the problem, 665 00:49:14,520 --> 00:49:19,350 to take this nearly singular matrix and push it up. 666 00:49:19,350 --> 00:49:22,370 You have to see, this matrix, A transpose A, 667 00:49:22,370 --> 00:49:27,280 never has negative eigenvalues, the problem is 0 or very small, 668 00:49:27,280 --> 00:49:28,730 small positive. 669 00:49:28,730 --> 00:49:34,210 So by adding on alpha*I, I push those away from 0. 670 00:49:34,210 --> 00:49:36,780 Over there we had the problem of eigenvalues 671 00:49:36,780 --> 00:49:41,000 being positive or negative, and by pushing in one direction 672 00:49:41,000 --> 00:49:44,680 we might just bring a negative eigenvalue right onto 0. 673 00:49:44,680 --> 00:49:47,720 That's not our problem in this upcoming lecture. 674 00:49:47,720 --> 00:49:56,250 In this upcoming lecture we know that nothing is negative, 675 00:49:56,250 --> 00:49:59,680 but we need to make it properly positive. 676 00:49:59,680 --> 00:50:00,590 OK. 677 00:50:00,590 --> 00:50:02,680 So that's what's coming. 678 00:50:02,680 --> 00:50:05,860 Can I just say: the application that I want to describe 679 00:50:05,860 --> 00:50:10,410 is the first one I've given in the life sciences, 680 00:50:10,410 --> 00:50:13,910 and I hope not the last. 681 00:50:13,910 --> 00:50:19,030 I don't know if any of you can suggest -- and if you can, 682 00:50:19,030 --> 00:50:22,120 an email will be very much appreciated -- 683 00:50:22,120 --> 00:50:27,020 applications of this whole framework in life sciences. 684 00:50:27,020 --> 00:50:29,180 So you'll see a particular application 685 00:50:29,180 --> 00:50:34,850 of this ill-posed problem in trying to determine gene 686 00:50:34,850 --> 00:50:38,450 structure from expressions. 687 00:50:38,450 --> 00:50:39,060 OK. 688 00:50:39,060 --> 00:50:42,980 So that's what's coming Friday. 689 00:50:42,980 --> 00:50:48,750 And I realize, here I'm talking about material 690 00:50:48,750 --> 00:50:54,050 that probably is not going to go in your projects, 691 00:50:54,050 --> 00:50:58,040 because of the timing. 692 00:50:58,040 --> 00:51:01,760 I'm thinking that your projects are 693 00:51:01,760 --> 00:51:06,790 more likely to build on what you had in project one, 694 00:51:06,790 --> 00:51:13,790 and the basic material about solving 695 00:51:13,790 --> 00:51:18,250 large systems and basic minimizations. 696 00:51:18,250 --> 00:51:18,890 OK. 697 00:51:18,890 --> 00:51:19,790 Thank you. 698 00:51:19,790 --> 00:51:21,530 See you Friday.