Chapter 9.5 - Finite Difference Methods: Overcoming Computational Difficulties

0:06 This is the last set of lectures on sort of really thinking about things you know computational difficulties we're going to face when we move to solving real problems. It is not enough to just get the solution. Your job is to get the solution most accurately and as fast as possible. If you use the generic ways of solving things, oftentimes you get yourself in a position where you're going to solve a x equal to b paying order and cubed. And if you're doing that professionally, you're going to get fired. You need to actually solve it in a smart way and and solve it in a way that's computationally tractable. So as we start thinking about these things, we also all think about difficulties that it can arise in in a variety of these sol complex PD systems that we want to address. Even not very complex PD systems, there's always issues we have to face and start thinking about handling ahead of time so that they don't sneak up on us and have problems.

1:02 In fact, one of the things that we want to do is already start doing the saying of the observations of troubles that are down the road for us even in terms of shallow water wave equations. So again, I love the shallow water wave equations because they already start even though it sort of doesn't seem that difficult, they produce very nice dynamics, there's already problems immediately that arise and start when you start a little bit deeper dive into what what am I actually doing with the solution technique. So remember this is just the poison equation but with periodic boundary conditions and if you don't think about it you think there's no problem with this just discretise solve. There's a huge problem with this and let me show it to you.

1:45 When you discretize this into ax b I'm going to ask you first the following question. If I had a solution to this let's call it f0 and if I had a constant to it. So here it goes. Here's a new solution. 5 0 plus constant. Is that still a solution to it? And in fact, if you plug this in, remember if I plug this in, this is two deriv, two deriv, which means it erases that constant C.

2:09 So the constant C has no impact on the solution here. And that constant C actually has no impact on the periodic boundary conditions either because that constant C is on both sides. In other words, what I'm telling you is there's an infinite number of solutions. Okay. So what does that mean? What is the consequence of that mean? Remember we took this del square s we we basically discretized it into this large ax= b system. What does it mean to have an infinite number of solutions? Well that means that the determinant of that matrix a was zero. And the one thing you always had to watch for, the only thing you had to be careful of in linear algebra when you're solving x= b is to make sure the determinant is not zero and the determinant is zero because you have an infinite number of solutions. In fact, the determinant of the solution being zero either meant you had no solution or an infinite number. And in this case, it's really easy to show there's an infinite number. And by the way, same thing happens if you have no flux boundary conditions in this as well. So there's two scenarios here where I have an infinite number of solutions and you've got a problem now of how you're going to solve this because when you solve this, you're trying to tell it solve a x equal to b computationally for something whose determinant is zero. This is a terrible idea and if you're not paying attention, you'll try to do this anyway, not realizing you're in this scenario.

3:33 Okay, so we want to figure out a fix to this. Well, here's a simple fix. If the determinant is zero, right, and that's what we have, then what we want to do is figure out how can I make the determinant not zero. So one easy way to do it. In other words, there's an arbitrary constant there. You saw this, right? I can add a constant to the solution, still a solution. So essentially, there's this arbitrary constant I have to pin down. So one way to do it is just simply say when I discretized on the diagonal let's say the first diagonal is supposed to be -4 I'll just change it to zero or three or 10 anything other than -4 and what I did by changing one equation of my system of equations is I took this matrix A which was determinant was zero and I made it now determinants not zero. So I I now it's much better conditioned and now I can solve the ax equal to b problem simply by this simple change. And notice this is in the corner of my domain. So has essentially what this is doing is pinning down a specific constant C.

4:41 That's it. So again you can do other kind of things that you want but this is a a very simple and quick and easy fix if you want it. Okay, that is these are the kind of prolevel tips you need when you solve things and also that you have to be aware of. The easiest way to check to see if this was a problem is to take the matrix A you built and do compute the condition number by doing com A and if it's big like 10 the 15 or 16 you know you're in the situ situation where the determinant is basically zero and so you need a way to fix that. Here's a simple way to fix it. Okay, so that's a really practical trick because now what you did is said instead of any constant C, it pins it down to one of them. Do I care which one? Because you might make make you uncomfortable. But notice when I solve for the stream function or vorticity, sigh is never alone. I'm always in those equations. It's always s of x or s of y. So the additional constant is always out of there. Know that when you take a derivative, the constant disappears. So doing this doesn't affect my solution at all at all. Okay? Partly because the only thing that shows up in those equations is sigh when it has derivatives.

5:59 Okay? Do you get around this if you do an FFT? Right? So that's a matrix that has a determinant zero. I still have the same issue if I do a fast poison solver like this with an FFT. In fact, if you go back to the FFT lecture on this, the way we did this is we took for transform in X, for transform in Y and that gives you a minus Kx^2 - K Y^2 and I can divide that through to get this representation here. So here is the solution of the free stream function in the for domain and all I have to do is invert this. So what's the problem? The problem is what happens when your kx and ky are zero?

6:38 Then you got a divide by zero problem. So the equivalent of a determinant zero problem in ax equal to b is a divide by zero problem in the fastfor transform. That's exactly what you have here too. And we need a fix. Otherwise as you write run your code you're having it divide by you're trying to make your code divide by zero, right? Which is a terrible idea. Okay. So how do you fix it? Same thing. You know, it's as easy as a fix as we did in the lelass equation with the matrix A, we just took one cell, the neg4, change it into something else. same thing we could do here.

7:19 Instead of the the zero wave number k 0, kx0 or k y0, in other words, the one that's creating the divide by zero problem, instead of being zero, which it is when you build it, you make it 10 - 6, something like this, something small. So, essentially, you just added a pertabbation. But all of a sudden this becomes not a divide by zero problem anymore and you've and it's a perturbatively close. Why 10 - 6? You know this is only be 10 - 12 when it shows up here. Your time stepping algorithm is usually 10 - 6. So this is going to make no impact practically on the solution. And what this does is again just like before it determines the constant C in the stream function making this wellpost.

8:05 Okay. So in other words that problem of an infinite number of solutions had to be handled if you're using the discretization of ax equal to b or in the 4 a domain and we found ways to handle it in both ways and their manifestation of a problem is it divide by zero here and in the other case it's the determinant being zero. super practical things to think about and you always want to understand your problem before you start solving it so that you don't think you have something that's correct when it's just wrong.

8:38 Okay, final things we want to do is talk about derivative matrices. So far we've talked about the A matrix before which was two derivatives X and two deriv. But we really want to do, you know, what showed up in our stream function equations and vorticity and stream function pairings is derivatives in x and derivatives in y. And so here's the derivative in x here, right? So if we take the derivative of the vorticity in x, this is just our essentially one derivative formula. This is a simple second order difference formula. But what we have to do is organize this also again like we did before. Let's write in our notation. Go to the mn position. This is the way it looks like.

9:19 And then this becomes a system of equations. And as we organize our vorticity, all of this just becomes if we arrange the vorticity this way, it's just a matrix multiply. So there's a B matrix that we've constructed here. And the B matrix heading this matrix here takes a derivative on X. And if I wanted to do that too with the Y direction, in other words, how do I do a derivative on Y? Let's call that a C matrix. So B takes this vector here produces a first derivative in X and the C matrix we could construct would produce a Y derivative of just with a matrix multiply as well. So what this allows you to do is once you've constructed this and by the way here's what the B matrix actually looks like.

10:07 It's n square by n squ matrix. And so you have, you know, figure out how to construct this. That's part of your homework to do is to figure out how did this actually result given that my data was stacked like this. It's not that hard to work out. It's a lot like building that a matrix. You have to go through that process once at least to make sure you understand the A matrix, which is two derivatives X and Y.

10:33 Understand the B matrix, which is derivative X. And then there's a C matrix which was one derivative in Y which looks different than this one as well. But once you have those matrices built then computing something like this this poison bracket between these derivatives just matrix multiplies all the way across. It's as simple as that. If you mess up anything in the ABC matrix your solutions are just wrong. So getting accurate ABC matrices is coming back to again the organizational part you're responsible for that is your responsibility that is entirely your code setup is how do I construct ABC to do the derivatives because once I have them everything else is set up now to go to time stepping and that's what we want to really engage in next is how do I do time stepping because now I have all the matrices set up to go about the process of solving that into in the future. So these are just practical chips and how we handle things. But all of these things have to be handled carefully and thoughtfully and you have to be aware of them. You can't just blindly solve these equations. These equations often will have things you know again when we solve that stream function with periodic boundaries. You have to understand that if there's an infinite number of solutions, your numerical code is going to break unless you take care of it.

11:54 So you can't just simply check your brain at the door. You've got to engage fully in the problem. Understand what are the things that you have to take care of before you just start computing away because once you start computing and get the algorithm right, it's really very fast from there. But you have to get everything right. So these are just lessons and how to overcome the computational difficulties that you're going to face in all kinds of problems all the time. And so it's your job to really kind of understand how to address them and the and the possibilities of addressing as you go forward.