ME565 Lecture 11: Numerical Solution to Laplace's Equation in Matlab. Intro to Fourier Series

alright hello everybody okay so this is a pretty exciting lecture because now we're able to start actually computing solutions of partial differential equations in MATLAB which is awesome so just a show of hands how many of you have solved a PDE before in on a computer see - any more than two okay so this is going to hopefully be pretty interesting let's see a couple of logistical things so I just posted solutions to homework 1 and homework 2 then I realized homework twos not due until this afternoon so I took homework two down but it'll be due I'll post the solutions like sometime this weekend so if you absolutely need an extension on this homework you have to tell me now like immediately or else I'm not going to take any late homeworks because I'm going to post the solutions today so you can like see them I was going to post a homework three but I think I'm going to do it after some day because I don't want anyone to do anything except focus their energy for the Seahawks winning all right and I've also posted all of the lecture notes up to date and everything through today's lecture is fair game for the midterm okay and the midterm will be handed out on Wednesday by handed out I mean it'll be posted here on this website on Wednesday after class okay let's take home exam same rules as last time any questions at all no okay all right so this is pretty exciting because I'm going to start to get into a few really important concepts today that I hope are things that you take with you for a long time so in the last two lectures we solved we solve Laplace's equation analytically by hand right and it was kind of messy alright we got a bunch of signs and cinches and an infinite series which seems pretty nasty and so what we're going to do today is we're going to show that the expression we derived is actually pretty simple to implement in MATLAB so we'll be able to the cool thing is we'll be able to solve this heat equation on a rectangular domain with arbitrary arbitrary boundary conditions so these could be any nasty function that you like of temperature on the boundaries I'm also going to show you how to solve this Laplace's equation without any analytic formulas like without any cinches or sines we're just going to do it purely numerically just like we solve numerically ODS using OD 45 we're going to cook up schemes to solve these PDE s numerically okay what else am I going to do yeah and this is really important we're going to start introducing the concept of Fourier transforms today and this is all that I'm going to talk about next week our Fourier transforms and how they are useful in PD ES and how they're useful for representing functions like our boundary conditions in a nice way in terms of sine waves okay all right so we're going to solve Laplace's equation numerically on the computer and there's a couple things we can do so I'm going to give you like five or six different patience for solving this we're going to try them all out in MATLAB and we're going to see which one's work best okay so numeric solution until applause equation okay so just for my review one is Laplace's equation del squared u equals zero okay so this is a function of two variables this would be partial squared partial X plus partial squared partial Y of U equals zero okay and this is a linear operator alright this is a linear operator in U so if I have two solutions u 1 and u 2 and their sum is also a solution okay so the first technique is actually kind of interesting so what we're going to do is we're going to use the heat equation UT equals some constant del squared U and we're going to try to solve this iteratively we're going to start with some initial condition some initial guess and we're going to iterate it forward in time and we're actually going to let the heat equilibrate using the heat equation write the Laplace's equation is the solution as time goes to infinity of this equation with alpha positive constants okay so how can we solve this kind of a this kind of a PDE on a computer what are what are what general approach should we take to solving these things on a computer something like a Newton Rapson for which part for this part or for this part the UFT part okay okay so we have this U of T part and this del square du part we're going to need some way of representing both of these maybe I should ask a more basic question how are we even going to represent the solution in a computer sorry but what do I want out at the end I want a movie I want to see my heat distribution evolving in time to send equilibrium and what is that movie going to have its going to have frames in the movie and each frame is going to be a rectangle right with colors so what I want is you to be a matrix okay you is going to be a two dimensional array in MATLAB it's going to be a big array of numbers okay and I'm going to try to update that array of numbers using this essentially differential equation this is like an OD e for every single little grid cell in this array of numbers okay so I'm going to have this two-dimensional array of temperatures at every point in X and Y and I'm going to figure out how at each point x and y how does that little grid cells temperature change as a function of its neighboring grid cells okay so the first crudest thing we can do so this U is like this big kind of array looking thing right like it's a big matrix or you know super tall vector if you want to think about it like that so the crudest approximation to this time derivative what do we learn from 564 what's the what's the simplest first derivative we can take sorry I heard boiler what else like a forward difference so we're going to say that this is U at time t plus delta t minus U at time T divided by delta T equals alpha laplacian of u and i'm basically going to take an initial condition U at time this is U at time T I'm going to take some initial condition at time 0 and I want to step it to time one time two times three and so on so what I want is U at time T plus delta T as a function of the previous time step so that's equal to u plus Alpha Delta T laplacian of U at time T so this defines an iteration if I have some initial guess u naught u know u at x 0 then I apply this formula and I get U at time 1 then I apply it again I get U at time to U at time 3 dot dot dot dot and eventually I get u equilibrium which is the solution of Laplace's equation so this is this isn't anything you know super fancy we're just doing what we've done before and 564 to solve this Oh de ok I'm thinking of this equation as a really big Oh de with you know M by n ordinary differential equations so let's say that there are maybe L by H so I'm going to have L by H ordinary differential equations and this is the formula to iterate all of those ordinary differential equations forward using forward Euler so forward Euler is the crudest thing we can try but it's going to work which is kind of cool okay does this make sense okay so what's the next thing we have to do in this equation so we've figured out this term what's the next hard part yeah how do we figure out what this laplacian of u is so what is the laplacian of u okay partial squared u with respect to x squared plus partial squared with respect to y squared okay so at every single grid point you know at every IJ grid point here I and J are indices for my grid points they go from index you know 1 to L 1 to M do we know how to compute these partial derivatives yes okay how do we compute these partial derivatives so first of all there's a really easy command in MATLAB called del 2 so del 2 will compute the laplacian of my matrix okay so this is the easiest and most accurate but I don't really know what del 2 is doing under the hood so I'm going to tell you how we would do this by hand if you had to so how would you compute this derivative by hand if you needed to right you take finite differences in the x-direction to get the second derivative and X and finite differences in the y-direction to get your partial derivatives in Y and this is essentially this method is called a stencil that's kind of a cool idea so let's say that I have my domain with all of my data points right this is a big array you just write all these down let's say that this is my large rectangle rectangular domain and so the idea is if I want to compute the partial derivative the laplacian at this point I can essentially do that by considering finite differences with its neighbors ok so this is what we're going to do we have so let's call this the I J point and if we want to compute the laplacian of u at the IJ point we're going to use this five-point stencil meaning we're going to use the data from the five points around it from it and the four points surrounding it its nearest neighbors okay this is called a five point stencil okay and it's actually super duper easy to figure out what this derivative should be you can also figure out how accurate the derivative is you can cook up higher order accurate schemes using nine point stencils and eighteen point stencils or there's all kinds of fancier ways of doing this if you take a larger and larger window to compute these derivatives just like we did in in 1d functions okay and though the func the formula is pretty simple so this is I J this is I plus 1 J I minus 1 J I J plus 1 and I J minus 1 okay so what is my partial derivative what's my first derivative okay I'm just going to write it down so first I'm going to do my partial derivative and X okay so partial u partial X is U I plus 1 J minus u I J divided by Delta X that's this partial u partial X here I'm using I'm finding the first derivative using forward difference and then I'm going to subtract the first derivative using backwards difference and divide by Delta X minus UI J minus UI minus 1 J divided by Delta X divided by Delta X and this is the second partial derivative I'm basically just computing two first derivatives here and then I'm using I'm taking the first derivative of those first derivatives that's what this is and you can work it out that this is really I mean this is we've done this before this is just the the second order central difference that we derived in the last quarter ok you can do a similar thing for partial squared u partial Y squared that's basically the same stuff except instead of varying the I index the X direction we're going to be taking this along the Y slice by varying the J index so pretty straightforward we have U of I J plus 1 minus u I J divided by Delta Y minus u I J minus u I J minus 1 divided by Delta Y divided by Delta Y and so let's just get rid of this and write out what the solution is so I'm going to assume that Delta X equals Delta Y and I'm going to add these two things up to get the laplacian so my laplacian equals 1 over Delta x squared I'm assuming Delta Y is equal to Delta x1 over Delta x squared times okay what do I have I have a positive I have a positive value of this point that's here I have a positive value of this point a positive value of this point and a positive value of this point so these four points I have a positive value of U I plus 1 J plus u I minus 1 J plus UI J plus 1 plus U at I J minus 1 and then I have four negative copies of you in the middle I have four negative copies of my central points and index I J so minus for you at I J ok so this is a really simple formula if I had a function defined on this array if I had my temperature distribution u naught or u1 or u2 or whatever then I could compute its laplacian using this formula and I would just write a huge for loop through I all the eyes and all the J's and I would have the laplacian at every point this only works for interior points right this function has I plus ones and J plus one so it only works for the middle the middle points okay should we code this up and try it I think we have all the ingredients we need we know how to compute this iteration in time this is time and we know how to compute this laplacian of u if we have an array a rectangular grid of points okay so we're ready to code this up let's try it again this seems like I should go farther down let's try that for now okay okay that works for me okay so what we're going to do we're going to solve the same basic setup that we had derived an analytic solution for in last class we're going to have a rectangle and this rectangle does not have to be square but I'm just going to make it square because why not so it goes from 0 to L 0 to H and it's going to have 0 temperatures on the top and bottom it's going to have we could do anything this is completely general but I'm going to have it have a temperature fixed at 1 on the left and the right okay we could do something else if we wanted okay so I made a color map for the Seahawks just this morning on my bus ride over people thought I was crazy and this is what the iteration looks like in time so I'll just show you what the solution is and I will walk through what all the math does and what all of this means but what we have is an imposed temperature on the left and the right wall is equal to 1 okay and we see that for long times our temperature is 1 we have an imposed temperature of 0 on the top and the bottom which is also satisfied and as time is evolving we have our temperature kind of diffusing into the middle where they're meeting okay which is kind of what we expect to happen it looks so much better on my computer screen that it's kind of irritating let me see if I can make this larger okay I can I can make a figure and I can set my position 100 600 by 600 see if this works ah but I lost my Seahawks color let's try again okay and I bet I can make it 800 by 800 maybe seven by seven I'm killing all the lights to help it all okay so you can see these kind of ISO contours that are in the middle and if you could see a little bit better you would see that there are these kind of hyperbolic looking level sets of this function which makes sense because our solution has hyperbolic sines and cosines that doesn't help okay all right so that's what the solution looks like now let's go through and figure out how we actually coded this up so I told you there were a bunch of options right the first one was to use okay let's let's just start from scratch so the first thing I do is I create a 100 by 100 square and I say that my temperature distribution starts out being 0 okay and what I'm going to do is I'm going to iterate my heat equation two thousand times for two thousand DTS I'm going to iterate my heat equation so this should come into equilibrium and every single one of these iterations the first thing I'm going to do is enforce my boundary conditions whatever happened to mess up my boundary conditions I'm going to rewrite over that I'm going to say no the temperature is fixed at zero on the top and the bottom and it's fixed to be 1 on the left and the right okay so I'm manually fixing my boundary conditions by just overriding those numbers in the array with the zeros and ones now technically I don't think I have to do this I think I could bring this whole thing out and it would work ok good so I wrote this code correctly okay so I'm not ever changing those boundary conditions okay so inside my iteration I do the following thing the first thing I do is you know every ten iterations I just plot the thing so we're seeing this it's doing ten to ten DTS and then it plots ten more DTS and then it plots it gets really slow if you do it otherwise so this is just plotting and option one I'm using the del to command to compute my laplacian okay I told you that there's this function del - and it computes remember we had you K plus 1 equals UK plus Delta t alpha times the laplacian of U at time K so there's a function in MATLAB called del 2 and it computes the laplacian of u now I call this L of U because this is a linear operator in u L of U is usually what I use for linear operator so this is L of U okay that's this line here 34 so I compute my laplacian using the slick MATLAB command and now what I do is I take all of my interior points here in the middle I get a little marker so I take all of my interior points and I inert them forward using this this formula okay so I evolve I know that my boundary conditions aren't changing so I don't need to change them in time I'm fixing them and so all of the interior points once I've computed this laplacian I update them with this formula that makes sense I take my old solution plus some constant times the laplacian of u so I say all of the interior points from 2 to L minus 1 & 2 to H minus 1 all the interior points are going to equal their old solution plus some number times my laplacian so this makes sense yeah does this not make sense to anyone I can go through it more carefully right this is a big array and this is a big array and those correspond to you at my last time step and my laplacian of you okay now what if I made this constant delta T times alpha larger right maybe I want this to converge more quickly I want my heat equation to be more aggressive so I'm either going to have more thermal diffusion or I'm going to have a larger time constant those are the two ways I could do it so let's say I make this 2 times the laplacian of u what do we expect might happen to this equation if I make DT or alpha too large say again well I could converge faster or it could go unstable those are the two options it's either going to converge faster until it doesn't and then it will diverge to infinity so let's see what happens when I make this larger ok so notice now that this is essentially this is 10 to the 300 right that was a huge number going to infinity it's still diverging let me make this a little smaller okay so what we notice is that very quickly the numbers get extremely large and we have this weird instability growing okay and so this iteration becomes completely unstable I have this weird checker pattern anyone know why this is a checkered pattern starting from the corners any guesses at all so the corner is the most tricky point that's true these corners are where you have almost an infinite you know kind of derivative what else could possibly be causing this so this checker pattern corresponds to the highest frequency sine wave that you could possibly have on this discrete grid okay if the highest frequency sine wave in both x and y and it turns out that that wave number that frequency if I took this partial differential equation and Fourier transform it formed it I would still have a differential equation and the differential equation for the highest frequency wave number is the most unstable so it's growing the fastest and when I plot it all I see is that one ripple a super high frequency pattern so we're going to understand this a lot better next week when we look at Fourier analysis but we're starting to see sine waves here this checker pattern to me says really really high frequency sine wave is unstable okay so we're just going to tone this back down to a reasonable constant of one that's as good as we can do okay well you can do a little better but that's basically as good as we can do okay so that was one option is just to use the del to command to compute this laplacian but I spent all the time explaining this five-point stencil so I want to implement this five-point stencil so instead of just computing the laplacian using del two I'm going to compute the laplacian for every point I and J using the stencil we talked about okay so for all of the interior points in both x and y for i to 2l minus 1 J 2 to H minus 1 this laplacian at point I J is going to equal this formula that I derived and I'm just going to assume that Delta X equals 1 because I don't want to think about it I don't want to normalize anything or divided by Delta X's so I just have minus 4 of you in the middle plus a positive U at each of its neighbors and the plus sign that's it and then I do the exact same numerical iteration scheme right I do the same forward euler time step to iterate my u solution forward in time now notice that I need a smaller smaller integration factor here because this is actually less numerically stable than this del 2 del 2 is more accurate than this crappy five-point stencil and so this I have to take really small time steps or else my error blows up we'll look at that in a minute okay so now we're doing this option 1b we're using a five-point stencil basically we get the thing we want doesn't look like it's changing much on the screen but a lot of interesting stuff is happening so I encourage you to actually run this code on a computer screen ok does that make sense how we're essentially now all we're doing is computing this laplacian ourselves and we're doing kind of something silly that's not as good as what MATLAB is doing MATLAB is actually probably taking the Fourier transform and it's laplacian is much much much much more accurate than ours which is why it's so much more numerically stable another option is to do something called an 8-point stencil which instead of looking at just the four nearest neighbors I look at my eight nearest neighbors and basically I say that my middle point you know minus eight of my middle point plus one positive value at each of these eight neighboring points so this is an eight point stencil there's all kinds of stencils there's a whole pretty interesting theory about these stencils there's books written Wikipedia pages but for our purposes we get the same solution every time and I like Dell too because I don't have to type any for loops for loops are terrible in that lab so Dell is way way way way way way better okay please interrupt if there's any questions at this point this basically makes sense okay there's another thing we can do so I told you that option one was to start with the heat equation and iterate it forward alright another thing we can do that's actually pretty it makes a lot of sense is I can take my stencil formula and I can just try to make Laplace's equation true at every time step so option two is try to manually enforce del squared u equals zero at each iteration so for example in my five-point stencil we know that our laplacian looks like the follow we have laplacian of u you know at the IJ grade point equals well it equals something like 1 over delta x squared times minus 4 u at IJ plus u at all of its neighbors right these are my neighbors so I could make this manually equal zero in a pretty simple way I could just set this thing equal to zero and solve for my function at I J I could say u I J and I'm going to get rid of Delta x squared in my simulation all these Delta X's are 1 so this is just equal to 1/4 times the sum of all of my neighbors now what is that this formula look like what is 1/4 of the sum of the values of all of my neighbors the average I'm literally taking my neighbors and I'm figuring out what their values are and I'm averaging them and overriding you I give you the same thing with my eight point stencil I could have one over eight times the sum of my eight nearest neighbors okay this is a perfectly valid way to solve Laplace's equation and we could iterate this forward and this will converge towards the true solution okay so that's option two I'm going to code this one up now really glad I pre coded all of this okay so option two is we're going to forget the heat equation we're going to iteratively try to make Laplace's equation equal zero with our five-point stencil so what I'm doing here now is I'm this is a I shouldn't have called this L of U but now I have my average function of U this is the average of all of my neighbors it's just one quarter of the sum of all of my neighbors and now all of my interior points at every iteration I just set them equal to the average of their neighbors this seems like incredibly simple but it works okay totally works okay so we are solving PDEs numerically in MATLAB right now okay we have the tools to solve laplacian x' and DDT's if I put in some nonlinear stuff you could solve those two this is not the best way to solve PDS numerically but this definitely works okay this is a starting place we can actually do something now okay and we can start to we can start to have more interesting boundary conditions too so for example I'm going to have this new boundary condition where on my right-hand portion of my square I'm going to impose a sinusoidal temperature distribution okay just a sine 2pi 1 to L over L that doesn't work let's see I have to do something I I commented out all of my options so I'm actually not solving anything I did this a lot of times in my office so okay so this is what the heat distribution looks like as it evolves for a sine wave distribution so you have plus on the top - on the bottom and you can see that these are just getting kind of wider and wider and wider as they encroach on the domain but the zero temperature boundary conditions out here are essentially bleeding away a lot of that energy okay kind of cool you can do any temperature distribution you like no problem okay okay that's what I want to show you for now in MATLAB there's more I want to show you if there's time so that's how to code this up using finite differences this is the crudest thing we could do we did forward Euler and central difference for the laplacian very very crude but it does work we're going to figure out much better ways to solve this use Fourier transforms so even even in computers we would like to solve this in the Fourier frequency-domain it's much easier but first what I want to do is I want to show you in our example from last class we actually derived an analytic solution to the temperature distribution for an arbitrary function of Y on the right hand side right we had this huge nasty expression we spent two lectures finding it and so what I want to do is show you that you can also code this up in MATLAB and it's actually a useful quantity okay so I believe what we had was U of x and y equals the sum from N equals 1 to infinity of a n sine n PI / hy cinch n PI / H X and this a n was pretty nasty we had a n equals the sum horrible expression like 2 over H cinch blah blah blah and PI L over this is terrible integral 0 to H F of Y times a sine wave in Y dy ok now this looks terrible but it's actually not terrible these are all numbers H is a number L is a number and and PI are numbers 2 is a number of these are all numbers okay in MATLAB we can solve this super easy now let's try that example that we just did where we have this sine wave so f of y equals and it's not just a sine wave its sine of 2 pi Y over H it has exactly one period from zero to H so it's 2 pi Y over H that's my function of Y let's try to plug this into here and see if we got something that's simple okay so what's the first thing I need to do sorry I need to find these eggs these eggs are like basically if I have AIDS then I can plot this you I have this you so I need to find all these aids now what do we know about if I plug in a sine wave in here what do I know about how it integrates with other sine waves so this was a homework this is why I gave you that weird homework about integrating sine waves this will integrate to zero for all frequencies and except and equals to so a n not equal to R all zero a of n is equal to zero for all N that's not equal to 2 and a 2 is actually kind of simple so it's equal to 2 over H cinch shove 2 pi L over H this is just a number integral 0 to H of sine 2 pi over H squared dy I can integrate this in MATLAB and I'll get a number out and then that number I can multiply it by this single frequency in space X&Y and I get my function you kind of cool so let's do that yes all right okay so I have this l11 lecture 11 analytic laplacian same deal I have a grid from 0 to 100 0 to 100 or maybe 1 to 100 1 to 100 l was 100 H is 100 I create a big mesh grid for X&Y I start out with zero temperature distribution and what I'm going to do here is I'm analytically solving for the temperature distribution this is no iteration there's no heat equation no look discrete laplacian I have my boundary condition function this is my F of Y just 2 pi Y over H is my variable Y it's a vector of Y variables I have my constant a to which I can solve exactly how I wrote it here and then my u function is just a 2 times sine of that frequency and Y times Cinch of that frequency in X and I'm going to apply it super simple and I get my temperature distribution like that no iteration this is the infinite time Laplace temperature distribution it's also there's no numerical error here this is completely mean except for the Machine precision of cinches and sines this is exactly the right answer and it was really simple to solve in MATLAB all of these nasty integrals and functions are easy they're just you know sums and signs and cinches does this make sense how we did this ok I want you to look at this code if it doesn't this is important stuff what if I had a nastier boundary condition what if I had a nastier boundary condition that wasn't just a nice sign then what would we do what if F of Y was not a nice single sine wave then I would have to just integrate over at a bunch of times I would have to I would have to compute this a n a bunch of times I'd have a computer for N equals one two three four five up to maybe 100 and I'd have a hundred of these signs and cinches that I'm adding up but I know what the coefficients are and I know what the frequencies are so I can do that so the last example I'm going to do and by the way this is exactly like a numerically taking the Fourier transform so the last thing I'm going to do this is my hard analytic laplacian is now I'm going to solve this for the boundary condition we did earlier which is zero zero zero one this is a hard boundary condition for first sines and cosines so this is my function of Y it's just one at all of these hundred points my boundary condition is just a big vector of one hundred ones and now what I'm going to do is I'm going to take my temperature distribution set it initially equal to 0 's and I'm going to add up a 1 times sine of the first frequency since of the first frequency plus a 2 sine of the second frequency since of the second frequency and so on for K equals 1/2 some number let's try ok I'll show you the final answer I'm going to add all hundred of these Fourier modes and I get exactly the right temperature distribution that I would expect we could verify this using our other code fact why don't we just do that take our other code and so our other code iterates eventually to the same solution except this has a much a numerical error and it's kind of crappy but basically this is what it's supposed to look like is some kind of weird shape that's bleeding towards the left but we can get it almost instantaneously here by adding up these hundred sine waves and this is the infinite time distribution this is actually much more accurate so let's see what happens if I just have the first mode not terrible what about the first five minutes interesting I wish you could see this a little clearer it gets bright then dim then bright then dim then bright it's beating weird let's try the first ten modes that's interesting it gets bright Jen right dim bright dim bright dim bright dim looks like ten switches I'm having 10 modes if I did 20 it would be a little bit better until eventually I have enough modes that that alternating pattern is smaller than my dy okay so the basic idea here this is this example is completely designed to motivate Fourier transforms the idea is that I can have an arbitrary function completely nasty function and what I want to do is I want to find these a of N equal you know all of this stuff I want to find the coefficients of that function in a sine wave basis I want to decompose my function I want to rewrite my function as a sum of sine waves of low frequencies to high frequencies of various frequencies so this is essentially where the Fourier transform comes from is writing my function Y my function of Y as a sum of sine waves it's now a function of frequencies okay these AAS and if we can do that we can solve all of these pcs numerically very quickly very accurately and we can solve them for nasty geometries and much harder equations than just Laplace's equations okay so this basic idea is what we're going to use to introduce Fourier transforms on Monday that's everything I have have a great weekend and I'll see you on Monday

Channel: Steve Brunton

Views: 14,697

Rating: 4.936842 out of 5

Keywords: Engineering mathematics, Mathematics, Applied mathematics, Matlab, Partial differential equation, PDE, Heat Equation, Laplace’s Equation, Fourier Series

Id: QWtNLqF2Omc

Length: 48min 58sec (2938 seconds)

Published: Wed Apr 27 2016