Chapter 9.1 - Finite Difference Methods and Numerical Discretization

0:05 We now want to move to talking about finite difference methods especially in context of solving partial differential equations. In other words, differential partial differential equations are essentially differential equations in multiple variables. And so what we have to think about is that we don't just solve for a time evolution. and it now depends upon space and time or it could depend on x and y and t or just x and y. In other partial differential equations give us a much ripp richer representation of what we have to really work with in physics-based models. So let's go back historically what this means. So we go all the way back to calculus and when Newton and Leeness both developed the idea of calculus, it was a way to compute this thing we called the derivative and its anti-derivative called the integral. So this gave us a new language of representing physics models. Up to that point, what we had was geometry and algebra.

1:07 And this really changed everything because once you have calculus, you have an entire new language to talk about physics. And in fact, one way to think about it is all of physics and engineering since has been really defining relationships among derivatives. So let me just give you very specifically what I mean by this. Here's Maxwell's equations. What this is is relationships among derivatives in three spatial dimensions and time. So what we have here is the electric field coupled to the magnetic inductance B which is related to the magnetic field itself and we have time dynamics here.

1:44 This is the rate of change of how the magnetic indues change with time and how the electric field is changing in time. But how they're changing in time is dictated by how things are changing in space with these curls and so forth. So this is one way to represent electromagnetic phenomena which is something that we have been working on for quite a long time because it dictates for instance all of our wireless communications all of our optics all of our laser physics here are the governing equations relationships among derivatives. Schroinger equation similarly is relationships among time and space derivatives and this dictates how we think about quantum mechanics.

2:27 Typically, we have a time evolution coupled to spatial evolution in 1D, 2D, 3D, depending on how we set up potentials. Here's just the 1D timeindependent relationship that we have or fluid flows. Here's the Navier Stokes equation. So, this is how we talk about how the evolution of fluids happens with let's say with incompressibility. So again all of these highlight that if we're going to model physics and engineering models correctly we have to understand how do these relationship among derivatives dictate how we're going to get solution techniques. Most of what we're going to do is start thinking about time stepping algorithms.

3:11 We're going to march forward in time. But of course that's not everything we have to do. There are many of these equations that where we have what are called elliptic PTEES where there's no time involved but there multiple spatial dimensions in which we have to compute steadystate solutions. Okay. So how do we do all of this? All of it's going to come down to discretization. All those models assume continuous variables. And what we're going to do and we've already seen this when we do time stepping for differential equations is that we're going to discretize. So discretization now in space and time becomes everything we need to start marching ourselves into the future with solutions. But everything's going to come down to discretization. So I want to spend a lot of time talking about this and setting up our framework for discretization because that's what we're going to have to work with as much as we can. Okay. So the context that I'm actually going to introduce this is a fairly complex system but it's just relationships among derivatives and this is what's called 2D shallow water equations. What we have here is the evolution of the vorticity where this poison bracket is just partial derivative of x of the stream function s which is related to the vorticity through this lelass equation here. Okay.

4:32 So we take the leloian of the stream function is equal to the vorticity. So the stream function itself we have to solve for and once we solve for it we can take the x derivative of it partial and multiply it by the y derivative partial y derivative of the vorticity minus the partial with respect to y of the stream function partial with respect to x with vorticity. So these go equations can be derived from first principles models. And you see this is our first PTE. We want to start thinking about how do we come up with a structure to solve it. So at first what we're going to do is just talk about what do we have to do here? What's going to be the consequences of discretization? And the reason I like this model in particular is that it's not just about discretise and solve.

5:21 There are a lot of things we have to start thinking about beforehand that can come up as problems that we have to think about ahead of time in order to even solve this set of equations here. So, we're going to walk through some of these and you're going to start seeing that mostly what your job is is to be an organizer of the discretization itself because once it's discretized and arranged, the actual solution techniques are going to be very simple as I'll show you.

5:50 Okay, so the hard part is the discretization. I'm going to talk about what else is required in this. First, we have to have an initial condition and what we're typically given in the shallow water equations is an initial vorticity. So, I give you an initial vorticity and you ask the question, how does evolve in time for here, we're going to assume periodic boundary conditions for both the streams function and the vorticity. So in other words, if we're solving on a domain in x from negative L to L and in YL L to L, that this is what it means to have periodic boundary conditions. So whatever goes off on the left side comes back in on the right. Okay? So these shallow water wave equations were typically derived for shallow water across a globe across the earth. And so the periodic boundary conditions make some sense because it would just go around the globe in some sense. Okay? So we're going to assume that as our boundary conditions and that is our initial condition. So the first thing we want to talk about is what is a reasonable algorithm if we're going to solve these partial differential equations. So the first thing is I give you an initial vorticity. And so when we come to the algorithm your first step is you now have to calculate the stream function. And the stream function here satisfies the equation. It's time independent. So what we're going to show is this is just a giant ax equal to b.

7:11 So the first thing we're going to learn is how do you solve how do you start thinking about x ax equal to b at scale that comes out immediately from discretization itself. Okay. But once we discretize this and find solve for c for the stream function then we can just treat the entire rest of the evolution of the vorticity. We're going to turn this into a differential equation. So that's mainly what we're going to do when we solve PTEES. The objective is either analytically, if you look at a lot of analytic techniques for solving PTEs, it's about taking the PD, turning it OD, turning the OD into some algorithm from there. And often what we did in OD's is try to turn it into algebra, solve it, and come back up.

7:55 Same thing here, except we're going to turn it into ODEs and we're going to march it forward with a time stepping scheme many of like the ones we've considered already and with numerical methods itself. Okay, with with OD solvers, so once you take a step forward into the future, now you loop it. So you take one step forward, you have to update that stream function. So the vorticity changes, which means you have to update the stream function, take another step, update the stream function, take a so forth. This is the loop that you want to be in. So what I want to do in this first set of lectures is just really talk about what's involved in this process of discretization. We're not going to yet solve this, but we are going to start thinking about what are the things that are going to arise for us. How what does what's the math that arises out of this and how do we start formulating the problem and organizing the problem? And part of what I wanted to do is show this in a in a context of scenario where you start to understand clearly that a lot of your job in the discretization pro problem is to be the manager of how you actually do this discretization.

9:01 So let's start off with the losses equation. Here it is. It's two deriv x two der in y is equals omega. So that's a partial differential equation. There's no time evolution in here. And the first thing we could do to try to solve this is just simply say well what we learn from from computing derivatives is finite difference methods. And finite difference methods for a second derivative. The way you compute a second derivatives is use a point in front minus twice the current minus point in behind divided by delta x^ squ. So this is an x derivative. So in the x direction I use a point in front point behind minus two of that point. This is a second order finite difference approximations of the second derivative.

9:44 I can do the same thing for the y derivative except now my point in front my point behind are in the y direction divided by delta y^2. So this here represents this equation in discretized form. So what you're doing is you're using your neighboring points to compute your derivatives. The x direction using its neighbors in x and the second d y uses it neighbors in y. So another way to visualize this is through this grid right here. So when you commit this is a 2D problem. It illustrates starts illustrating to you the management you have to do for this. So first of all in this x and y domain you're going to chop it up into a finite number of points in x which is along this direction finite number of points in y in this direction.

10:33 Right? So I've labeled these here as m that's the total number of point or actually in this direction we're going n and this is m. What I'm going to actually start off by doing is just say I chop it up equally in x and y. So capital n and capital m will be the same. The delta x and the delta y's will be the same. But in either case, you start thinking about how do I represent each one of these grid points. Okay, a collocation point. Here it is. If I'm at one position, let's say the nth position in x, the mth position in y, then when I compute derivatives, I use my neighbors.

11:11 Okay, so this is the point one in front, one behind, one in front, one behind. Okay. And so this organization becomes really important when we start to get to practical solution techniques because the first implication of this is that when we solve that string function equation, we're going to get a giant ax equal to b problem to solve. Okay. So let's go ahead and represent this as our point. S of n mm will tell you which x point you're at and which y point you're at. Okay, so that's kind of how we represent this shorthand.

11:52 And now what we're going to do is if we take delta x and delta y to be little delta to be the same. Then this here is the second derivative in x second derivative in y. Remember the second derivative in x takes a point in front point behind which is these two here. Second der of y is point in front behind but now in the y direction. and then you have the -2 of the term you're on and that gives you two -2s from the x direct from the x d and y ds you get a -4. So this here is essentially the equation we have to solve for each n and m. If I organize all of them at every single grid point together, what I end up with, oh, by the way, here are the boundary conditions. At the boundary, since it's periodic, the first point is equal to the last point in X and in Y. And if I organize all these together, I get a giant system of equations. For every single point I go to, I use my neighbors to compute the derivative. And what this represent is my collection of governing of equations I have to solve for. So essentially what I have is a system of equations that's n by n^2 by n^2 that is an ax= to b. This is what we actually end up getting here. So what is a? So s itself is going to be the vector that we arrange and a itself is this. This is my matrix that tells me on the diagonal it's minus4 that comes from the second derivative x and y each will have a -2.

13:26 So that's -4. And then I use my neighboring points point front point behind for X point front point front point in front point behind and Y and it creates this matrix A. And this is a really important matrix that you should work through to make sure you can compute this and know how it's structured. The organization of the data is here. Here is my stream function and how I've arranged the data into a flattened vector. In other words, I've taken my stream function as two dimensional and I've written it as a as a vector. Same thing with my vorticity. So these are the two quantities I have. I'm advancing the vorticity in time. But every time I take a step in time, I have to also advance the stream function. And remember when we go to a solver where we're time stepping, what we step into the future are vectors. So the important thing is is how do I organize my vector and once I've organized that vector what does it mean to take two deriv x and y it means multiplication by this matrix this matrix here is the lelassian in this coordinate system that's given by where is it here okay anyway u sorry I wanted to show you the picture yes this picture here So this coordinate system here where I've taken all my data and I've decided how to stack up all these values of the stream function, okay, and the vorticity in this 2D domain. I'm going to turn it into a vector. And what I end up getting here is this set of equations here. This is ax equal to b. Now it's hard to work.

15:10 You just have to actually work through this a little bit. Take some time. Do this calculation. Make sure you can reconstruct this matrix. and this matrix. Now when you multiply a vector times this matrix here what you're getting is two deriv x 2 deriv y. Okay that's a really important concept to come in. This is a natural consequence of discretization. Under discretization, taking derivatives becomes matrix multiplication and arranging your data becomes your number one job because once you have this then you just solve a x equal to b. So that part's easy. The organization of the data and keeping it straight in your head that becomes the most important part of what you have to handle. Presumably you know how to solve a x equal b or certainly computationally it's easy to do. So then this is the hard part is getting this out. So spend some time understanding how to construct this. So once you have that now we have this first part done. So essentially what this is is ax equal to b that del square that lelassian is a matrix multiply. And the same thing is going to happen here when you start to advance it in time is that when you look at now the rest of the problem now that I have my stream function and vorticity I've now got to march forward in time. And of course now we have a a new complication which is before we were just trying to solve for the for the stream function.

16:38 Now we have to advance it in time. And notice that when we write this all out it's again a bunch of derivatives everywhere. So let's look at this term here. del^ squ times this the the vorticity del squar is a matrix multiply that matrix A that we constructed all I have to do is multiply the vorticity by a and I have two deriv x 2 deriv y and so what's left for us to do is also compute matrices what does it mean to take what's the matrix that takes one derivative in x what's the matrix that takes one derivative in y so under discretization what you have defined is first of all under the organization of the data you made.

17:23 Taking derivatives is going to just be matrix multiplication. So that part's easy. The organizational part is what you have to focus on in this step two to advance it forward in time because this is what the discretization looks like overall. And this is using just second order schemes. And this is what I'm trying to show you here is this is the implications of discretization. Once you commit to discretizing, you have to commit to organizing your data because all the derivatives are going to be matrix multiplies throughout. And what you end up doing is you advance your solution in something like this. So you start off with your solution at time t equals z. And when you take a step forward in time, each point in the future depends not just on the current point but on the neighbors.

18:11 So when you look at this, this is your stencil. Use a second order scheme which means you use your neighbors in X and Y. And so you get five points that update your next vorticity. And then if I update the vorticity again requires its neighbors again there. So you have to update all take all these points to update this entire set of points on this plane which then updates this plane and so forth as you march forward. We're going to talk a lot about time stepping itself which has its own complications regarding especially stability and accuracy that we have to do this way. But right now we want to just say like yeah but what how do I even solve for instance my string function which has no time in it and that's going to be a large ax= b problem and all the derivatives that I have to compute for these relationship among derivatives become matrix multiplied. So we have to figure out how to do that in a practical way. And by the way this is what the solution evolution of this looks like. Here are some just tiff time snapshots of this 2D shallow water wave equations. one I've given you is an elliptic vorticity.

19:18 And you can see what it does. It spins into sort of this vortex-l like structure as it propagates along. And of course, the vortex starts to die away because there's a little bit of diffusion in that problem. But it's kind of a nice hurricane looking behavior, which is exactly what it should look like because it's 2D shallow water waves, which we're trying to map or model surface waves of the ocean. Okay. So that's where we end up right now. And again, part of all of this is starting down this road of looking at finance different methods. And this first let's say set of lectures are really going to start to dig down on what are the implications. What are the things you have to start thinking about before you even start stepping off into the future. There are a lot of organizational issues you have to take care of and you have to start thinking about the fact that now you're going to get into very large scale problems because the axe a b problem we're talking about is now a very large scale system. We have to think about how to solve that in efficient ways which we will do as we continue on.