Chapter 8.2 - Scientific Computing: Time-Stepping Errors and Stability

0:08 We're going to continue talking about time stepping schemes for differential equations. And specifically now we're going to get into some of the the details of both the error analysis as well as the stability analysis of the scheme. And what I'm going to show you are sort of things that people often look at in terms of being able to understand the scheme that they're going to apply for modeling their system. Obviously, what we want is the most accurate scheme that is stable and of course we want it to be fast as well.

0:41 But the main thing you have to ensure upfront is that you have a stable and accurate scheme. The speed is something we often work on. And with today's modern computation computational architectures, you know, we've really improved the speed on on our processors. So often people think less about that and they think more about making sure you have accuracy, making sure you have stability. So let's go after this and start thinking about it. It's all going to come down again to tailaylor series expansions. And those tailaylor series of expansions are going to allow us the theoretical framework to be able to make statements about stability and about accuracy.

1:19 Okay. So again this falls broadly under scientific computing, numerical differentiation, integration, differential equations and also what follows after this would be partial differential equations. All these are sort of these traditional scientific computing schemes which allow allow us to solve science-based problems which are typically based upon figuring out how relations of the derivatives space and time relate to each other to produce some dynamics that might we might observe in practice.

1:50 All right, so let's get back to it. We're focusing on differential equations. That's the generic form. dydt equals f. So f is like your governing equation. It specifies the evolution of the system. This could be for instance F= ma. Force equals mass times acceleration. so all of our many of our standard models fit into this form. And many more fit into generalizations of this which are partial differential equations quantities that don't just depend on time but they depend on time and space. In either case, this is what we're trying to solve. Our basic structure in solving it as we showed last time was say well I can do a tailaylor expansion and now we're going to get to the error act error analysis.

2:33 So I want to come back to this error analysis because we didn't dive into it too much last time. We more painted a picture of how you would construct schemes but now we're going to dive into how you would do a formal error analysis. So let's start with a tailaylor series expansion. Again yft plus delta t is y oft plus delta t. Here this is an exact formula and I truncated at the second derivative with some c there which is between t and t plus delta t. So this is an exact formula. So tailaylor series expansion with the term and what I'm actually going to do with this this is going to be my truncation term. I'm going to throw it away but notice it has a delta t ^ squ in front and that's going to sort of start dictating the size of the error.

3:16 So here we go. So the easiest thing to do is you can manipulate this equation into the following which is the oiler scheme. This is what we did in the very first thing in the last lecture which is my solution in the future y of n plus one. So in other words at t of n plus one is right here is equal to the current solution y ofn plus delta t plus f and then notice here this is my error term order delta t squ that comes from that truncation right there. So this is the oiler scheme. Now in the oiler scheme before we just simply ignored this term and we just said okay here's a stepping scheme but now we're going to bring it back and start to look at the consequences of having that delta t^2 term around.

4:01 So the error analysis begins with the following you're saying look my error itself right and by the way you know I'm approximating the error is going to be what y of tk is minus the actual value. So this is my estimate of what the future is versus what it actually is. Okay? So I'm going to look at that difference. So remember my time stepping scheme is going to walk me to the future. Okay? And I'm going to then look at the difference between that and what the actual solution is in the future.

4:35 Okay? So that's what I'm going to estimate here. So let's talk about what this error is. So actually y of tk + one this is the actual solution at t of k + one. So I actually say okay this is the truth. This is the ground truth. This is my estimate of it. So I start to produce the future state. I take my current value at this time and add a little bit of something. So this is my time stepping scheme. This is how I approximate the future state which is supposed to be this. So I'm looking at my approximation to the future versus the actual future and let's call that epsilon of k + one. Okay, from the time step t of k to t of k + one.

5:19 So that is my error. And so now I'm going to actually break up the error into two pieces. I know that every time I take a step I'm going to make an error. And that's called the local error. The local error is just basically how much of a mistake did I do in taking one step into the future. Now I actually know exactly what this is from the tailaylor series point of view which is this is what I dropped. I dropped this correction term which had this delta t squared in front. Okay. So in other words when I take a step with oiler the oiler term I drop this piece of the tailaylor expansion and so this is of size delta t ^ squ and that is then the error I make in one step.

6:06 But I don't really care about one step. I want to remember solve this differential equation, walk it out into the future. And I would like to know is well, how much error do I accumulate into the future? Because if each step makes a mistake, that mistake is made upon another mistake and so forth. So what's the cumulative error? So that's the global error. Let's call it a capital E of K. And the capital this E of K is the sum of all the errors you're making. Now remember every step I make I'm making a an error of delta t ^ squ. So I just sum all these up where I don't know what the cfjs are. But this is in some sense I replace all of this there's k steps by some average value of the second derivative times k steps. So this sum each one has some value. I replace it by its average. There's k steps. I take the average value average times k. That is my approximation for the error in the future. Okay.

7:02 But then I can also ask the questions. Okay. Well, so I have this delta t ^2 over 2. I have the second derivative which is the average value. I have k steps. But k steps are equal to well if the interval goes from a to b and I have steps of delta t. So delta t * k is equal to b minus a. So I'm going to take this k and replace it by b minus a over delta t. In other words, this here is a number of steps related back to delta t.

7:37 Now notice what this does. It means there's a delta t on the bottom here which cancels out one of the delta t's here. So I had a delta t ^ squ here, a delta t there. When I do this math out, I get delta t here. which means that this whole error term is order delta t. So the local error one step is delta t ^ squ but the global error is delta t. Okay. So when I walk this out into the future my error is deltat t which is quite large for the oiler scheme and typically why nobody uses an oiler scheme. It's not even mentioning the stability properties of boiler, right?

8:20 So if delta t is 10 minus 2, my air is 10 minus 2. What you really want is something much more accurate. I mentioned previously in the last lecture fourth order runa which allows me to take steps of delta t and it gives me an accuracy of delta t to the 4th. So if I take steps of minus 10 minus 2, my accuracy is 10 minus 4. So this is actually not a very good scheme but we can work everything out with it in terms of the local and global error. And like I said mostly you're concerned with is the global error.

8:53 Okay. But it's not just error due to tailaylor series expansions which I just showed you. One of the things we swept under the rug and which something that we brought up when we were doing numerical differentiation is this idea of numerical roundoff which is when you take a step you might be just tempted to say take smaller and smaller of delta t's but just like we did in numerical differentiation that's actually not a viable strategy once you take delta t small enough so let's talk about this if I do a slope approximation which is rise over run and I have this error here the error error is not just the truncation error but now the error includes numerical roundoff and we're going to see what the impact of numerical roundoff is okay so the impact of numerical roundoff so first of all whenever I evaluate the y of n plus one it's actually numerically evaluated some approximation right so remember finite precision in your computer is 10 - 16 precision for let's say double double precision so it gives you some value of n plus one plus the roundoff term here. So when I think about marching forward when I take this derivative estimate I'm actually estimating its numerical value with roundoff and the error now includes both numerical roundoff and truncation error.

10:17 So that's there's two two things compounding my error estimate. So the total L is right there. Roundoff truncation. Here's the roundoff error. There's the truncation error. So this is looking a lot like what we did when we did numerical differentiation, which is you get some formula like this. And you can start to see now numerical roundoff is going to have a role because it's being divided by delta t. So as delta t gets really small, this term is going to get really big past a certain point.

10:50 So let's put error bounds on this. So the error bounds are pretty simple. We can just say that roundoff is bounded by some epsilon r. Typically this is 10 - 15 10 - 16. That's the biggest it can be. We'll also say that second d has a max size m. So this is in some sense looking at what's the biggest error I could make and those are the estimates for the biggest error. And so I can put them back into the error formula for the absolute value. And here's what I get.

11:18 So this thing here is the approximation that comes from the numerical roundoff from the truncation of the derivative which gives you this formula here. And notice again the divide by zero is going to happen if you go to delta t goes to zero. So past a certain point numerical roundoff is going to make your error go bigger. So you can't just simply take smaller and smaller and smaller delta t's because past a certain point the error is going to grow. Okay.

11:47 So how do we find this minimum delta t I should be taking? You take the derivative of this with respect to delta t. Here it is. And what you find is the optimal step size is given right here. Okay. So simple calculus. This is a minimizer. This delta you find the minimum and delta gives you the smallest step size. Now let's go back to e of r being let's say 10 - 15 10 - 16. This would give you something like a step size of 10 - 8. So you just shouldn't take step sizes below 10 minus 8 because your error will actually increase. Okay, not because of the truncation but because of the roundoff. So this is important information for you to keep hold of which is yes you can take smaller and smaller delta t's and often this is what's done in practice but make sure you understand there's a limit of how far you can go before it actually starts to backfire on you as a strategy in some sense for gaining accuracy.

12:42 Okay, so that's one issue. This is the accuracy issue sort of detailed in these calculations. Right? So the accuracy issue is I want to go ahead and take smaller steps. What's my error that I incur? the nice thing about a scheme like a runa cutter fourth order which I advocated for the last lecture is that you know the error the global error goes like delta t to the fourth. So it's a fantastic scheme because it really pushes the error down. Obviously, it's a more sophisticated scheme than oiler, but you get so much benefit both in terms of accuracy, but also you get benefit in terms of stability. The stability analysis is very difficult to do for runata, but we're going to do stability analysis just for the oiler to show you how it works.

13:25 Okay, so let's take a simple OD to illustrate this point. Here's a OD. It's a first order OD. dydt equals lambda y and y not is this is the initial condition at time zero it's y z so what you're seeing here this is a a simple od I can actually write down the solution very easily for this okay but what I want to do is now think about how do I if I were to numerically approximate it with oiler what would happen okay so here's what the oiler method would say the solution in the future is what it is now plus delta t time lambda times the function f of y. The function f of y is lambda y. There it is. So I can put that in here explicitly and then I compute this and here's what I get for this 1 + lambda y lambda delta t is y ofn. So here's my stepping scheme and in fact if I take n steps into the future it's just this is this is exactly what the oiler will give me for that differential equation. Okay, so here is the term that's multiplying to the nth power the initial condition as I walk forward with an oiler scheme.

14:36 So the question comes of course well what does the stability of scheme mean? Stability of the scheme means that I'm going to actually look at this from the point of view of saying all right I'm going to go step into the future but now suppose my initial condition isn't just y not but there's a numerical roundoff. So this is what the epsilon represents. So I have a little bit of numerical roundoff just on the initial condition.

15:01 I'm just going to ask what happens to that little bit of numerical roundoff as I march out to the future. Okay. Well that little bit of roundoff is going to hit this term right there. raised to the nth power. So what's that going to do? Well, your error then just from the roundoff is going to look like this. It's going to be e 1 plus lambda delta t to the nth power times that roundoff. And of course we want to watch this far out into the future. So n's going to get large. So you have to ask the question in the end large limit.

15:33 You can already see there were going to be some problems maybe with this with this error as we walk out to the future. poop. All right. Missed a slide here. Okay. All right. So, here are the error options. Option one, option two. What is the absolute value of 1 plus lambda delta t? If it's less than one, then the error actually goes to zero. So, even if you have a numerical roundoff, this numerical roundoff in this system, there's no issue. the numerical roundoff goes to zero. However, if 1 plus lambda delt is bigger than one, then the roundoff is going to go to infinity.

16:19 So this is a problem because now you're walking forward and just numerical roundoff itself is going to blow up to infinity which is a problematic piece. So we normally think about this as that's an unstable scheme. If the error goes to infinity, right? I have this little tiny error and it goes up to infinity. It's unstable scheme. Okay. And that's the way we illustrate it in this sort of mapping. Remember that all a time stepper is is a is a map between t to t plus delta t. And you march it far into the future. So it's an iteration scheme. And if that iteration scheme blows up, your numerical method is unstable.

16:55 And this is what it's telling you. Let's contrast that oiler method which was unstable which can be unstable to what's called backward oiler. And I'm going to contrast it because we're going to make this very simple change to the stepping scheme and it's going to have a profound impact on stability. So my future is what it is now plus delta plus delta t * f. But now notice what I'm do with the f. Instead of saying lambda and f is evaluated at the current time, it's going to be evaluated at the future time. Often times this is called an implicit scheme because you actually can't evaluate the future very easily because you don't have it. So you have to solve a set of equations. So it can be costly. However, people will often do this for one specific reason because implicit schemes tend to be much more stable than explicit schemes.

17:48 So let me show it to you. If we do the same kind of analysis as before, what we can find is the future iterations are given by now by this formula 1 over 1 minus del lambda delta t to the n power y0. So you just get a sl slightly different formula here than you had previously. All because of the change of this guy right there. So that's the idea. You walk it to the future and this is telling you how that is going to go into the future for you. All right.

18:17 So the again we can ask the same kind of questions. What are the implications? Well, here's what the error is doing now. So if I have a little bit of numerical roundoff on that y0 that numerical roundoff is going to be multiplied by this factor here raised to the nth power. And then I have two things happening now. I'm going to compare these two schemes. I have this scheme here for backward oiler and I had this scheme here for forward oiler. So if I call lambda deltat t the forward oiler has this condition for stability.

18:48 The backward oiler has this condition for stability. Okay. So 1 plus z is greater than one. one over one minus E greater than one. And if I plot these, you'll see that it has a huge impact on the stability of these two algorithms. So let me plot the picture here. So this is the forer oiler scheme. The forward oiler scheme is only stable in this little ball right here. Everywhere else you step, your scheme will blow up. Okay? So you have to be aware that the only a little region which allows your scheme to actually be stable.

19:26 In contrast, backward oiler is stable everywhere except for that little ball right there. So as long as you stay out of that ball, your scheme is stable. And so these are the kind of things that you have to start thinking about. And this is by the way why people typically if people are going to commit to implicit scheme, it typically tends to be because of this. It's all about the stability. I take steps and I don't have to worry about the solution blowing up because I've got these big range of areas where I can step and it's going to be completely stable. Whereas if you have some of these exploit schemes like rungakuta sorry like oiler then it's problematic because you have a very small area where you can step otherwise any numerical roundoff is just goes to infinity and then your solution is really no good anymore. Okay. And so stability is balanced against accuracy.

20:21 And so this kind of analysis that we do here is actually very difficult to do with most schemes, right? Generating these pictures about where you can stably step or not stably step based upon the differential equation, especially nonlinear differential equations is very hard to evaluate. But it is something that sits there underlying all these schemes is this idea of stability and having to take care of stability. So that's one thing we want to think about and in in all these schemes and I was just going to emphasize that when we start thinking about again implementing a scheme stability and accuracy rung of fourth order is should be your go-to because it has great accuracy properties great stability properties retains good speed properties as well. So it's kind of a balanced profile of a stepper and that's why people use it so commonly. and so you always have to think about these schemes and again why do people use it? Because of this kind of analysis I just showed you. Everybody's trying to balance accuracy, stability, speed. So you want is an algorithm that can do well on all three of the criteria and you know the the field itself has developed so many techniques to be exploited in different application regimes where sometimes for a specific application regime they have really great steppers that are sort of in sense bespoke for that community because it takes advantage of some properties which allow them to solve code faster. Okay. and retains accuracy and stability at the same time. So those are things to be thinking about and I think what we want to do next is actually implement this in practice. But these are things you think about in your design for a stepping scheme.