Lec 17 MIT 18.085 Computational Science and Engineering I

so this is lecture 17 we're starting the second half of this series of lectures and today is kind of finite differences day first lecture about steady-state equations Laplace being the outstanding example we're in 2d of course we're doing partial differential equations and it's that step-up in dimension that makes these serious numerical problems so if I just look here so I'll describe the standard replacement the Laplace finite difference equation the five-point scheme has the basic accuracy we could upgrade it to a nine-point scheme with higher accuracy the issues are the solution of that system how much would it cost by direct elimination methods and those are by no means ruled out first cuz they're so simple and you know we're a lot of sub practice with Gaussian elimination but then the we can often save an order of magnitude in the computing time by some iterative method well-chosen so I'll speak about some of these I've just at this moment added the last one their conjugate gradients let me say that that's these are all described in section five point three of the of the text so section five point three in the chapter on numerical methods and including conjugate gradients are described there I'm not going to tackle a conjugate gradients in the lecture I've learned that there they're a sort of I mean that it's they're not difficult to program they're very powerful and quite an interesting history actually so this is the sort of area that I did my thesis in well more the afternoon lecture the initial value problems but there was a lot happening this was I don't want to say how many years ago but the the these issues were really coming to the front and conjugate gradients was an idea which my one of my advisors Heston s predict provided by the way I was highly interested to see in Eric missouri's lecture last week I don't know if you notice caught the name of the chap whose papers he had read which had transformed his life and transformed a lot of science teaching actually he's the outstanding science teacher in the country one might say anyway they were papers by David Hesse teenis at Arizona State son of my hero Magnus Heston s and by at the same time conjugate gradients were suggested by man named Blanche aus in Ireland another outstanding person anyway they're not an obvious idea they're built on three term recurrence relations and but oh what I was going to say about the history Sarris I hope it's okay if I talk about conjugate gradients without actually writing the formulas up because that would fill the board and we'd lose the lecture anyway it was they give the exact answer in in n steps oh I can't even tell you a little about the problem of course when anyway it was discovered that they didn't sort of they at the end of n steps numerical errors had spoiled it and they had to take more than n steps and it seemed like it wasn't competitive with elimination which also finishes right at a finite time so conjugate sort of as a direct message to get the answer at after a fixed number of steps didn't win even though it was a heck of a smart idea and then so some years passed and then a guy tried it I thought well maybe I should think of this as an iterative method do contra gradients get me close quickly because that's what we're looking for an iterative method we're not looking to get to the exact answer we're looking to get real close fast and the answer was they did and so they're back in favor very much conjugate gradient methods I can even tell you it with a picture what what the problem is suppose we're well you know that I'm gonna be solving K u equal B write this K is sort of like the interesting issue now what is it this is a K and this is a matrix K in 2d two dimensions and the question will be what does it look like and of course we hope it's going to be symmetric positive-definite and it will be I mean Laplace's equation is symmetric positive-definite and so we would expect the same out of the finite difference scheme okay so that's our job how to solve k you week will be when K is this two-dimensional matrix and I'll I can hardly resist writing down because that's like the main point here but to describe conjugate gradients maybe I could connect this do we remember that this is also a minimization problem here we have an equation of linear equation or a minimization a quadratic minimization right that those are the two that go together and the minimum we minimize some energy expression which is 1/2 of U transpose K u that gives the second degree terms whose derivative is K U and then oh I usually put def on that right inside don't like funny I knew as I wrote be that I wasn't something I'm so used to writing ax equal B but now I'm say you equal F is the right thing for our square symmetric positive-definite 1 so you remember that you subtract the Ute the linear term is U transpose F that was in an earlier lecture and probably not bad idea to see it again that that minimizing this quadratic is what gives that linear set of equations all good ok so what does this quadratic look like well tell me if K was the identity matrix what would the graph look like what would the graph of P of U so this is the u 1 u 2 multiple new directions of course I'm only on this blackboard I'll stick with 2 and now I'm gonna graph this polynomial and and it sort of looks like this yes it's a bowl remember K is positive definite so it's going upwards and and the minimum here what's that minimum value of U where's the minimum in this in this quadratic if I'm graphing this is a graph of the surface P of U it's kind of bowl shaped where does the minimum occur well according to us it occurs at K inverse F right so locating the minimum is the same as solving the system that's the thing we have a choice we could think of it as a minimization problem or we could think of it as linear system and like every step we take is intended to we could we could say we could discuss it in these in this language or that one and every step we take would make it simpler closer to the solution now I asked what happens if K is the identity for example not very interesting problem of course because we know the answer right away but if K is the identity matrix what is what's the shape of the graph what's the shape of that cross-section circle right it's a circle and we have no trouble I mean if it's a nice circle I mean it's fantastic of course one step will just go downhill and boom we're at the bottom no problem but now a real K has what sort of a cross-section well I mean what's the more J okay that's not just a scalar multiple of the identity or it's an ellipse it's an ellipse and the axes well of course it's an ellipsoid soit's oeid and you know in in in dimensions but but we really interested in the shortest and the longest axes that's what really matters and of course those are given by well they're related to lambda min and lambda max so the ratio of these axes is well as Bobby the condition number lambda max over lambda min or maybe it's a square root or something okay now why does the condition number make it harder ideas for you to see geometrically that we're trying to get to the bottom of a very narrow Valley and you narrower the valley little narrower extended valley the tougher the problem is because if I try to start downhill I let me let me kind of iterate do do what imagine that we're sort of like we've reached some point on the valley and we want to what happens is a typical iteration bounces back and forth instead you know you can't get you can't get an iterative method it's gonna go directly to the answer I mean that would be too much to ask for so what actually happens is you bounce back and forth across this valley and it takes a heck of a long time and this ratio is is the problem well so the idea of the conjugate gradient method is that unlike Newton's method which just keeps trying if Newton's method goes downhill like steepest descent that's the most obvious algorithm go down yeah wherever you are in this thing go straight down until you you see what will happen if you go downhill you'll you'll be on some kind of a parabola so you'll loop up again and you'll say okay that's the best Newton step stop there they go downhill from there it'll loop up again down here can you I don't know my I'm drawing this but I'm not really seeing it then you go down again it loops up again and then you get in this in this thing where you're struggling just back and forth you struggling to get to the center okay well so just in a word conjugate gradients is a very good idea that's chief point and the main point about it is that it works a little harder to choose the direction and it chooses a direction that really is not just one degree away from the direction it just took so it takes a fact it takes a perpendicular direction in in the K inner product and that's perfect so so you if this so that's why you get there exactly if this was really u1 u2 you'd get there in two steps to conjugate radius up shoot you take this steepest direction and then conjugate gradient you'd grind through the the formulas and it would give you a really the right direction to travel for the second step and you'd be there okay so that's that hasn't revealed the algorithm for kinds of great but it's mentioned this highly important name conjurer gradients and also mentioned the central idea and it gave us an idea of this another look at this connection between the two but my talk lecture now will really focus on this linear system okay so let's come to what's K that's our question what's the matrix K I'm ashamed to say that the book didn't do a very good job of describing the matrix K it discusses the solution of the system in in that section 5 3 and quite brilliantly but well forgive that but but it doesn't but the it's just a couple of pages of description of the matrix so I'll be doing a little more back in Section 5 1 so 5 I have to put plus five point 1 for the matrix it just talks a little about it so let's let's do the same now we're replacing the pluses equation and I'm gonna put minus signs in it why do I put minus signs in it because I want positive definite and the second derivative is negative so I put a minus and make it a positive thing okay and and I shouldn't say Laplace I should really say plus all because I've got a right-hand side of course the the source term might come in the boundary conditions no not in the interior either way okay so what's the what's the system this is it some region ok so I create a grid already is the limitation of finite differences if we're in a square when the differences are wonderful I mean the grid is just right if we're in a general curved region we're gonna have a lot of thinking to do about these points where the grid hits the boundary there'll be a headache you can just predict it and I don't I don't think we need a headache this afternoon so are these that would take more time and and away from the central issue of the of the lecture which is the interior which is nice and clean and simple let's suppose it's a square grid the distance H both ways H H is the H is Delta X and also Delta Y that's suppose of course I could deal with unequal grids okay then what's the replacement for the Laplace equation what's the discrete equation well we've talked about very long about how to replace a second derivative right it's it's with a minus 1/2 minus 1 right in the X direction a second difference replaces the second derivative now for you why why a second difference will replace the second Y derivative but it's a vertical derivative maybe I'll suggest that with a minus 1/2 minus 1 writing them vertically ok and then I put those together combine them and I have a 4 in the center and a minus 1 in the east west north and south nodes so that's the five-point scheme that's the little molecule five-point molecule and it gives I think second-order accuracy because these are centered so centered usually gives guarantees you that one step above the minimum accuracy so the second-order accuracy and that's very and it doesn't spread out very much and it's so it's a small molecule and and it's gonna be the same except at the boundary so we're certainly would like to the model case of course would be for when the region is a square then we'll have no boundary problems well we still have to think what are the boundary conditions so I guess I have to there will still be a question at the boundary and the answer to those questions will be depend on what boundary conditions are imposed if we're imposing Dirac lay conditions fixed conditions you remember like u equals zero then let's let me impose that here u equals zero Wow okay and then maybe at this end I'll impose so this is the ekhan in the x-direction maybe I'll impose in this direction the Noi Monde condition the free condition the derivative condition which would be D u the N in the normal direction but that's x equals zero say at this point okay and let me Oh let me impose u equals zero a top and bottom yeah that's a good so U is zero on all these sides U is zero on three sides but not just to give you some variety the fourth side is a is a free condition a condition on the flow rather than on the temperature a condition on W in other words oh do you remember that we have we expect to have the same U matrix a coming to V maybe I called it in or or e either way and going to W and to F and this will be a transpose of course and this will be in this case the identity and you remember the continuous version this was the dive this was the gradient so it was the gradient of U and C was one yes otherwise I'd multiply by something and then this was minus the divergence and that equals F that was that's that's our equation in our standard framework and now the problem is make it discrete well we've already said indicated how we're gonna make it discrete and then solve it okay so I could I don't think we want to spend a long what well okay look like I mean all right this is our pattern but but here I've directly said what are the what are the numbers so how many let's see what size will K be here I guess it's the number of in here we know what u is so those are not unknowns and here we have a boundary condition so we can relate that to this and we can relate that to this so actually I think it's only gonna be four by four and maybe we could even write down what K looks like this four by four K let me I'd have to number the nodes of course to know what I'm doing so let me number by rows Row 1 Row Row 1 and 2 these are the these are the numbers of the mesh points u u 1 u 2 u 3 u 4 okay so what's the replacement for the X derivatives for the X derivatives so so here here's my matrix K and it's gonna be 4 by 4 because it's gonna multiply u1 u2 u3 u4 okay all right let me look at you one here look at that point around that point in the x-direction oh well maybe we can put the x and y together probably we can did you do the whole four point thing at once so here we're in we're looking centered on point one where we would normally have we would normally have a four there let me let me here hang on a moment let me make you equal zero on all sides to make it all clear at the start and then I'll fix it it'll be the difference between our fixed fixed and our fake free going back to the very first lecture I'm guess I'm hoping that you're keeping in touch with that first lecture okay fixed fixed everything fixed I guess I should say fixed four times here because all four boundaries are fixed okay then what then we have around here we have four minus one is multiplying its neighbor minus 1 is multiplying the neighbor above and that's not involved you force is no there's no diagonal and similarly you two around you two there's a four and its neighbors are number one and number four and then around you three at you three there's a for its neighbors are four and two and this guy is neighbors are that and that okay that's my matrix K in this small example what are its properties suppose I had sorry did I mention did I mess up here this this row you're freaking out for youth three four and yes and it's no sorry it's did I get the wrong neighbor oh yeah oh I certainly did better come out symmetric right thank you thank you absolutely that's right now yeah yeah symmetric absolutely thanks catch that on camera now it's this is the right one okay one of the properties of that matrix well what do you figure symmetric we just made it that right positive definite you bet you bet okay it's coming from our old they transpose a we could figure out what a a is here I could or I'll put it in the notes for the lecture because we here's a case where we can go directly to the guy here so it's positive definite for sure and and by the way if I just showed you that matrix and you didn't know any thing about where it came from would you know that this was positive definite yes because what's so one test would be eigenvalues well usually I like in values it would be a messy test in fact I don't I've forgotten offhand what the eigenvalues of that matrix are so I'm sure they're quite neat but I know they're positive how can I tell that looking at that matrix that it has positive eigenvalues because the diagonal is dominant right the the I so here's a fact that all the eigenvalues of a matrix if the diagonal dominates the sum along the rest of the row and I'm willing to take absolute values of all those other guys so I would say the rest of the row gives me two but the diagonal is bigger it's four and if every diagonal is right is there then no question the the matrix is positive definite in fact I could even say something about where the eigenvalues are there's this sort of I've well in this case let me do this case first and then speak more generally I have a four times the identity so let me mark the four here that's that's where the eigenvalues of four I would be and then the eigenvalues of this off diagonal stuff how large could they be just just look at those minus ones do you want to make a guess on how large or small they probably have plus minus diagonal actually plus minus eigenvalues any guess on what the eigenvalues how big could they be if I if I ignore the four is now just look at the minus one one or two I think could be as big as two actually I can think of an eigenvector for eigenvalue 2 there can you see it I think Oh No maybe not yes I'm just going to take the all ones vector yeah I would say if I took the old ones vector I'd probably see that it was an eigenvector and it's eigen value would be 2 so from this for I could go like 2 in either direction from 2 to 6 the eigen values are between 2 & 6 I think it happens the 2 is an eigen value but the I know the are you measure between 2 & 6 so I know they're positive so I know the matrix is positive definite I know Gauss elimination will be great conjugate gradients once once and all these all these other ideas work well for a positive definite symmetric matrices yeah do you want me to take another second on that little location of eigenvalues because it's quite a neat result that is in the text somewhere it's gershgorin some Russian guy observed that if I took a matrix then he has a bunch of circles and gershgorin circles contain all the eigenvalues that's its so where are those circles of course well they're circles around the diagonal entry let's call that a I I so here are gonna be the circles this is the complex plane because I'm not George Corinne isn't only talking about symmetric matrices so here's the real axis here's the imaginary axis and here's let me suppose that suppose it's a real matrix so there's a I I that's the center of it's of circle number I now what's the radius then I've told you there circle the radius is that sum of stuff on that row other than this Center guy it's the absolute where you add up the absolute values of all the others and that gives you the radius so you add up the absolute values of all these entries it gives you a radius and suppose it's that beat and then there's another circle for Row one and two and all so we have a bunch of circles centered at like a 1 1 a 1 2 A 2 2 a that might have been a 3 3 maybe here's a 4 4 what could we learn if those were the gershgorin circles well I guess all we can learn what we learn immediately is that all eigenvalues are between what and what if they're real in other words all eigenvalues are in that so in one of them in one of those circles so certainly all eigenvalues are somewhere between that end and that end if they're real or in this funny shape if they're not real in a symmetric case we know they're real so we are only concerned with intervals of the real axis then could could this be could this matrix be Paulo here's an interesting question bring out the moment the key a little bit extra about gershgorin could this matrix be positive definite it couldn't and do you know why not sorry there's a circle there that doesn't even have a chance to give us a positive one and and the gershgorin result is that every does every circle have an eigenvalue I'm pretty darn sure that that would rule out positive definite yeah yeah yeah I'm pretty sure I've lost absolute confidence here because I don't not true that every circle must contain an eigenvalue well anyway your scoring is a it's kind of a little fun it's very elementary and it gives quick information about eigenvalues so in this case all the circles would be centered at four and they would all have radius two so they're all the same circle in this example and and symmetric matrix means that I only have to think of real eigenvalues they're between two and six okay good okay now I guess the next step should be to try to visualize what K becomes in other cases this was a small case what would K look like in the let's let's change that back to this other boundary condition well what's the effect of that yeah can you tell me what the effect of that is so what is this little molecule look like at a yeah maybe this is a good way to do it what happens to our molecule when when this guy when the boundary condition here this is not an unknown anymore it's not unknown anymore what's known is a derivative so it's known that this is the same as that that's all at least that's our approximation so approximation is that that's the same as that so what happens to our molecule it loses a atom Oh kind of cute but that's as being the same as that left only three there is that right I think that would be so I think that if we if next to one and three we have flow condition instead of a fixed condition then I think next to one and three this is only a three so next to one and three so let me look at number one so this had this had originally a four remember there was a four there and there were minus ones but this was a known so it didn't show up a minus one there and the minus one there but now this guy is not a known it's a it's only known that it's the same as that so that becomes a three and also number three becomes a three now is that matrix I was good to ask you this is that that matrix is still symmetric is it so positive definite yes what can you say about its eigenvalues what I can tell right away that all eigenvalues are larger than what one good their score in circle I don't believe one is an eigenvalue now but yeah now okay so but that's still fine to work with we totally have to deal with those we might have the same we might have zero slope here we might have zero slope here what would happen if we had a few more zero slopes than these fours start dropping the threes and even twos in some cases Oh tell me what happens if I have this condition at all a normal derivative zero at all all edges free so free free free free square here then what are the entries of K what are the diagonals in K so if it's free free free free or what's what's gonna be at for example here they're all down to two right they're all down to two now now I'm this is I'm not gonna allow that the name K anymore why not because it's like our B matrix from the past it's what all right it's now singular it's now singular physically it's singular because of course if I'm giving the derivatives on all the edges then I could only know the answer up to our constant right if I a constant would solve the the zero of the homogeneous equation the constant is in the null space of course a constants in the null space the vector 1 1 1 1 is in the null space of that matrix yeah so we don't you need something to hold the whole system down in other words if it was like a mattress I mean it's kind of a line of Springs is now jumped into a two-dimensional mattress is that the best word to pick from imagining Springs in both directions or cot maybe then then if there's nothing fixed at the edges of course it's completely out of control so I won't we won't we'll go back to the good one all right okay okay so the good one was something well the very simplest one was all fours okay no one more minute on this before I speak about iterative methods because we just have to visualize the equation jump up to a bigger system hundred by hundred say a hundred by a hundred interior unknowns so so hundred and one steps here hundred in the middle and a guy on the left and a guy on the right okay what's it what's the order of K so so we have now so this is like a thought experiment because I don't plan to write down the matrix so we have a hundred by a hundred interior use okay so that so the make a has size what ten to the fourth right hundred n squared and squared let me let me just call it and buy it n times in unknown so so the size is N squared okay now N squared by N squared now I'm in a minute I'm gonna have to ask you what are the numbers in K oh yeah I already have to begin to ask you what does k look like yeah could you tell me something about what K looks like I'll order those unknowns just as I did here a row at a time then the next row then next row than the next row but do you see the problem that we have a 2d problem with we have a two dimensional arrangement of nodes and we've got to put them in a line we've got a number of them from 1 to N squared and it's sort of the curse of dimension we were our matrix a row of our matrix is got to deal with not only the row of the thing but the guy above and below so anyway what what can you tell me about K here so big big matrix what's a typical entry in on the diagonal for still for sitting on the diagonal I'm somewhere in the middle of the matrix let's say for is a typical guy on the diagonal and what else is on that diagonal on that row what's on the row along with the 4 well it's just coming from the molecule is I'm in the middle here of I'm in the middle of this big square so I have the true molecule nope no boundary stuff is messing me up so what do I see then in on that row I see just to the left of this 4 I see a a minus 1 and just to the right I see a minus 1 so that's like tri diagonal then quite a distance how many how many spaces over to the guy who's in the the you that's in the row above yeah n like a hundred and a hundred in this direction sorry if I'm at if I'm right on the edge of a boundary if my if this if this guy was real new I mean I've drawn it sort of as if it didn't quite hit a boundary if if it was if my Ford my center node was next to a boundary then that might change to a3 or I might lose I would certainly lose the minus one but there there is a row that's sort of well into the middle of things now so you see the main idea of the matrix is we've got fours on the diagonal we've got tri diagonal but then we've got these diagonals which are far off that distance is exactly as you said and the bandwidth is and the matrix is a band of width M because it took n counts before we got to the end of the row started the new row and got to the node above or the node below this this is n steps to get from here to here okay so suppose we use elimination which we certainly can what will happen what will happen and in the in my message about this lecture in advance I created the matrix K for some boundary conditions and I think I cried I looked at the inverse just out of curiosity but I couldn't see anything special I should have looked at L&U and what tell me about L&U so if I do elimination if I do elimination to the fact the L&U factors look like this the triangular part looked like that and the upper triangular part looked like that that's the key question imagine doing elimination on this five diagonal matrix and as you've said you might have a break here where where a row ends and the next row starts there might be a boundary anyway but it's it's got five diagonals but they're not all minus ones I mean we have to hit the end here somewhere if they were all minus ones forever well of course then matrix would only be semi definite but the boundary condition saves the positive definiteness but come back to this point what about elimination because people certainly do elimination on this and what will be the problem the coffee what's the cost of elimination and where do we see it what's the cost elimination cost so size N squared then wit's end or half bandwidth you might say is that's the this is M this is N squared what's the what's the cost first of all for a matrix if we know if we know it's a banded matrix and we know the size is N squared and we know that's so that's the order of the matrix and then we know the bandwidth the general rule is order times bandwidth squared that's the that's the elimination for eliminate elimination count that's why I try diagonal matrix is so fantastic its bandwidth is only one I mean we're proportional to its proportionate there's a constant here we're proportional just to the size of the matrix for a tri diagonal now ours is not tri diagonal for sure and we get n squared times the bandwidth square time so that's end and we get into the fourth operations for elimination and and in three dimensions oh can you imagine three this is in 2d this is 2d for an N by n grid do you want to tell me what 3d would be for an n by n grid well what's the size of them 3d matrix this how many unknowns in 3d and cubed and what's the bandwidth and squared right because you you haven't finished that row until you've got normally to the one above but the one in the other dimension so it's N squared twice ATM seventh so now you're getting to like serious serious tax cut is and and what's what's coming back to what's happening at elimination as I do elimination the band is preserved I'll never have to do any elimination down there right the band will stay there and nothing will ever show up there but all these this big band of zeros between right between that minus one and the next one about 99 zeros that will all fill in that's the problem you get tremendous fill in if you just do gas and you see what happens when you subtract multiples of one row from another it starts filling in these otherwise zero things and of course the symmetric thing the L transpose has got the same and and you'll see it if you if you would just take the matrix in the in the message about this lecture and and do Lu which I should have done before I send it you would see that it fills in not so you could hope that some of the fill-in numbers are not too big oh and I can tell you right away the key iteration idea compare it with exact tell you so I'm gonna compare exact Lu with incomplete Lu so exact tell you is what so this is perfect to make this point in and this incomplete Lu Lu as I believe has is a favorite its people like to like to work with it so so if we don't have time to discuss the subtleties of over relaxation which was the hot so David Young wrote a thesis at Harvard in the 50s with Garrett Birkhoff in which he studied this over relaxation idea I haven't explained its it's a iterative algorithm and that was the big success for a while then there were other competitors conjugate gradients looked like it was dead got revived so it's in there and then a guy named van der Vorst in Holland came up with this incomplete oh you simple idea and so that's what I'll describe here because it's so in to ative and it you know it connects directly without any fancy formulas maybe you're seeing though from this from this brief account of the history a rather important point that everybody's got a shot at this thing like you don't have to be gauss to devise an iterative method which could be quite a good one I mean van der voice is a sharp guy but he's and Gauss and it's made a major difference so that's let me tell you what it is then okay so complete Lu this is this is elimination this is Gauss okay just do you do Lu you solve you factor in you solve Lu UE black right so this is our au equal it complete is a equals Lu of course and then how how is it you solve this system if you knew the factors and you could forget a you only need to remember Ellen you and how would you solve that system just reviewing basic idea of elimination how do you what what good does it do you to have done the work of elimination and I've got to Ellen you how would you solve that system there by L first right it's this you would solve for this unknown first say L C equals F and then C is this so they did solve u U equals C this is forward elimination this is back substitution forward elimination as it changes you as you go forward down the matrix you're changing the original F into some C so the C is what's left when you're ready to start back substitution back substitution is back up with the new matrix okay in other words two triangular systems that's fast the cost was finding L&U the cost was this factorization because you had you see why it's why this number might come in this N squared time this size times bandwidth squared why is it sighs x times W squared because this thing has this row has length W and you have W work W possible things to eliminate below so how much does it cost to do the B get to begin elimination on that first column you're taking multiples of this row of length W from these W rows so W squared to get this guy and then you're down to the next size and there's a W and it's got W and close to you W squared to get that column right and so on so it's size times W squared of course you and L have that same bandwidth W but back substitution would just be size time W naught squared and okay so what can you guess what incomplete Lu is seem so natural well of course not it's natural to think of it it's not clear that it would work what you do is you go forward with Ellen you but for entries that are small you set them to zero so that hopefully there's a lot of small stuff in this fill in and you say the heck with that so incomplete Lu is a is a proxy this is an L approximate you approximate in fact it equals this plus some error some stuff that you set to zero and there's a whole lot of freedom clearly you have to decide okay what's small what's not small okay so so what's how do you use this well iteratively so the idea now is I'll solve l approximate you approximate these are nice ones with with lots of zeros times you my n plus first guess well to me use a different letter I plus first guess is I'm going to put this part of the the e part on the other side as minus e you on the old guess and then I'll add in so now we're seeing iterative method and we're seeing one one iterative method in particular and let's suppose it converges does it converge to the right answer if it converges do we get the answer right sure if this approaches if you approach a some limit as I goes to infinity then this approaches you and this approaches you and and this approaches l approximate you approximate you and this is minus e U and this is the F and I bring that over and sure enough I have I have K u equal F yes good so now you see every iterative method there's two two requirements two conflicting goals here one goal is that this part should be easy to deal with we want this to be easy and the second goal is we want this iteration this convergence to be fast and those are in conflict because what's what's going on up there the convergence is immediate I'll put everything on the left side I did it all exactly right in one step I get the answer but that step isn't easy isn't quick it's the end for Thor and sevens now this step doing it this way it'll take some convergence in fact we aren't gonna run it to the limit either we want to get with him some tolerance of the answer that'll stop our iterations but do you see that it is now easier because these L approximate and you approximate are are trying you're there maybe maybe just just the three diagonals that we would have hoped for maybe I'll put those in yellow so so for the L approximate we might just keep those numbers thinking well those are probably the biggest or if we have to keep maybe one more because it's significant to have four diagonals but we don't have any diagonal we don't have a hundred diagonals NL or it and similarly mu it's easy it's got all the other diagonals in it but so he I could rewrite as a - L approximate you approximate you see that you see the idea the idea is to try to take as much of a as you can as much of cake I should be saying cash and I I've fallen in my trap again ferment of thinking of everything is really ok my litter I want to put as much of K I'm gonna put the heart of K on the left side without without making it hard to without making the each each iteration step hard because I'm gonna iterate maybe let's see if if N is a hundred I would maybe hope to do it in like ten iterations so there'll be ten iterations like square root of N or something this is what this is where you get in the calculation and and what I'll just end with the question what controls the number of iterations what can or whether or convergence at all you I think you're answered that question before we if we subtract this equation from this one we get the error equation so this was if we subtract this from this we would get the same matrix L app view app times the error at step I plus 1 equals minus e times the error at step I and the F so would cancel and I could put a I put I could bring all this I could put a I plus 1 as the inverse of this times the e every iteration is multiplying by some matrix which we never actually compute but but we know that's what we're doing the error is getting multiplied by that matrix and so what's the deal on convergence what would make that convergence fast if I if I'm repeating a multiplying over and over and over by this same matrix to bring the error down what's the requirement on that matrix the eigenvalues yeah the powers have to be getting smaller and we hope that the eigenvalues they're our best guide to those powers so the largest eigenvalue had better be below one or we are not we don't have a great idea here and we hope it's if it's down as low as 0.9 that's quite good Gauss elimination and and or sorry gauss seidel and over relaxation might produce in a typical problem my Gauss Seidel might produce an eigenvalue of maybe let's say 0.99 so it takes quite a few iterations to get down to something over relaxation might give an eigenvalue of 0.98 or even 9 7 that would be a tremendous improvement and incomplete Lu I don't know depends how close to complete it is so I'm I'm gonna stop there with the with the with this recap this is the type of matrix that we're looking at those are the some of the ideas that have gone into iterative methods plus conjugate gradients and finally let me mention one other idea if it really is a square then you would expect expect some special methods Fourier related methods FFT related methods to be successful and they are so there are special algorithms if you have a nice a square rectangle because because you can take Fourier transform them so those are the fastest of all by a longshot but of course they're only for squares and then you might try to build up your region out of squares so that's called domain decomposition and there's a whole there are conferences every year on domain decomposition where you take your domain and take nice pieces and then of course you have to make Corrections because if you saw up in one piece and hold the other piece constant and there's a correction to make so I could add domain decomposition to our list of ideas well so that's a sort of give you a sense of some of the thinking that's going on in in the numerical analysis community numerical PDEs really but it's matrix on this type of steady state problem and my second lecture yeah 3:30 will be quite different because the equations are very different we're moving forward in time there are time dependent problems and stability becomes the big issue here speed was the issue but there's the question is stability can we well and speed of course can we take big time steps and get where we want or will big time steps give us a blow-up situation how big could the time step be that's so that's the central question for us for lecture 18 okay let's take a break and and any questions and if anybody has any success with that curious matrix problem please please tell me about it I'm I'm looking for a a good good argument there okay