Flash and JavaScript are required for this feature.
Download the video from Internet Archive.
Description: This session introduces finite volume methods in two dimensions and Eigenvalue stability, then reviews the advantages and disadvantages of the methods covered thus far in the course before the midterm exam.
Instructor: Qiqi Wang
The recording quality of this video is the best available from the source.
NOTE: Videos for the next two classes, Sessions 13–14, are not available.
Session 12: Numerical Metho...
The following content is provided under a Creative Commons license. Your support will help MIT OpenCourseWare continue to offer high quality educational resources for free. To make a donation or to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.
QIQI WANG: So today, I kind of prepared not for the whole lecture, but a little bit short of that. So I really expect you to ask questions on this material. It's supposed to be a lecture that helps you review the material we have already covered and prepare you for the midterm.
So instead of me just going mechanically through the material, I want to ask you to initiate what do you think is more confusing, or you'd like me to clarify again, and things like that. If you feel something is confusing, it's probably confusing for the whole class. So please raise it so that I can spend more time on it.
OK, so before I do that, I just first want to finish up the finite volume scheme we have been working on for the last two lectures. We have been discussing finite volume schemes in one dimension. But applying the same concept to two dimensions or even three dimensional is surprisingly straightforward. So let me first give out what we did in 1D.
We started out in finite volume schemes-- we started out with the integral form of the differential equation. We started out with d dt of the integral of a conserved quantity inside a control volume, omega, it is equal to minus of the flux out of the control volume. The minus sign is because the time derivative is really the flux into the control volume. But the definition of our n, normally, is usually out of the control volume. So we have a minus sign to reverse the normal so that it points into the control volume.
Times the flux as a function of rho, ds. So this is like the integral form of the conservation law, right? And I can write the same thing in 1D that is the specific case of omega being just an interval. Now, this is the flux at the left minus flux at the right. Or it's really rho at left and rho at right.
Right, so in 1D, what we did was we set rho bar of k is equal to the size of the control volume of k. So it's the integral of the control volume divided by the size of the control volume. And basically, by plugging in the definition into the integral form of the conversation law, what we get is 1 over the size of the control volume, the flux at the left side of the control volume-- I'm just going to say F of k minus 1/2 minus F of k plus 1/2. That is what we use to denote the flux at the cell interfaces.
Now, let's go back to 2D, or the multi-dimensional form of the conservation law. We're going to define rho bar of a control volume, k, being the same thing. It is really the integral over the k-th control volume, rho dx, divided by-- in 2D, that is the area of that control volume. All right, let me draw a typical mesh in two dimensions.
So let's say this is a triangular mesh in 2D. And this is like the k-th control volume. Right, for example, this is part of a mesh in two dimensions. Right? And we are computing-- we found the volume scheme. We are tracking the average of the solution inside that control volume.
So this is my omega k. That's the control volume. And the area of that control volume is my Ak. And rho bar of k is the cell average there. So what is the time derivative of that cell average value? Because the area of the control volume does not change. It is equal to 1 over the area times the time derivative of the volume integral, right?
And the time derivative of the volume integral, by the integral form of the conservation law, which you can get by applying divergence theorem, right? We did that last lecture by applying divergence theorem to the differential form of the equation to get this integral form. And just by plugging in, we get 1 over Ak times the integral over the boundary of the control volume. We get a minus sign here. And the outward normal dotted with the flux, ds.
Now, let's look at what this is on the control volume. So we found the integral over the boundary of the volume. If the control volume is this angle, what are the boundaries? We have [INAUDIBLE], right?
Now, we can express that as minus 1 over Ak of summation of ni dotted with Fi times the length of each of the sides. So i is, I'm going to say, like the sides of k. So Sk is the sides, the three sides of the triangle, all the boundaries of this triangle, k. So Sk is a set of these integrals that gives you these three boundaries.
And the i the normal of each edge. So for example, in this direction, so this is the ni's pointing in this direction. These are the ni's. And the graph, the Fi is the flux on these edges, which we have gone through, again, approximately, that by applying the graph function of these two neighboring cells.
It's the same thing as we did in 1D. In 1D, we approximate the flux at the cell interface by looking at the value of the function at the neighboring two cells. In 2D, it's the same thing. We approximate the flux of the cell interfaces, so yield between two cells, by looking at the flux, by looking at the solution of the cell averages over the different sides.
If we have to do [INAUDIBLE], we do [INAUDIBLE]. We look at which direction is the [INAUDIBLE] over this edge. Is it from this side or this side? And choose the value on either side of this edge. If we go to three dimensions, it's the same. Instead of edges, we have spaces between control volumes. That is the case where we really become looking at interfaces, interfaces at each place in three dimensions instead of an edge in 2D or a [INAUDIBLE]. Now, we approximate using a numerical flux that is rho bar of k, rho bar of-- I'm just going to say neighbor. And then, we're done, right? We have approximated the time derivative.
So you know the [INAUDIBLE] of the cell average in the cell k. And the function of the cell average in k and the cell average in the neighbors of the k-th control volume. [INAUDIBLE] other point is the area, just like normal, the length of the k-th [INAUDIBLE], the n-th normal in this space. And the length of the Ses are all known quantities from the mesh.
So this is [INAUDIBLE] flux, a function of the cell averages. So essentially, we have turned a PDE into an ODE. What is that ODE? That ODE is d dt of rho bar 1, rho bar 2, et cetera, to rho bar of the last control volume is equal to some joint function of all these rho bars. OK, so this n is the function that connects the time derivative of one cell average to the value of a cell average of that control volume and its neighboring control volume.
Is it clear? How we do the same thing in 2D? There is a lot-- there is going to be a lot of bookkeeping just because of the mesh and for every cell, we need to be able to find its neighboring cells, and find which edge is connected to which cell, and who is whose neighbor. Conceptually, this is exactly what we do.
OK, so we have started finding the difference, finding the volume, and all these methods are methods that are intended to turn a partial differential equation into an enormous ordinary differential equation, right? And once we turn a PDE into an ODE, we can apply exactly what we have been doing in ODEs to study these PDEs, to discretize, solve, and study the behavior of these PDEs.
OK, so this is going back to the review part. I'm going to really both the rule of the mesh and the outcomes and put particular emphasis on several of the points. The first point is order of accuracy. So the order of accuracy in an ODE scheme is found by looking at the equation, du dt equal to f of u. And we discretize this into some operator. We discretize the d dt into a finite difference using a finite difference operator.
So we discretize the d dt into some kind of a delta over delta t. I mean, this is going to be finite difference operator, operating on the u. Now, this is not going to be-- if you know the analytical solution, if you know what u is, for example, if f of u is equal to lambda u, then you know what u is, right? And you plug this into the finite difference operator. This is not going to be exactly equal to 0. Right? This is not going to be equal to 0.
OK, so for example, if you're looking at Forward Euler, if you're looking at Forward Euler, u of k plus 1 minus u of k divided by delta t, if you subtract by f of u k, it is not going to be equal to 0. If instead of u k you plug in the [INAUDIBLE] solution, you're going to get 0. But if you plug in the real, analytical solution, you're not going to get 0. And the order of that approximation is going to be the local order of accuracy. So if this is equal to O delta t to the k-th power, then k is the local order of accuracy.
Now, here, the usual confusion that happens here is that when you do the truncation error analysis in another way, if you do the truncation error analysis by writing down the tau equal to u of k plus 1, if you plug in the scheme, say, OK, u of k plus 1 minus the right hand side of the scheme, which is u k plus delta t times f of u k, I'm going to get-- tau, for example, in Forward Euler is O of delta t squared. Which, for any scheme, or for Forward Euler, it is because k is equal to 1. Delta t squared is actually delta t to the k plus 1 power.
So why is there a plus 1 when I compute the local order of accuracy by figuring out the truncation error of the update scheme? Well, I get k if I look at the differential. If I look at the approximate, the time derivative could be [INAUDIBLE]. Yeah?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Right. We still need a manipulation. You can base the finite difference approximation of the n-th derivative and the update formula. So how do you go back and forth from this and this? How do you go back and forth from the approximation or the PDE to the update scheme? You have [INAUDIBLE] to multiply these two equations by delta t in order to get to here. You need to multiply this by delta t to get the update scheme.
In other words, the approximation error of this update scheme is equal to the approximation error of the time derivative multiplied by delta t. Therefore, of course I get one more delta t in the truncation error. So that is the reason we ask you to subtract 1 when we figure out the local order of accuracy. If you figure out the order of the truncation error, it's the single step upward. You can go from the truncation error of the single step updates to the truncation error of the time derivative. You divide this whole thing by delta t. And dividing [INAUDIBLE] by delta t, you get O delta t to the k, right?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, yeah, sorry. Thank you for that. Let me use p for the order of accuracy. It is not by purpose.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Ah, sorry. This is k. Those are k's. Right, I confused myself.
AUDIENCE: [INAUDIBLE]
QIQI WANG: All right. Yeah, so the k is the time step number, and p is the order of accuracy. Yeah, thank you.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Any other questions? So that is, I think, the one order difference in the local order of accuracy in that case is a big source of mistakes in figuring out the order of accuracy local order of accuracy analysis. So try to remember that.
OK, so one thing, in the PDE context, if we approximate a PDE by a set of ODEs and solve it using ODEs by [INAUDIBLE], what do you think is going to be the order of accuracy we'd get at the end? How much area are we making in that approximation?
All right, what determines the accuracy if I approximate a PDE by an ODE and solve the ODE using, let's say, a Forward Euler?
AUDIENCE: [INAUDIBLE]
QIQI WANG: The accuracy of the ODE scheme?
AUDIENCE: [INAUDIBLE]
QIQI WANG: And the accuracy of the discretization, actually, both. Actually, both. So if that happens even for some advanced graduate students when they look at the order of accuracy of a PDE discretization, for example, it's quite usual for them to have a plot of the area versus the grid. And often what they find out, as I refine the grid further, I no longer improve my accuracy. So what happens?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yes, because then I'm replacing my grid spacing. And if I'm not also decreasing my cell step size, if I'm not decreasing my delta t, even though my spatial discretization is extremely accurate, has no error, I feel that a truncation error would be proportionate to delta t. So if I don't decrease my delta t, I'm not going to reduce my discretization error further.
So the ODE discretization scheme makes an error. And that error is actually sometimes overlooked when you do PDE analysis. So remember, there is still going to be an ODE component, always, even if you look at PDEs, time dependent PDEs.
OK, so this is the local order of accuracy. How does that relate to the global order of accuracy? So what if I do global order of accuracy? How does-- I mean, it's a very simple answer, so please answer me. How does local order of accuracy relate to the global order of accuracy? Global is one more than local?
AUDIENCE: They're the same.
QIQI WANG: Or are they the same? If it is consistent, yes.
OK, global is equal to local if the scheme is zero stable. That's the Dahlquist Equivalence Theorem. Yes. Right, so local and global order of accuracy are the same, right? There is no plus one or minus one. The plus or minus one happens over here when you figure out the order of the one-step truncation error. Then you can get local order of accuracy, you need to subtract one from the delta t. But there is no plus or minus one in the global and local order of accuracy.
AUDIENCE: [INAUDIBLE] one or greater [INAUDIBLE]. You can't have like zero order of local accuracy?
QIQI WANG: Oh, yes, this is provided we have a consistent scheme, provided we have at least equal to 1 here. Right, yeah, thanks for that. So the whole thing is under the assumption that p is greater than or equal to 1. I mean, that's usually the case. We usually only look at vertical schemes that are consistent, right?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, the definition of consistency is p greater than or equal to 1, local. That is a consistent scheme. Yes?
AUDIENCE: And the practical way of finding that is the sum of the--
QIQI WANG: The practical way of finding what?
AUDIENCE: If p is greater than or equal to 1 [INAUDIBLE] discrete surface [INAUDIBLE] equals 0.
QIQI WANG: A practical way of finding out if a scheme is consistent is by doing the global truncation error analysis. You have to look at the Taylor series expansion to find out if the scheme is consistent or not. I don't think there is an easier way to find out if a scheme is consistent or not.
AUDIENCE: There was something in the notes about the coefficients of the non-derivative terms sum to 0. Might be a little bit [INAUDIBLE].
QIQI WANG: Yeah, OK, so there is a condition that the coefficients of the constant terms sum to 0. But I think that is a necessary but not sufficient condition.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Right, right. So if I flip a coin of [INAUDIBLE], we're going to get a sum of 0, that is not a consistent scheme.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, the sum to 0 only means that you get consistent approximation of the d by dt equal to 0 equation.
AUDIENCE: [INAUDIBLE]
QIQI WANG: OK, again, if you look at a PDE, the global order of accuracy is going to be determined both by the spatial and temporal discretization. So even though if you look at how the area decays as you refine your grid. You still need to be careful with the time derivative term. You need to refine your time step at the same time. And maybe you have to refine your time step even more if your order of accuracy in time is less than the order of accuracy in space. All right, any questions on accuracy?
And accuracy, if you need to figure out what a scheme is, the order of accuracy is a good way to check what the scheme is. Because different schemes have different orders of accuracy. Forward Euler, Backward Euler are first-order accurate. Trapezoidal rule and midpoint are second-order accurate. And there are more advance schemes. Many of them have an even higher order of accuracy. So that's a good distinguisher of different schemes.
OK, eigenvalue stability. OK, can somebody tell me? We looked at zero stability, right? That is, what makes a scheme have a global order of accuracy equal to local order of accuracy? How does eigenvalue stability, how is that different from zero stability? Yes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Hm? If the general case--
AUDIENCE: [INAUDIBLE]
QIQI WANG: Right. So eigenvalue stability, in some sense, is an expanded version of zero stability. So it basically says that my solution, Vk, is bounded for the question of du dt equal to lambda u. So this eigenvalue stability, while zero stability says that the solution is bounded for du dt equal to 0.
So when you talk about eigenvalue stability, you have to give me two things for me to determine if a scheme is eigenvalue stable or not. You have to give me the scheme. What scheme are you are using? Are you talking Forward Euler, Backward Euler, midpoint? You have to also give me something else for me to determine if I'm going to have eigenvalue stability or not. What is that?
AUDIENCE: The governing equation?
QIQI WANG: Hm?
AUDIENCE: The governing equation.
QIQI WANG: The governing equation, or more specifically, lambda times delta t. Right? You have to give me these two things in order to determine if a scheme is eigenvalue stable or not. And if you search for MIT math links, eigenvalue stability, the first thing you're going to get on Google at least is this thing. Always run, OK, run. You have to do a bunch of Java things. Come on.
Yeah, that is going to give you the eigenvalue stability. So you need two things. One thing is the scheme, right? For different schemes, you are going to get a different plot of eigenvalue stability. So if I get Backward Euler, that is my stability region. And this plot is a plot [INAUDIBLE] depending on the distance lambda delta t, right? That's the second thing you need to tell me in order to [INAUDIBLE] eigenvalue stability. You also need to tell me lambda delta t, and that determines one point in the complex set.
And depending on [INAUDIBLE] and eigenvalue stable [INAUDIBLE] or the eigenvalue unstable [INAUDIBLE]. Right? Any questions? So questions over there? OK, so yes, tell me a scheme, and tell me lambda delta t. [INAUDIBLE] I'm going to tell you I'm stable or not. And one thing is for Backward Euler, as your delta t decreases, whatever is happening, for the same lambda, as my delta t decreases, I'm more and more zooming into small regions and in all regions. And the smaller I get, the more I'm approximating the stability of the true ODE.
The value of the true ODE is the [INAUDIBLE]. If lambda has a [INAUDIBLE], then that's unstable behavior. If lambda has a [INAUDIBLE] value, then it's stable. Must be stable, right? If I zoom in [INAUDIBLE], that's going to be more and more of what I get.
And same thing for Forward Euler, right? As I zoom more and more into the [INAUDIBLE], that's what I get. And the same thing for trapezoidal. It's really depending on if I zoom in or not, but the same thing for RK2. As I zoom in [INAUDIBLE], I get [INAUDIBLE] stable, right term being unstable. And same thing for [INAUDIBLE]. I get the behavior that on the left hand side is stable, on the right hand side is unstable. Yes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: For the Backward Euler, I can easily [INAUDIBLE] delta t, and it actually is going to make an unstable equation stable, yes. That actually happens with some PDEs. Now, you can even sometimes, for example, [INAUDIBLE]. And a lot of the fluid flows are actually unstable. So if you're looking at, for example, water shedding behind the [INAUDIBLE], what you see there is an unstable flow field. But you can actually get a stable flow field if you use Backward Euler. You can use Backward Euler with a really [INAUDIBLE].
If you do that, you can actually converge, force yourself to converge to an actually unstable closed solution.
AUDIENCE: [INAUDIBLE]
QIQI WANG: The question is, are you still consistent? The question of consistency is for ODEs. If we find that lambda delta t goes to zero, consistency means the [INAUDIBLE] behavior of the ODE as lambda delta t goes to 0, in this case, that's always consistent because if your lambda delta t goes to 0, and approximating the true behavior [INAUDIBLE]. Consistency has nothing to do with if I give you a really big delta t. If I give it a really big delta t, then it has nothing to do with consistency. Because consistency is behavior of the linear [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, OK, yes. I said converges in a different sense. So the convergence we are talking about here, as I decrease my delta t, as the numerical approximates the analytical solution, what I mean over here is that when [INAUDIBLE] delta t is low, as I've [INAUDIBLE], grows closer and closer to a solution of the time independent equation.
So that's what I was saying when I spoke about convergence earlier. It has a-- I think that this is not a very good equation, but convergence in [INAUDIBLE]. One thing is when I decrease my [INAUDIBLE] time and space, I might converge in the analytical solution.
The second thing is in an iterative scheme, I guess you're going to learn more when you go through more advanced classes. When you apply an iterative scheme, trying to compute the solution to a steady state differential equation, these states are [INAUDIBLE]. The steady state never stops closing. And as I increase the iterations, do I get closer to the solution of the steady state equation? It's more like a convergence in the sense of we can apply Newton-Raphson. We can apply Newton-Raphson to solve the founding equations. That happens to be more [INAUDIBLE] of your iterative solution. Convergence moves, as I do more Newton-Raphson steps, do I converge to the solution for the [INAUDIBLE] equation?
So there are two completely different concepts of convergence. One is, as I decrease both t and delta x, I get closer to a solution. The second is as I why iterate more, do I converge to a solution?
AUDIENCE: [INAUDIBLE] What about [INAUDIBLE] So what it's saying is that [INAUDIBLE].
QIQI WANG: Right.
AUDIENCE: [INAUDIBLE]
QIQI WANG: OK, I'm--
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yes, so I'm going to draw what [INAUDIBLE]. So if you look at one where there's [INAUDIBLE] over here, that means a lambda that is that type of property, the real part-- and if you look at the analytical solution, the analytical solution here is something that oscillates sinusoidally while growing in magnitude exponentially. Right? So that's an analytical solution.
If they used Backward Euler with a small time step, so that is like when your lambda delta t is going to be 0.1, with a small delta t, you're [INAUDIBLE] lambda to somewhere close to the origin. So you may get a solution that maybe you wouldn't get exactly the right behavior because unless your delta t is infinitely small. But you are also going to get something that grows exponentially larger.
So that is when you have a small delta t. Now, if you use a large delta t, like you are scaling the lambda delta t to somewhere larger along the same line, because the delta t is real, what happens is that we get a stable solution. So although analytically the solution grows larger, you are expected to get a solution that looks more like this. So it is going to get to the wrong answer qualitatively even. And of course, this is because you're using a super large delta t. You're using a delta t that is actually much larger than 1 over the magnitude of the eigenvalue. Right? Does that make sense? Yes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yes.
AUDIENCE: [INAUDIBLE]
QIQI WANG: So you're asking a very good question. So for a system that is analytically unstable, what is a good way of telling my numerical scheme is doing a good job or not? This is a much deeper question than I can answer. Yeah, there is a [INAUDIBLE]. If the solution analytically is unstable, that means to approximate using numerical methods is extremely difficult. If you make a small error in the beginning, even if you have a small error let's say in terms of [INAUDIBLE], even though we'll assume you are doing everything exactly starting from time step 0, you're still going to get a [INAUDIBLE] error as you come to here. Just because the equation itself is unstable. The equation itself being unstable means you can make a small perturbation over here. That perturbation will [INAUDIBLE] grow larger and larger as we integrate more and more.
So treating a case like that numerically is very difficult. And if interested, let's talk more about that after class. All right, any other questions? OK. And by the way, that's something I'm actually looking at in my research right now. So very good question. I'm very impressed.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, that's why people can't predict the weather, right? I mean, they try to solve a PDE to get the weather seven days later. But you know it's going to be a not very good solution by experience. So that's exactly what they're trying to do, solve an unstable system forward in time for something like seven days.
All right. OK. I got a question of, what is the advantage, what is the main advantage and disadvantage of explicit versus implicit methods? Let's do a comparison here.
You all did the project. So you all kind of know the disadvantage of an implicit method. What is that?
AUDIENCE: [INAUDIBLE]
QIQI WANG: A lot of coding, right? Why lots of coding?
AUDIENCE: [INAUDIBLE]
QIQI WANG: You have to solve a non-linear equation. Solve non-linear equation, right, in every time step. And the way we solve it is using Newton-Raphson. So we have to apply Newton-Raphson equation within every time step. That means a nested loop within every time step. So you have an outer loop that is looping through the time step. Within the outer loop, you need to have the inner loop that does Newton-Raphson iteration. So of course, it's much more complicated than explicit schemes, right? Where you don't need to solve any non-linear equations. That's why it is explicit, right?
OK, but now, what is the advantage of implicit schemes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: It's way more accurate. But I wouldn't hold that as a rule. I mean, in the problem, yes. We used an implicit scheme that turns out to be way more accurate. But the reason may just be our implicit scheme is eigenvalue stable, right? The explicit scheme we were using was [INAUDIBLE]. That turns out to be eigenvalue stable only along the imaginary axis. So accuracy is actually not the main driver of people adopting implicit schemes over explicit schemes.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yes, the main driver is a larger stability region. Larger stability region, as we were just looking at. OK, just for example, compare Forward Euler with Backward Euler. Forward Euler, tiny, right? Backward Euler, the region where it's unstable is tiny, right? That's kind of an extreme comparison, but gets the point through. Yes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: OK, good point. What happens if Newton-Raphson doesn't converge in the implicit method? And now, when you say converge, you don't mean as delta t and delta s go to 0, right? You mean as my iteration goes to infinity, doesn't converge?
So it's a very good point. Because yes, if you have a very non-linear problem, if you use a super high delta t, it's quite possible your Newton-Raphson doesn't converge. So if you use implicit methods, there is actually an implicit restriction on delta t in which you cannot get through this eigenvalue stability analysis. It is actually a delta t that is going to enable it to converge rapidly with Newton-Raphson equation.
So one thing I think you can-- by just a little bit of analysis you can find out is that as I decrease my delta t, my Newton-Raphson is going to have a much easier time to converge. In fact, if my delta t is very small, then my Newton-Raphson, I'm going to have a very good initial guess of my Newton-Raphson. Because if my delta t is very small, then my next step solution is going to be pretty close to my current time step solution. Right?
The change of the state over the two steps wouldn't be that large. And Newton-Raphson will always converge if your initial guess is close enough. So that's the nature of Newton-Raphson because it uses a linear approximation to get the root of that linear approximation. I can talk more about that.
But if you have a close enough initial guess, Newton-Raphson will always converge. Therefore, by decreasing, if Newton-Raphson doesn't converge, a very straightforward recipe is decrease your delta t. And that is going to make the change between one state and the next time step state closer, and therefore give you a much better initial guess to Newton-Raphson. And that is going to allow you to converge much easier. So yes, Newton-Raphson actually can diverge if you have a super non-linear problem and you use a super large time step.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, yes. By the way, I can make the same argument for [INAUDIBLE]. If I get an unstable system, I can always decrease the time step so that I get into the stability region of the implicit scheme. That is still except for the time step restriction for Newton-Raphson, we said over here, is set by a different mechanism as the largest delta t I can handle up here.
The largest delta t I can use in the explicit step is done by how large, by a larger eigenvalue. It is really-- we can really obtain it from a linear analysis. And linearizing the equation [INAUDIBLE] the largest eigenvalue, I can find out what is the largest time step I can take here. So [INAUDIBLE] the implicit method is governed by if Newton-Raphson converges. And Newton-Raphson will always converge in one step if I have a linear equation. Newton-Raphson will converge in one single step if I have a linear equation.
So the time step restriction here is set by the convergence of Newton-Raphson, not by the eigenvalue step, right? How non-linear the equation is. So there are some integration methods where people actually separate out the linear part that has super large eigenvalues and the non-linear part.
So for a linear part that has super large eigenvalues, they do it implicitly. And for the non-linear part, they do it explicitly. And there's a class of methods for the [INAUDIBLE]. Fancy name, but short for implicit, explicit methods that creates different parts of the equation, either explicitly or implicitly. And as you can guess, the part they treat implicitly is going to be the linear part that has super large eigenvalues. So that you converge in one step, you can do Newton-Raphson, and you avoid the time step restriction set by the large eigenvalues.
All right, OK. Right. So this is really especially in stiff problems. And what is a stiff problem? A stiff problem is actually defined in terms of if you use an explicit integration scheme.
If you find yourself having to use a super small delta t not because you want super high accuracy but because it's still going unstable if you use a larger delta t, then you know you have a stiff problem. That is really the easiest way, I find, to define a stiff problem. That is like, you are forced to use a small delta t by the stability region, not for accuracy reasons, not for wanting higher accuracy, but simply trying to avoid [INAUDIBLE].
OK, so then you get a stiff problem. Any questions? Yeah, so any question, please raise it. Otherwise, I'm going to review Newton-Raphson a little bit, and there is no more things I'm planning to cover. So I'm kind of expecting questions from you guys.
AUDIENCE: Finite difference [INAUDIBLE]
QIQI WANG: Yeah, finite difference, and finite volumes I include in the scope of the exam. As you saw from the previous year's questions I posted, sometimes we include it in the actual exam, sometimes we don't. So it is included in the scope, but because we just did it, I think it's a good idea to review some of the earlier stuff that you may have forgotten. And another benefit is that we are also going over this now that we have [INAUDIBLE] to give you a sense that all the stuff we learned applies to PDEs because what we did in PDEs is just to approximate the PDE using a big ODE. Right?
And you can apply implicit methods on the PDE also, except for the Jacobian. You get it's going to be something close to the matrix form of finite difference, we did in the finite difference vectors. Yeah?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Oh, OK, the example format is that it's closed notes, right? It's closed everything in the period from you get the exam to when you come to our office. So it's closed everything, closed computer, closed cell phone. Closed Wikipedia.
We are actually letting you work out, really think about the problem during the time you are looking at the problem. But even though you don't get the answer, or during our face time, we are going to also interact with you when we ask you questions that you may not have expected, or we may actually help you going through some of these. So that's kind of how [INAUDIBLE] goes. Do you have something else you want to ask?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, right. You can write things down and bring whatever you have written into our office. And in the office, you're expected to use a whiteboard or blackboard and explain to us like you're the professor, the other students, explain to us what you got. Any other questions? No?
The last thing I want to go through is Newton-Raphson method. And it's another sort of confusing-- as I said, Newton-Raphson method is a method that simply solves regular non-linear equations. And a non-linear equation can appear in 16.90. It can appear anywhere else in your future career.
So what you're learning about solving non-linear equations really goes very far, even if you don't deal with numerical methods later on. So as I said, it's a method of solving non-linear equations if the set of non-linear equations is relatively small. Say you have two equations or three equations, you can go to Matlab and use F solve to get it solved like brute force.
But if you have a large set of differential equations you have to solve simultaneously, like what we're going to have if you use an implicit time integration method, apply it to a discretized PDE, then you're going to get at least 100 ODEs you have to-- which, if you have to solve using an implicit scheme, you're going to get hundreds, maybe thousands, maybe millions, or maybe trillions of algebraic equations you need to solve.
OK, so imagine you have to solve f of u equal to 0 where u, instead of writing down a simplified version of what you would get in an implicit scheme, your implicit scheme you would get u minus u k, which is the stuff you already know, over delta t is equal to some right hand side of u, and uk, and maybe something else, right? And instead of writing the same f here, I'm just going to write a big F of u equal to 0. So that big F in this case is really defined as u minus u k over delta t minus f of u, u k, et cetera.
So this is what happens if we solve, let's say, a discretized PDE [INAUDIBLE]. U is the next time-step solution we want to get. And u k is the previous time step solution you already know. Now, if F is non-linear, you have to find the root of its big F, which is essentially the left hand side minus right hand side of this implicit update [INAUDIBLE].
Right? You need to find a u which is a vector of the solution that makes this f 0 for all the components of f. And f has the same dimension as u. If u is a 3D time step, f is going to be a 3D time step.
And you need to find these 3D numbers that make the [INAUDIBLE] f equal to 0 simultaneously. That's not an easy task. And I think F solve is going to have a hard time dealing with this if u is a high dimensional vector.
Now, what does Newton-Raphson do? Newton-Raphson starts with an initial guess. So u, for example, I'm going to use parentheses to denote the information [INAUDIBLE]. Parentheses 0 is the initial guess. I'm going to set it to u k. I'm going to set it to the solution at the previous time step.
Then, I'm going to approximate this non-linear function using a linear function. Can I approximate F of u by something linear? Let's do it one by one. Let's approximate the first component of F. Think of F having components. We're just going to approximate the first component and have all the other components follow.
OK, so if you just approximate the first component, I'm going to use Taylor series. And I'm going to use a Taylor series of a function of the [INAUDIBLE] variables, all right? So a function-- the periodicity of a function of a [INAUDIBLE] variable is pretty complicated. But the lucky thing is that I don't need to keep all the terms. I'm going to throw away anything that involves more than the first derivative.
In other words, I'm only going to keep the zeroth order term, which involves no derivatives, and the first order term, which involves only the first order derivatives. OK, now what is the zeroth order term?
The zeroth order term is F1 at my initial guess, right? What is the first order term? Actually, in this case, because it's a multivariate function, I have more than one first order terms. If I have as many first order terms as there are many u's, I'm going to be summing from i goes from 1 to the dimension of the system, let me call big N, partial F1, partial u of the ith dimension times ui minus u0i. Right? Does that make sense? That's the Taylor series expansion of F1.
We are only keeping the zeroth order terms as opposed to all the terms. If I start to write all the other order derivatives, it would get too complicated. But I'm not going to write them, so I'm going to truncate them. And I'll approximate [INAUDIBLE] method, it's a linear function. It's a linear function because this F1 of u is probably a linear function of u. But now, I'm truncating. I'm removing everything that happens over here, so if it's quadratic, or cubic, or anything. I'm only keeping the constant part that is independent of u, and this part, which is 0.
Now, I get a linear approximation of F. We all know how to find the root of a linear function? Even though it's a million dimensional? Right, that's the reason we have linear algebra, right? That's why we have matrices. That's why Matlab is called Matlab. That is because we can write functions like this, we can write linear functions in matrix form. And to solve linear equations like this, even ginormous ones, we just use linear algebra, right?
OK, so we are also approximating all the components of F, all the way to the n-th component using Fn u0, plus summation i goes from 1 to N, partial FN partial ui, ui minus ui0. So there are big N equations. We have big N of these variables. How to solve them? How to solve these coupled, N-coupled linear equations? Yes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, it's just write it into a matrix form, right? First, write it into a matrix form. So I want to try to say, OK, if I want to set all these things to 0, I'm trying to set like 0, 0, 0, 0, is equal to F1 at u0, et cetera, to FN of u0. This is this term. And how about this term? This term is just a joint matrix vector multiplication. The matrix is other derivatives, which is called the Jacobin when I put that into a matrix form.
The first row is my partial F1, partial blah. The first column is my partial blah, partial u1. The last column is partial F blah, partial uN. So each row corresponds to one [INAUDIBLE] the residual. Each column corresponds to one thing in the independent variables, in the u's. And multiplying that by order u1 minus u1 0, et cetera, to uN minus uN 0. Do you all see these linear equations are the same as this matrix?
All right? So what I did is the Taylor series expansion for all these non-linear equations and write them all into a matrix form. And what I'm trying to do is I'm trying to solve for these u1 to uN. How do I do that? I just invert the whole Jacobin. If I call this as J, then my u1, et cetera, uN is going to be equal to u1 0, uN 0 minus J inverse times my F1, FN evaluated at the initial guess.
Right? So this is the solution [INAUDIBLE]. This is one step in the Newton-Raphson, right? In Newton-Raphson, we compute the residual [INAUDIBLE] at the initial guess. And we do the Jacobian. [INAUDIBLE] the Jacobin inverse times the residuals.
Then, you stop at the result from the [INAUDIBLE]. And get your next step solution. And what is the next step solution? The next step solution is the zero of the linear approximation? All right? Which hopefully gives you a better initial guess than u0. And what you do next is you call this set of u's to be u parentheses 1. We call these as u parentheses 1. And then, you linearize again on these.
So the one dimensional analog is this. You try to find F of u. OK, you want to find F of u. You start over an initial guess. This is my u0. I construct my first order Taylor series approximation, which is what in this case? If the zeroth order term is going to be this, this is the zeroth order term. It's a constant, right? It's just F of u0.
The first order term is a derivative at this point, which is going to give you this. So this is my first order term is going to give me the tangent line. My u1 is the solution of the Taylor approximation equal to 0. So this is my u1. And then, I'm going to linearize around this again. I'm going to construct another first order approximation, and go to this point, et cetera, et cetera, until I converge to the zero at the blue line?
Now, it's hard to draw this picture when both u and F are million dimensional. But the concept is the same. I construct a linear approximation. And the linear approximation can be arranged into matrix form using the Jacobin. Jacobin is just a matrix containing all the derivatives of Fu's. And I keep iterating. I keep iterating. I solve the 0 of the linear approximation, linearize again at that new point, and then solve, get the zero of that linear approximation. I linearize again at the new point, I get the zero of that linear approximation, and then linearize again at the new point. Any questions on this?
General technique of solving large systems of non-linear equations. Using the fact that we are already good as solving large systems of linear equations. OK, now I'm all out of things. Just time to answer your questions. Yes?
AUDIENCE: Will we be expected to [INAUDIBLE]
QIQI WANG: No, you don't need to quote anything. You don't need to have your hands on the keyboard.
Yes?
AUDIENCE: [INAUDIBLE]
QIQI WANG: What does big F represent in this case? Big F represents the function of u if we want to get the zero. If I want to solve [INAUDIBLE], if I use an implicit scheme and I want to evolve from this case to the next u, to the next time step, then my F of u would be the left hand side of the scheme minus the right hand side of the scheme, right?
Imagine if I have Backward Euler. Then my big F would be u minus uk divided by delta t minus F of u. That is the Backward Euler. If I trapezoidal, my big F would be u minus uk over delta t minus half of F of u minus half of F of uk. All right? And the difference is implicit schemes are going to get a different value. I have a Backward Euler, [INAUDIBLE], half of that F of u plus half of that uk. And that is in the derivative.
Whatever the scheme is, F is the thing along with [INAUDIBLE] zero that big F is [INAUDIBLE].
AUDIENCE: [INAUDIBLE]
QIQI WANG: It's a [INAUDIBLE], yes.
Questions? After this, we are going to [INAUDIBLE] the midterm.
AUDIENCE: [INAUDIBLE]
QIQI WANG: All right. Anything else I can talk to you about? All these materials? Oh, and just to remind, we are really looking at the measurable outcomes. And don't forget to go [INAUDIBLE] something that [INAUDIBLE]. Analytical solutions [INAUDIBLE] what is the solution of u [INAUDIBLE]? And things like that. And if the solution is stable, when it is unstable, you need to know that and things like that. So this is important. Yeah?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Should d of dt be equal to lambda u? You need to solve it analytically. Huh?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah. For d of dt over the lambda u, you need to be able to solve it analytically, yeah.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Huh?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yes, yes. D of dt equal to au, yes, right. Any matrix, A, right? You need to be able to solve it.
AUDIENCE: [INAUDIBLE]
[LAUGHTER]
QIQI WANG: Yeah. My A can be a trillion by a trillion. And you need to know how to solve it. You are not expected to--
AUDIENCE: [INAUDIBLE]
QIQI WANG: Hm?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah, of course, you won't be able to do the eigenvalue problem of A in your head if A is a 1,000 by 1,000 matrix. But you need to know how to get the analytical solution of the ODE.
AUDIENCE: [INAUDIBLE]
QIQI WANG: Hm?
AUDIENCE: [INAUDIBLE]
QIQI WANG: Yeah.
AUDIENCE: [INAUDIBLE]
QIQI WANG: All right. Anything else? And homework here. If you didn't [INAUDIBLE] homework for [INAUDIBLE], the previous homeworks, you can always come to my office and figure it out.
