1 00:00:00,140 --> 00:00:02,560 Welcome to chapter 12 of the course 2 00:00:02,715 --> 00:00:05,495 Numerical Methods for Computational Science and Engineering. 3 00:00:05,650 --> 00:00:09,990 This chapter is still concerned with numerical integration. 4 00:00:10,145 --> 00:00:13,185 As you know, this means the approximate numerical solution 5 00:00:13,340 --> 00:00:16,420 of initial value problems for ordinary differential equations. 6 00:00:16,575 --> 00:00:20,175 Well, you may wonder why we need another chapter. 7 00:00:20,330 --> 00:00:25,050 In the previous chapter, we learned about explicit single step methods, 8 00:00:25,205 --> 00:00:28,845 in particular, the powerful family of Runge Kutta single step methods. 9 00:00:29,000 --> 00:00:33,560 And this family provides methods of arbitrary order. 10 00:00:33,715 --> 00:00:39,375 And it's even amenable to a smart local in time adaptive time stepping strategy. 11 00:00:39,530 --> 00:00:41,150 Why do we need anything else? 12 00:00:41,305 --> 00:00:45,825 Well, there are initial value problems for which the methods that we know so 13 00:00:45,980 --> 00:00:48,140 far perform really poorly. 14 00:00:48,295 --> 00:00:53,935 And the most widespread, the most notorious of these families 15 00:00:54,090 --> 00:00:56,930 are the so called stiff initial value problems, 16 00:00:57,085 --> 00:01:00,245 about which you will learn in this chapter. 17 00:01:00,400 --> 00:01:05,580 So for stiff initial value problems, we essentially have to start from scratch. 18 00:01:05,735 --> 00:01:11,835 We have to figure out new families of single step methods that do a good job 19 00:01:11,990 --> 00:01:16,350 about solving them accurately and efficiently. 20 00:01:16,505 --> 00:01:18,945 This is what the current chapter is about. 21 00:01:19,100 --> 00:01:23,460 And now let's first see how explicit single step methods fail 22 00:01:23,615 --> 00:01:27,855 and analyze the reasons for this failure. 23 00:01:28,010 --> 00:01:33,850 So let's look at a simple initial value problem for the logistics scalar ODE 24 00:01:34,005 --> 00:01:39,645 and apply our embedded Runge Kutta single step method with adaptive time stepping 25 00:01:39,800 --> 00:01:45,920 that we constructed in the previous chapter and see how it performs. 26 00:01:46,075 --> 00:01:49,635 This plot gives you an overview of what we observe. 27 00:01:49,790 --> 00:01:52,870 Well, this green dashed line presents the solution. 28 00:01:53,025 --> 00:01:55,305 And in the beginning, everything looks OK. 29 00:01:55,460 --> 00:01:58,980 The method very well resolves the steep part of the solution. 30 00:01:59,135 --> 00:02:02,535 But then the exact solution becomes essentially constant. 31 00:02:02,690 --> 00:02:07,490 But what you see is that in this part, the adaptive Runge Kutta single step 32 00:02:07,645 --> 00:02:09,165 method does a very poor job. 33 00:02:09,320 --> 00:02:13,360 The approximations are, well, rather far off 34 00:02:13,515 --> 00:02:16,715 despite the use of really small time steps. 35 00:02:16,870 --> 00:02:21,030 What we see here is an outright failure of time step control. 36 00:02:21,185 --> 00:02:25,785 And you see this even more blatantly on the right. 37 00:02:25,940 --> 00:02:29,500 The red line is the time step size, dashed green line, the exact solution. 38 00:02:29,655 --> 00:02:30,935 The exact solution is constant. 39 00:02:31,090 --> 00:02:34,410 Actually, a huge time step should be OK here. 40 00:02:34,565 --> 00:02:38,445 However, the time step here lingers around 10 to the minus 3, 41 00:02:38,600 --> 00:02:40,960 which is really unreasonably small. 42 00:02:41,115 --> 00:02:42,355 So what's going on? 43 00:02:42,510 --> 00:02:47,510 Let's try to understand this and switch to an even simpler example. 44 00:02:47,665 --> 00:02:54,225 Now, this is the simple scalar linear, well, growth or decay ODE here, 45 00:02:54,380 --> 00:02:56,940 decay because lambda is minus 80. 46 00:02:57,095 --> 00:02:58,495 That's a large negative value. 47 00:02:58,650 --> 00:03:05,090 Again, we apply our standard workhorse adaptive explicit Runge Kutta time 48 00:03:05,245 --> 00:03:07,365 step is ODE 4.5. 49 00:03:07,520 --> 00:03:10,600 And well, you make about the same observation 50 00:03:10,755 --> 00:03:12,195 as in the previous example. 51 00:03:12,350 --> 00:03:16,030 The exact solution after a short time is essentially constant. 52 00:03:16,185 --> 00:03:17,465 It's even 0. 53 00:03:17,620 --> 00:03:22,060 And nevertheless, the adaptive method does a very bad job. 54 00:03:22,215 --> 00:03:25,575 Time step size stays small, unreasonably small, 55 00:03:25,730 --> 00:03:28,090 considering that the solution is essentially 0. 56 00:03:28,245 --> 00:03:32,365 So same observation as before, time step control fails. 57 00:03:32,520 --> 00:03:37,840 It does not realize apparently that the solution is extremely simple. 58 00:03:37,995 --> 00:03:38,995 Here's even 0. 59 00:03:39,150 --> 00:03:45,070 Well, maybe this has to do with this complicated explicit integrator ODE 4.5. 60 00:03:45,225 --> 00:03:50,705 So let's apply a much simpler time stepping scheme to the same example. 61 00:03:50,860 --> 00:03:55,300 So we again consider this very simple scalar linear decay ODE, 62 00:03:55,455 --> 00:04:00,375 decay rate lambda fairly large and apply explicit Euler, 63 00:04:00,530 --> 00:04:03,410 the simplest explicit Runge Kutta method we know. 64 00:04:03,565 --> 00:04:06,325 And then we make the following observation. 65 00:04:06,480 --> 00:04:11,560 So what I've plotted here is the error at final time for various values of lambda. 66 00:04:11,715 --> 00:04:13,035 And look at the scale. 67 00:04:13,190 --> 00:04:14,790 The scale is 10 to the 20 here. 68 00:04:14,945 --> 00:04:20,825 So when lambda is not too large in modulus, this is the, say, the blue line. 69 00:04:20,980 --> 00:04:25,940 Oh, well, you observe something that could be called convergence. 70 00:04:26,095 --> 00:04:29,895 But when lambda becomes larger in modulus, though negative, 71 00:04:30,050 --> 00:04:32,970 then here in this region, nothing is okay. 72 00:04:33,125 --> 00:04:36,565 We observe a blow up, at least for large time steps. 73 00:04:36,720 --> 00:04:40,680 Well, then finally, when the time step becomes small, 74 00:04:40,835 --> 00:04:44,755 then the error drops to very small values close to 0. 75 00:04:44,910 --> 00:04:47,990 So here this is okay, but here we have blow up. 76 00:04:48,145 --> 00:04:50,945 And there's a geometric explanation for the blow up. 77 00:04:51,100 --> 00:04:56,340 We know that Euler has a geometric interpretation as following the tangent to the curve. 78 00:04:56,495 --> 00:05:03,055 The blue arrows in this plot give you sort of the tangent sizes and directions. 79 00:05:03,210 --> 00:05:07,890 And what we see is that for this decay ODE, the Euler method goes too far. 80 00:05:08,045 --> 00:05:10,125 So in the beginning, it drifts off. 81 00:05:10,280 --> 00:05:15,440 And then at the next state, the tangent is already very steep. 82 00:05:15,595 --> 00:05:18,595 And the tangent vector is fairly long. 83 00:05:18,750 --> 00:05:22,670 So it goes in the other direction, but also again too far. 84 00:05:22,825 --> 00:05:25,345 And there's an oscillatory divergence. 85 00:05:25,500 --> 00:05:26,940 You could call this overshooting. 86 00:05:27,095 --> 00:05:32,175 And this is the geometric explanation for this, well, exponential blow up 87 00:05:32,330 --> 00:05:36,970 that you observe when applying the explicit Euler method to this decay ODE. 88 00:05:37,125 --> 00:05:39,245 Well, this is just a single example. 89 00:05:39,400 --> 00:05:40,760 Let's look at another one. 90 00:05:40,915 --> 00:05:48,355 So the next example is the logistic ODE that we have already studied earlier. 91 00:05:48,510 --> 00:05:51,710 But now we apply the explicit Euler method again. 92 00:05:51,865 --> 00:05:55,905 So simplest explicit time stepping scheme that we know. 93 00:05:56,060 --> 00:05:57,900 Now, guess what will happen? 94 00:05:58,055 --> 00:06:00,055 Same story as before. 95 00:06:00,210 --> 00:06:01,650 Again, here, blow up. 96 00:06:01,805 --> 00:06:04,125 Well, it's an exponential blow up. 97 00:06:04,280 --> 00:06:06,440 And there's nothing to be added. 98 00:06:06,595 --> 00:06:12,355 Once the time step becomes small, then the error suddenly drops to reasonable levels. 99 00:06:12,510 --> 00:06:13,470 It becomes very small. 100 00:06:13,625 --> 00:06:15,785 Same phenomenon, same explanation. 101 00:06:15,940 --> 00:06:21,700 Well, when you would plot the direction field associated with this logistic ODE, 102 00:06:21,855 --> 00:06:27,535 then again you see that Euler is just carried far off the true solution. 103 00:06:27,690 --> 00:06:34,010 And then it goes into exponential oscillatory mode. 104 00:06:34,165 --> 00:06:37,685 Again, due to overshooting, same phenomenon here. 105 00:06:37,840 --> 00:06:42,640 Well, explicit Euler is simple enough so that we can make an attempt 106 00:06:42,795 --> 00:06:45,355 to understand this rigorously. 107 00:06:45,510 --> 00:06:46,470 Let's try this. 108 00:06:46,625 --> 00:06:53,585 And the approach that will turn out to be successful is sort of linear model problem analysis 109 00:06:53,740 --> 00:06:56,220 based on what's called diagonalization. 110 00:06:56,375 --> 00:06:58,535 So let's see what happens here. 111 00:06:58,690 --> 00:07:05,410 So we consider, again, the decay ODE, which has exponentially decaying this solution. 112 00:07:05,565 --> 00:07:08,925 Tending to zero, ST tends to infinity. 113 00:07:09,080 --> 00:07:14,840 And if lambda is negative with a large modulus, then the decay is really fast. 114 00:07:14,995 --> 00:07:18,515 Now let's see what happens when we apply explicit Euler. 115 00:07:18,670 --> 00:07:21,470 So let me write down the recursion. 116 00:07:21,625 --> 00:07:24,265 This is the familiar one. 117 00:07:24,420 --> 00:07:30,740 And now the f of y of k is just lambda. 118 00:07:30,895 --> 00:07:33,615 So we can write this in this form. 119 00:07:33,770 --> 00:07:38,810 Single step and from it, well I should add that h is the step size. 120 00:07:38,965 --> 00:07:47,365 And from it we can read off a formula for the states in the sequence generated by the explicit Euler method. 121 00:07:47,520 --> 00:07:52,400 It's just the powers of this one plus hk term. 122 00:07:52,555 --> 00:07:54,475 And here's the initial value. 123 00:07:54,630 --> 00:07:58,390 Okay, oh not hk, h lambda. Sorry, lambda negative. 124 00:07:58,545 --> 00:08:01,825 Now what will happen as we send k to zero? 125 00:08:01,980 --> 00:08:04,940 This means as we consider large times. 126 00:08:05,095 --> 00:08:07,495 Well, we have to distinguish two cases. 127 00:08:07,650 --> 00:08:15,570 If the modulus of one plus h lambda is bigger than one, then it goes off to infinity. 128 00:08:15,725 --> 00:08:22,845 If, well, this is smaller than one, this modulus, then the sequence converges to zero. 129 00:08:23,000 --> 00:08:30,520 Actually this convergence to zero is what we want because this matches the behavior of the exact solution. 130 00:08:30,675 --> 00:08:34,515 However, this is tied to a condition on h. 131 00:08:34,670 --> 00:08:37,070 Now remember that lambda is negative. 132 00:08:37,225 --> 00:08:43,865 So in one plus h lambda, the modulus can only be smaller than one if h, time step size, 133 00:08:44,020 --> 00:08:47,380 is smaller than two over the modulus of lambda. 134 00:08:47,535 --> 00:08:54,175 And if this is violated, h is larger than that barrier, then we get this exponential blowup. 135 00:08:54,330 --> 00:09:00,970 Summing up, we incur blowup unless the time step size is small enough. 136 00:09:01,125 --> 00:09:08,245 And this, without observing this barrier, we get a, well, useless solution. 137 00:09:08,400 --> 00:09:14,440 So we could call this stability induced time step constraint because the occurrence of 138 00:09:14,595 --> 00:09:17,755 blowup is also called an instability phenomenon. 139 00:09:17,910 --> 00:09:28,550 Now, now you've understood why explicit Euler fails if you use a large time step and apply it to a decay equation with rapid decay. 140 00:09:28,705 --> 00:09:35,105 In the first example where we studied the behavior of the adaptive explicit Runge Kutta time stepper, 141 00:09:35,260 --> 00:09:36,740 we did not observe blowup. 142 00:09:36,895 --> 00:09:38,735 So there was something else going on. 143 00:09:38,890 --> 00:09:45,410 So now let's revisit adaptive time stepping and revisit it for this simple decay ODE. 144 00:09:45,565 --> 00:09:51,765 Here we use a rather simple time stepping scheme where we understand the components. 145 00:09:51,920 --> 00:10:00,800 It was a toy adaptive time stepping method based on the explicit Euler method as a base 146 00:10:00,955 --> 00:10:06,955 scheme and the explicit trapezoidal method of order two for local error estimation. 147 00:10:07,110 --> 00:10:11,790 But in any case, this involves explicit Euler and so we expect problems. 148 00:10:11,945 --> 00:10:13,225 Let's see what we get. 149 00:10:13,380 --> 00:10:15,500 Same as before. 150 00:10:15,655 --> 00:10:17,975 In fact, we do not observe blowup. 151 00:10:18,130 --> 00:10:20,130 Time step control avoids blowup. 152 00:10:20,285 --> 00:10:21,605 But at what price? 153 00:10:21,760 --> 00:10:26,000 The price is that time steps have to stay tiny. 154 00:10:26,155 --> 00:10:31,835 More precisely, the time step have to stay so small that the stability induced time step 155 00:10:31,990 --> 00:10:34,530 constraint is always satisfied. 156 00:10:34,685 --> 00:10:41,665 So step size control, in a sense, is not a complete failure because it manages to avoid blowup. 157 00:10:41,820 --> 00:10:49,000 It keeps the time step below the limit or below the threshold that is important to ensure stability. 158 00:10:49,155 --> 00:10:51,335 So it's moderately successful. 159 00:10:51,490 --> 00:10:57,690 But the price is that many, many tiny time steps have to be executed. 160 00:10:57,845 --> 00:10:59,725 And this, of course, not desirable. 161 00:10:59,880 --> 00:11:05,400 Now after we have sort of understood what's going on in the case of explicit Euler, let's 162 00:11:05,555 --> 00:11:10,795 go beyond that very simple method and look at a more complicated one. 163 00:11:10,950 --> 00:11:17,570 And that more complicated method is the explicit trapezoidal single step method that we constructed 164 00:11:17,725 --> 00:11:20,905 in the previous chapter and that enjoys order 2. 165 00:11:21,060 --> 00:11:23,280 We already derived the equations. 166 00:11:23,435 --> 00:11:29,035 So these are the formulae written in the form of a Runge Kutta method with increments. 167 00:11:29,190 --> 00:11:32,630 And now we apply it to the decay equation. 168 00:11:32,785 --> 00:11:33,545 No problem. 169 00:11:33,700 --> 00:11:36,280 Now f of y is just lambda y. 170 00:11:36,435 --> 00:11:40,155 And so now we can elaborate the formulas. 171 00:11:40,310 --> 00:11:42,110 So let me write down. 172 00:11:42,265 --> 00:11:46,465 So this K1 is then just lambda y naught. 173 00:11:46,620 --> 00:11:49,140 So we look at the first step. 174 00:11:49,295 --> 00:11:57,055 K2 from above is lambda y naught plus a lambda K1 lambda y naught. 175 00:11:57,210 --> 00:11:58,210 Plug it in. 176 00:11:58,365 --> 00:11:59,505 This is what we get. 177 00:11:59,660 --> 00:12:07,080 And so we find that y1 is 1 plus lambda h plus 1 half. 178 00:12:07,235 --> 00:12:08,555 This comes from here. 179 00:12:08,710 --> 00:12:10,750 Lambda h squared and y naught. 180 00:12:10,905 --> 00:12:12,025 I'm doing this quickly. 181 00:12:12,180 --> 00:12:15,100 This is elementary manipulations. 182 00:12:15,255 --> 00:12:18,495 And now let's abbreviate this number. 183 00:12:18,650 --> 00:12:22,730 Let's call it for reasons that we will become clear in a moment. 184 00:12:22,885 --> 00:12:27,165 S of lambda h is a sum function of lambda h. 185 00:12:27,320 --> 00:12:32,640 With this abbreviation, we can now write all the members of the sequence generated by the 186 00:12:32,795 --> 00:12:34,835 explicit trapezoidal method. 187 00:12:34,990 --> 00:12:41,030 It's this number s of lambda h to the power k when applied to the initial value. 188 00:12:41,185 --> 00:12:42,665 And this goes on forever. 189 00:12:42,820 --> 00:12:46,100 And now let's look at this function s. 190 00:12:46,255 --> 00:12:51,615 Because it's clear, in order to avoid blowup, we need as a necessary condition, necessary 191 00:12:51,770 --> 00:12:57,130 and sufficient condition, that the modulus of s of lambda h is below 1. 192 00:12:57,285 --> 00:12:59,685 And now let's look where this is the case. 193 00:12:59,840 --> 00:13:00,840 Sorry. 194 00:13:00,995 --> 00:13:03,355 This gives you the graph of this function s. 195 00:13:03,510 --> 00:13:07,870 Now I have abbreviated z for lambda h. 196 00:13:08,025 --> 00:13:12,065 And so you see that z has to be in the green region. 197 00:13:12,220 --> 00:13:16,660 And this means that lambda h has to be between minus 2 and 0. 198 00:13:16,815 --> 00:13:19,735 Lambda is negative, so the sign is not a problem. 199 00:13:19,890 --> 00:13:25,290 But the minus 2 bound is a problem because it forces h to be smaller than 2 over lambda. 200 00:13:25,445 --> 00:13:33,125 So this s function tells a lot about the occurrence of blowup, the violation of stability. 201 00:13:33,280 --> 00:13:36,480 And this is why it was dubbed stability function. 202 00:13:36,635 --> 00:13:41,475 And this is the stability function for the explicit trapezoidal rule. 203 00:13:41,630 --> 00:13:48,070 So now it seems that even for more complicated methods, it's possible to understand the onset 204 00:13:48,225 --> 00:13:54,465 of instability and to understand what kind of stability induced time step constraints 205 00:13:54,620 --> 00:13:55,620 we will face. 206 00:13:55,775 --> 00:14:01,775 So now let's tackle, well, already rather generous Runge Kutta single step methods. 207 00:14:01,930 --> 00:14:02,930 All right. 208 00:14:03,085 --> 00:14:08,205 So we study an s stage Runge Kutta single step method given through its Butcher tableau 209 00:14:08,360 --> 00:14:11,440 with an s by s, so called Butcher matrix. 210 00:14:11,595 --> 00:14:12,595 This is the a here. 211 00:14:12,750 --> 00:14:17,090 We call the definition of the method based on the Butcher coefficient. 212 00:14:17,245 --> 00:14:19,225 So these are the formulas. 213 00:14:19,380 --> 00:14:23,620 And now we apply this method as before to the decay ODE. 214 00:14:23,775 --> 00:14:28,775 So we apply to the case n equals 1, so one dimensional. 215 00:14:28,930 --> 00:14:34,970 And then we just plug f of y is equal to lambda y into the increment equations. 216 00:14:35,125 --> 00:14:41,525 And what is written here is a direct consequence of these formulas for the special case of 217 00:14:41,680 --> 00:14:43,640 the scalar decay ODE. 218 00:14:43,795 --> 00:14:47,875 Well, I've written bold one for this vector of all ones. 219 00:14:48,030 --> 00:14:51,430 Then there's a set of s before. 220 00:14:51,585 --> 00:14:56,125 I abbreviate lambda times h as set. 221 00:14:56,280 --> 00:14:59,000 And there's this k vector. 222 00:14:59,155 --> 00:15:04,775 And this is just a column vector of all the increments, the increments in this simple 223 00:15:04,930 --> 00:15:07,210 case n equals 1 are numbers. 224 00:15:07,365 --> 00:15:09,965 So collect the increments here. 225 00:15:10,120 --> 00:15:17,640 ks, this is a vector with s components, and well, there's a 1 over lambda here too. 226 00:15:17,795 --> 00:15:27,235 So well, this arrow is, say, not really tedious, elementary matrix vector manipulations. 227 00:15:27,390 --> 00:15:33,430 I wanted to rewrite it in the form of this linear system of equations. 228 00:15:33,585 --> 00:15:40,145 And the reason is that for linear systems, we have a so called determinant formula for 229 00:15:40,300 --> 00:15:41,020 the solution. 230 00:15:41,175 --> 00:15:43,055 And this is the so called Cramer rule. 231 00:15:43,210 --> 00:15:46,410 You may have heard about it in linear algebra or even in school. 232 00:15:46,565 --> 00:15:48,685 I'm going to give you the formula here. 233 00:15:48,840 --> 00:15:56,960 So the Cramer rule tells us this formula for the last component y1 of the solution of the 234 00:15:57,115 --> 00:15:58,115 linear system. 235 00:15:58,270 --> 00:16:03,390 And now the Cramer rule is convenient here because the denominator here is particularly 236 00:16:03,545 --> 00:16:07,825 simple, because this determinant is just equal to 1. 237 00:16:07,980 --> 00:16:08,500 Why? 238 00:16:08,655 --> 00:16:15,135 Because for an explicit Runge Kutta method, the Butcher matrix A is strictly lower triangular. 239 00:16:15,290 --> 00:16:21,130 So we compute the determinant here of a matrix which is lower triangular and has all 1s on 240 00:16:21,285 --> 00:16:22,285 the diagonal. 241 00:16:22,440 --> 00:16:27,140 And you know that the determinant of a triangular matrix is the product of its diagonal entries. 242 00:16:27,295 --> 00:16:28,295 So we get 1 here. 243 00:16:28,450 --> 00:16:33,750 So all that matters here is the numerator and this is this determinant, as one can see 244 00:16:33,905 --> 00:16:38,445 by expanding the determinant as you have learned it in linear algebra. 245 00:16:38,600 --> 00:16:43,700 So this gives us the following theorem that I have already displayed here about the so 246 00:16:43,855 --> 00:16:49,655 called stability function of an arbitrary explicit Runge Kutta method given by its Butcher 247 00:16:49,810 --> 00:16:50,790 scheme. 248 00:16:50,945 --> 00:16:52,545 It's this determinant. 249 00:16:52,700 --> 00:16:58,220 And as before, the stability function is the function that when you plug in lambda times 250 00:16:58,375 --> 00:17:04,555 h gives the number that, well, a state has to be multiplied with to get the next state. 251 00:17:04,710 --> 00:17:08,390 So what can we tell from this formula? 252 00:17:08,545 --> 00:17:15,585 The determinant of a matrix can be computed as a product of some of its entries, a sum 253 00:17:15,740 --> 00:17:18,180 of products of its entries to be precise. 254 00:17:18,335 --> 00:17:22,615 Here the entries of this matrix, they are all linear in set. 255 00:17:22,770 --> 00:17:28,490 So this is a sum of products of terms that are linear in set and the products in the 256 00:17:28,645 --> 00:17:35,705 case of an S by S matrix, which is the case here, will never contain more than S factors. 257 00:17:35,860 --> 00:17:42,120 So with this argument, the determinant formula immediately tells us that this is a polynomial 258 00:17:42,275 --> 00:17:44,195 of degree smaller equal S. 259 00:17:44,350 --> 00:17:50,790 So corollary, S state explicit Runge Kutta single step methods have a stability function, 260 00:17:50,945 --> 00:17:54,425 which is a polynomial of degree smaller equal S. 261 00:17:54,580 --> 00:18:00,700 Well, now you may be curious how these polynomials look like, these stability functions. 262 00:18:00,855 --> 00:18:06,435 Well, we already saw the stability functions for the simplest Runge Kutta method, the explicit 263 00:18:06,590 --> 00:18:07,090 Euler method. 264 00:18:07,245 --> 00:18:09,245 It's just 1 plus Z. 265 00:18:09,400 --> 00:18:14,520 And for the explicit trapezoidal rule, we computed it, 1 plus Z plus 1 half Z squared. 266 00:18:14,675 --> 00:18:20,655 And now, of course, the formula that we have just derived would have been a great shortcut 267 00:18:20,810 --> 00:18:22,490 to arrive at these formulas. 268 00:18:22,645 --> 00:18:28,245 We can apply now the determinant formula to any Runge Kutta method. 269 00:18:28,400 --> 00:18:33,620 For instance, the classical fourth order Runge Kutta method, and then we get such a polynomial. 270 00:18:33,775 --> 00:18:38,915 Then you see a pattern, 1 plus Z, 1 plus Z plus 1 half Z squared. 271 00:18:39,070 --> 00:18:42,050 And then for fourth order, this polynomial. 272 00:18:42,205 --> 00:18:43,265 Maybe maybe not. 273 00:18:43,420 --> 00:18:44,420 I'm going to tell you. 274 00:18:44,575 --> 00:18:47,215 There is a connection with the exponential function. 275 00:18:47,370 --> 00:18:53,750 And this connection is by no means surprising, because the stability function was derived 276 00:18:53,905 --> 00:18:58,665 by applying the Runge Kutta method to the decay ODE or growth ODE. 277 00:18:58,820 --> 00:19:01,380 So what does the Runge Kutta single step method do? 278 00:19:01,535 --> 00:19:04,135 Well, the first step, it computes y1. 279 00:19:04,290 --> 00:19:09,150 Let me highlight the role of the stability function applied to y0. 280 00:19:09,305 --> 00:19:13,665 And now let's track the exact solution over a time step of size h. 281 00:19:13,820 --> 00:19:18,020 So y of h, y of t, the exact trajectory. 282 00:19:18,175 --> 00:19:24,655 And now this is again multiplication of a number onto the initial value. 283 00:19:24,810 --> 00:19:32,990 And now you see that in fact, if our single step method gives us good accuracy, then the 284 00:19:33,145 --> 00:19:39,405 stability function should be a reasonable approximation of the exponential function. 285 00:19:39,560 --> 00:19:47,560 So we in fact expect that S of lambda h approximate the exponential of lambda h. 286 00:19:47,715 --> 00:19:50,535 And now recall the exponential series. 287 00:19:50,690 --> 00:19:55,230 So lambda h to the power L. And now replacing lambda h with Z. 288 00:19:55,385 --> 00:20:01,965 Well, what we see above is that the polynomials are truncated exponential series. 289 00:20:02,120 --> 00:20:07,880 And this is not a coincidence, this is even necessary once you know the order of the single 290 00:20:08,035 --> 00:20:09,035 step method. 291 00:20:09,190 --> 00:20:11,790 Because there's a rather general relationship. 292 00:20:11,945 --> 00:20:18,385 When you have an order q single step method, then the stability function, the first q plus 293 00:20:18,540 --> 00:20:24,100 is a polynomial, and the first q plus one power monomials are just the terms in the 294 00:20:24,255 --> 00:20:25,695 exponential series. 295 00:20:25,850 --> 00:20:27,270 Then there may be other terms. 296 00:20:27,425 --> 00:20:31,605 For the methods we have seen above, these terms are actually zero. 297 00:20:31,760 --> 00:20:34,520 And this is of course not surprising. 298 00:20:34,675 --> 00:20:38,595 It's even necessary for a particular order. 299 00:20:38,750 --> 00:20:43,830 And let me explain you why and state this explanation in the form of a lemma. 300 00:20:43,985 --> 00:20:50,705 Because we know that for order q method, the one step error decays as q plus first power 301 00:20:50,860 --> 00:20:52,400 of the step size. 302 00:20:52,555 --> 00:20:57,175 And what we have here, this is the one step error for the decay equation. 303 00:20:57,330 --> 00:21:04,850 So it must decay like order q plus one of the step size. 304 00:21:05,005 --> 00:21:08,565 So this is the connection to the order of the method. 305 00:21:08,720 --> 00:21:16,160 And now since we know that for a s stage explicit Runge Kutta method, s of z is a polynomial 306 00:21:16,315 --> 00:21:21,195 of degree s, and that a polynomial of degree s can approximate the exponential function 307 00:21:21,350 --> 00:21:28,390 in zero only up to order s plus one, then it's clear that an explicit s stage Runge 308 00:21:28,545 --> 00:21:32,425 Kutta single step method cannot go beyond order s. 309 00:21:32,580 --> 00:21:38,420 Well, there's no polynomial fitting the exponential function with a better accuracy. 310 00:21:38,575 --> 00:21:45,255 Now this was a bit of an excursion into the meaning of stability functions and its relationship 311 00:21:45,410 --> 00:21:47,390 with the exponential function. 312 00:21:47,545 --> 00:21:55,105 But our original intention with which we started the stability function was to understand stability 313 00:21:55,260 --> 00:21:57,560 induced time step constraints. 314 00:21:57,715 --> 00:21:59,835 Now let's return to this goal. 315 00:21:59,990 --> 00:22:07,490 Well, since the stability function evaluated for h times lambda is this multiplying factor 316 00:22:07,645 --> 00:22:14,485 in every time step, of course, now it's clear when exponential blow up will occur and when 317 00:22:14,640 --> 00:22:16,180 it won't. 318 00:22:16,335 --> 00:22:22,515 It will occur whenever the modulus of that value of the stability function is larger 319 00:22:22,670 --> 00:22:23,670 than one. 320 00:22:23,825 --> 00:22:29,045 And now recall that for an explicit Runge Kutta method, the stability function is a 321 00:22:29,200 --> 00:22:33,720 polynomial of degree s and a non constant polynomial. 322 00:22:33,875 --> 00:22:36,795 And this is necessary for the sake of consistency. 323 00:22:36,950 --> 00:22:41,410 And then whenever you have a non constant polynomial and you plug in large arguments, 324 00:22:41,565 --> 00:22:45,205 then it will become arbitrarily large in modulus. 325 00:22:45,360 --> 00:22:51,700 So the set of arguments for which the stability function has the good value, so a value with 326 00:22:51,855 --> 00:22:56,275 modulus smaller than one, which does not lead to exponential blow up, will necessarily be 327 00:22:56,430 --> 00:22:57,430 bounded. 328 00:22:57,585 --> 00:23:05,805 And this means, recalling that set is h lambda, that we can avoid exponential blow up only 329 00:23:05,960 --> 00:23:09,440 if the time step is sufficiently small. 330 00:23:09,595 --> 00:23:11,875 We cannot tamper with lambda. 331 00:23:12,030 --> 00:23:13,870 This is given in the OD. 332 00:23:14,025 --> 00:23:16,185 All right, so let's summarize. 333 00:23:16,340 --> 00:23:22,620 So exponential blow up for explicit Runge Kutta single step method applied to the simple scalar 334 00:23:22,775 --> 00:23:29,275 linear decay OD can only be avoided if lambda times h is sufficiently small. 335 00:23:29,430 --> 00:23:33,230 So since lambda is fixed, if the time step is sufficiently small. 336 00:23:33,385 --> 00:23:39,265 And this is, well, just to repeat it again, a stability induced time step constraint. 337 00:23:39,420 --> 00:23:45,280 So in particular, when lambda is negative and very big in modulus, then, well, our time 338 00:23:45,435 --> 00:23:49,535 step has to stay tiny just to avoid blow up. 339 00:23:49,690 --> 00:23:51,830 And this will mean many time steps. 340 00:23:51,985 --> 00:23:57,585 So it will have a really negative impact on the efficiency of the method. 341 00:23:57,740 --> 00:24:01,820 So now what if lambda is positive? 342 00:24:01,975 --> 00:24:05,555 Then the exact solution already shows exponential blow up. 343 00:24:05,710 --> 00:24:10,250 So an exponential blow up of the numerical solution is harmless. 344 00:24:10,405 --> 00:24:16,005 It's even desired because we want to capture the behavior of the exact solution. 345 00:24:16,160 --> 00:24:20,280 This is why here the focus is on negative lambda on the decay case. 346 00:24:20,435 --> 00:24:27,135 Nevertheless, you may object that this simple linear decay equation is not really interesting 347 00:24:27,290 --> 00:24:31,830 because, well, we have an analytic solution, the exponential function. 348 00:24:31,985 --> 00:24:39,305 Let me convince you that understanding this very simple scalar problem is important nevertheless 349 00:24:39,460 --> 00:24:43,260 because it tells us everything about more complicated cases. 350 00:24:43,415 --> 00:24:49,295 So let's look at a more complicated case, namely systems of ordinary differential equations, 351 00:24:49,450 --> 00:24:50,450 still linear. 352 00:24:50,605 --> 00:24:56,885 So y dot is my, y now vector with n components and m an n by n matrix. 353 00:24:57,040 --> 00:25:02,400 And let's again start with an example, example from circuit modeling. 354 00:25:02,555 --> 00:25:08,435 This is a very well known simple circuit, a so called RLC circuit. 355 00:25:08,590 --> 00:25:16,230 And while electrical engineers have an approach called nodal analysis and it yields second 356 00:25:16,385 --> 00:25:22,305 order differential equation for the voltage in this central node. 357 00:25:22,460 --> 00:25:27,480 We know how to convert this into a first order equation, introducing the first derivative 358 00:25:27,635 --> 00:25:30,415 of u as an auxiliary variable. 359 00:25:30,570 --> 00:25:32,410 This is what we do now. 360 00:25:32,565 --> 00:25:34,285 Auxiliary variable here called v. 361 00:25:34,440 --> 00:25:43,000 So now we have a 2D state space and we have this simple here non autonomous linear differential 362 00:25:43,155 --> 00:25:44,155 equation. 363 00:25:44,310 --> 00:25:50,830 And, well, from electrical engineering, we also know that the parameters that we have 364 00:25:50,985 --> 00:25:56,305 introduced, parameters in the model, well, they will be very large both and beta will 365 00:25:56,460 --> 00:25:57,740 be even larger than alpha. 366 00:25:57,895 --> 00:25:59,775 So this is a very usual setting. 367 00:25:59,930 --> 00:26:07,410 So now it can be written like this as an affine linear ODE and then we can apply our workhorse 368 00:26:07,565 --> 00:26:16,525 explicit adaptive time stepper ODE 4.5 and the result is as disappointing as in the beginning 369 00:26:16,680 --> 00:26:21,560 when we saw its well poor performance for the logistic ODE. 370 00:26:21,715 --> 00:26:22,715 Same here. 371 00:26:22,870 --> 00:26:28,230 This circuit is highly damped with the values for alpha and beta that were chosen. 372 00:26:28,385 --> 00:26:34,065 And so after no time it reaches a steady state so nothing happens anymore. 373 00:26:34,220 --> 00:26:37,700 Nevertheless the time stepper needs a lot of time steps. 374 00:26:37,855 --> 00:26:43,495 So time step control again enforces very small time step and now we already suspect that 375 00:26:43,650 --> 00:26:46,130 this has to do with some stability problems. 376 00:26:46,285 --> 00:26:52,685 Now let's try to understand this rigorously what's going on as we manage to understand 377 00:26:52,840 --> 00:26:55,800 it for the scalar case using the stability function. 378 00:26:55,955 --> 00:27:01,375 The key technique is diagonalization and this is what you've already seen. 379 00:27:01,530 --> 00:27:07,590 You've seen diagonalization when you were taught how to compute the solution of a linear 380 00:27:07,745 --> 00:27:10,945 ODE of a more precisely a system of this form. 381 00:27:11,100 --> 00:27:15,420 So we assume that the matrix M can be diagonalized. 382 00:27:15,575 --> 00:27:21,655 There's a regular matrix N by N which though N is real may be complex. 383 00:27:21,810 --> 00:27:27,170 So that many matrices can only be diagonalized in the complex numbers. 384 00:27:27,325 --> 00:27:29,685 So we don't rule out complex matrices here. 385 00:27:29,840 --> 00:27:36,800 V can be a complex matrix and so can be the eigenvalues that are on the diagonal of that 386 00:27:36,955 --> 00:27:44,755 matrix D, lambda 1 up to lambda N. And now the diagonalization technique is fairly standard. 387 00:27:44,910 --> 00:27:46,730 Let me just repeat it. 388 00:27:46,885 --> 00:27:48,745 You've seen it in linear algebra probably. 389 00:27:48,900 --> 00:27:50,420 You start from the equation. 390 00:27:50,575 --> 00:27:53,855 You do a substitution, so introduce a new unknown. 391 00:27:54,010 --> 00:27:57,610 Z of t is V inverse y of t. 392 00:27:57,765 --> 00:28:04,005 And well, when you do the math, simple manipulations, then you end up with an ODE for the new unknown 393 00:28:04,160 --> 00:28:05,880 function Z of this form. 394 00:28:06,035 --> 00:28:09,035 Matrix is now D. And of course D is diagonal. 395 00:28:09,190 --> 00:28:11,910 So now here we have decoupled scalar ODEs. 396 00:28:12,065 --> 00:28:20,465 Of course this is now when you denote by a subscript index the components of this Z function. 397 00:28:20,620 --> 00:28:24,780 And this is Zm dot lambda m Zm, m from 1 to N. 398 00:28:24,935 --> 00:28:26,695 We know how to solve this. 399 00:28:26,850 --> 00:28:31,870 So Zm t is e to the lambda m t Zm 0. 400 00:28:32,025 --> 00:28:38,285 And then we can extract the original solution by undoing the transformation. 401 00:28:38,440 --> 00:28:43,800 So this is just a reminder that how this diagonalization technique works. 402 00:28:43,955 --> 00:28:47,075 Now we collect all the decoupled solutions. 403 00:28:47,230 --> 00:28:50,230 And this is from m1 to N. 404 00:28:50,385 --> 00:28:52,465 So this is diagonalization. 405 00:28:52,620 --> 00:28:58,220 I've presented it to you in some detail so that you now can appreciate that the very 406 00:28:58,375 --> 00:29:01,835 same idea can be applied to the Runge Kutta method. 407 00:29:01,990 --> 00:29:07,610 But before we do this, let's first tackle the simplest case, the explicit Euler method. 408 00:29:07,765 --> 00:29:10,025 Because then it will be very transparent. 409 00:29:10,180 --> 00:29:15,080 So we start from the ODE and I write down the explicit Euler recursion. 410 00:29:15,235 --> 00:29:16,235 This one. 411 00:29:16,390 --> 00:29:19,430 Again, I assume that the matrix M can be diagonalized as above. 412 00:29:19,585 --> 00:29:22,305 Complex matrix V is admitted. 413 00:29:22,460 --> 00:29:26,220 And we now introduce transformed states. 414 00:29:26,375 --> 00:29:30,455 We multiply the Euler recursion by V inverse. 415 00:29:30,610 --> 00:29:38,390 And we see that this matrix will become D. And so, well, in analogy to the substitution 416 00:29:38,545 --> 00:29:42,825 made above, I now call this transformed state Zk. 417 00:29:42,980 --> 00:29:48,040 And I end up with the formula for the Z recursion. 418 00:29:48,195 --> 00:29:55,315 And as before, this is not surprising, thanks to D being diagonal, these are just decoupled 419 00:29:55,470 --> 00:29:56,470 scalar recursions. 420 00:29:56,625 --> 00:30:01,305 So I can now write this decoupled scalar recursion. 421 00:30:01,460 --> 00:30:07,780 And well, there are scalar recursions for the explicit Euler methods applied to scalar 422 00:30:07,935 --> 00:30:10,355 linear growth or decay equations. 423 00:30:10,510 --> 00:30:15,410 So this is explicit Euler for y dot is lambda L y. 424 00:30:15,565 --> 00:30:20,885 Now this is interesting because by this simple transformation, we can reduce the explicit 425 00:30:21,040 --> 00:30:28,160 Euler method for the linear system to explicit Euler method applied to scalar equations. 426 00:30:28,315 --> 00:30:37,195 So and now it's clear, since V is a constant matrix, if the Zk sequence experiences a blowup, 427 00:30:37,350 --> 00:30:38,970 so will the yk sequence. 428 00:30:39,125 --> 00:30:43,165 And so the following sentence can summarize what we have observed. 429 00:30:43,320 --> 00:30:49,300 The explicit Euler method generates uniformly bounded solution sequences for this linear 430 00:30:49,455 --> 00:30:55,495 autonomous equation, assuming that m can be diagonalized with eigenvalues lambda 1, lambda 431 00:30:55,650 --> 00:31:02,250 n, if and only if the same explicit Euler method, while there's only one, generates 432 00:31:02,405 --> 00:31:10,065 uniformly bounded sequences for all the scalar linear ODEs, scalar linear ODEs induced by 433 00:31:10,220 --> 00:31:11,320 the eigenvalues. 434 00:31:11,475 --> 00:31:19,715 So and we already know when this uniform boundedness of the sequences will occur for Euler equation 435 00:31:19,870 --> 00:31:22,710 applied to the scalar linear autonomous ODE. 436 00:31:22,865 --> 00:31:31,225 This will occur if 1 plus h lambda i has modulus smaller than 1, or modulus 1 would be okay 437 00:31:31,380 --> 00:31:32,460 for boundedness. 438 00:31:32,615 --> 00:31:37,735 Well now we have discovered the stability induced time step constraint for systems, 439 00:31:37,890 --> 00:31:40,690 and it's the same as for the scalar case. 440 00:31:40,845 --> 00:31:44,565 Remark, this is also true if the eigenvalues are complex. 441 00:31:44,720 --> 00:31:45,720 So same argument. 442 00:31:45,875 --> 00:31:53,075 I'd like to refer you to remark, say 12, oh sorry, I think I have this already. 443 00:31:53,230 --> 00:31:58,390 Now let's do this after the explicit Euler method for general Runge Kutta method. 444 00:31:58,545 --> 00:32:02,105 And in fact the reasoning follows exactly the same lines. 445 00:32:02,260 --> 00:32:07,100 So let's recall the formulas and let's assume the Runge Kutta method is given through its 446 00:32:07,255 --> 00:32:08,255 Butcher scheme. 447 00:32:08,410 --> 00:32:11,090 And now we do the same as for explicit Euler. 448 00:32:11,245 --> 00:32:17,545 We try to understand what happens when we apply the explicit Runge Kutta method to this 449 00:32:17,700 --> 00:32:20,520 linear autonomous ODE. 450 00:32:20,675 --> 00:32:25,035 We plug f of y is equal to m of y into the increment equations. 451 00:32:25,190 --> 00:32:26,270 This is what we get. 452 00:32:26,425 --> 00:32:28,765 And then do diagonalization. 453 00:32:28,920 --> 00:32:35,100 Of course assuming that the matrix M can be diagonalized as before with eigenvalues lambda 454 00:32:35,255 --> 00:32:39,095 1 up to lambda n, the eigenvectors collected in a matrix V. 455 00:32:39,250 --> 00:32:41,590 So, and now same idea. 456 00:32:41,745 --> 00:32:47,165 We transform the states and also the increments through applying V inverse. 457 00:32:47,320 --> 00:32:51,760 So we apply V inverse to the increment equation and get this. 458 00:32:51,915 --> 00:32:56,015 And well, you see that D will come into play, of course. 459 00:32:56,170 --> 00:32:59,110 And then we also apply the inverse to the second equation. 460 00:32:59,265 --> 00:33:02,065 That gives us the next state. 461 00:33:02,220 --> 00:33:07,740 Then introduce just notations for the transformed increments, transformed states. 462 00:33:07,895 --> 00:33:12,435 Again of course a set for V inverse applied to a y state. 463 00:33:12,590 --> 00:33:16,090 And then this equation becomes this one. 464 00:33:16,245 --> 00:33:21,245 And well, this is now decoupled because D is diagonal. 465 00:33:21,400 --> 00:33:24,800 And now we have again decoupled scalar equations. 466 00:33:24,955 --> 00:33:25,955 Same as before. 467 00:33:26,110 --> 00:33:27,150 Let me write them down. 468 00:33:27,305 --> 00:33:32,785 Say for the top equation in component wise form looks like this. 469 00:33:32,940 --> 00:33:37,560 The diagonal entry of D lambda m, the eigenvalue, comes into play. 470 00:33:37,715 --> 00:33:40,975 And of course the second equation is automatically decoupled. 471 00:33:41,130 --> 00:33:43,190 Just a vector operation. 472 00:33:43,345 --> 00:33:49,185 And what we have here, what these decoupled equations represent is just one step. 473 00:33:49,340 --> 00:33:55,000 So the first step of the same Runge Kutta single step method now applied to a scalar 474 00:33:55,155 --> 00:33:56,255 linear equation. 475 00:33:56,410 --> 00:33:59,230 The scalar linear equation y dot is lambda m y. 476 00:33:59,385 --> 00:34:03,305 So a decay or growth equation with a rate lambda m. 477 00:34:03,460 --> 00:34:05,320 The eigenvalue shows up here. 478 00:34:05,475 --> 00:34:06,655 So what does this tell us? 479 00:34:06,810 --> 00:34:12,070 Well, the same message holds as we have seen for the Euler method. 480 00:34:12,225 --> 00:34:19,685 So the Runge Kutta single step method will manage to avoid blow up for this autonomous 481 00:34:19,840 --> 00:34:20,840 linear ODE. 482 00:34:20,995 --> 00:34:25,635 We assume m to be diagonalizable with eigenvalues lambda 1, lambda n. 483 00:34:25,790 --> 00:34:34,190 If and only if it generates uniformly bounded sequences when applied to the scalar ODE's 484 00:34:34,345 --> 00:34:39,225 z dot is lambda i z where lambda i can be any eigenvalue. 485 00:34:39,380 --> 00:34:45,320 And we have seen that explicit Runge Kutta method will avoid exponential blow up for 486 00:34:45,475 --> 00:34:51,175 the scalar equations only if the time step is sufficiently small and the size of the 487 00:34:51,330 --> 00:34:55,890 time step is determined by the stability function, the threshold. 488 00:34:56,045 --> 00:35:00,125 So everything can in fact be reduced to the scalar case. 489 00:35:00,280 --> 00:35:07,120 So once we understand this behavior of Runge Kutta methods for the scalar linear autonomous 490 00:35:07,275 --> 00:35:14,515 ODE, we essentially understand what's going on also in the case of linear systems of ODE. 491 00:35:14,670 --> 00:35:20,730 And we can draw this conclusion I would like to point out because the y states will blow 492 00:35:20,885 --> 00:35:29,985 up if and only if these set components, if one of them will blow up for of course k to 493 00:35:30,140 --> 00:35:33,420 infinity, k to infinity for one m. 494 00:35:33,575 --> 00:35:39,695 So blow up in one component will lead to a blow up of the whole solution. 495 00:35:39,850 --> 00:35:44,450 So let's summarize this again and I do it in the form of a theorem. 496 00:35:44,605 --> 00:35:49,405 So an explicit Runge Kutta single step method when applied to the linear autonomous ODE 497 00:35:49,560 --> 00:35:55,720 y dot is m of y can avoid blow up only if its stability function has modulus smaller 498 00:35:55,875 --> 00:36:02,755 than one for all values lambda i h where lambda i runs through the whole set of eigenvalues 499 00:36:02,910 --> 00:36:03,910 of m. 500 00:36:04,065 --> 00:36:10,265 Of course for our derivation we assume that m is diagonalizable but this result also carries 501 00:36:10,420 --> 00:36:12,420 over to arbitrary matrices. 502 00:36:12,575 --> 00:36:18,215 And so key to understanding the stability induced time step constraint is again here 503 00:36:18,370 --> 00:36:19,990 the stability function. 504 00:36:20,145 --> 00:36:26,845 And now let's finally characterize the values for h for which we can be sure that no blow 505 00:36:27,000 --> 00:36:28,040 up will occur. 506 00:36:28,195 --> 00:36:32,195 And this is connected to the concept of a region of stability. 507 00:36:32,350 --> 00:36:39,430 And well what I'm telling you now is after the whole discussion very natural, very clear 508 00:36:39,585 --> 00:36:40,585 I hope. 509 00:36:40,740 --> 00:36:46,980 So we call region of absolute stability of a single step method the complex numbers for 510 00:36:47,135 --> 00:36:51,095 which its stability function has modulus smaller than one. 511 00:36:51,250 --> 00:36:57,250 And then well it's clear that if lambda h belongs to this region of stability then we 512 00:36:57,405 --> 00:37:04,445 can be sure that no blow up of solutions while of course computed with the Runge Kutta single 513 00:37:04,600 --> 00:37:08,680 step method for the ODE y dot is lambda y can occur. 514 00:37:08,835 --> 00:37:14,515 So now we also know that for explicit s stage Runge Kutta single step method the stability 515 00:37:14,670 --> 00:37:16,410 function is a polynomial. 516 00:37:16,565 --> 00:37:23,745 For polynomial s of psi will inevitably be bounded and this means that the above condition 517 00:37:23,900 --> 00:37:30,580 implies a time step constraint in the form that well lambda m h is has to belong to the 518 00:37:30,735 --> 00:37:39,615 stability region and for all indices m and for all eigenvalues of the matrix that defines 519 00:37:39,770 --> 00:37:49,290 this, well that defines for all eigenvalues that and the matrix defines the linear ODE. 520 00:37:49,445 --> 00:37:53,205 And then time step constraint is necessary to avoid blow up. 521 00:37:53,360 --> 00:37:59,720 One remark of course these eigenvalues can be complex but complex numbers can be treated 522 00:37:59,875 --> 00:38:05,635 here like real numbers when it comes to their impact on blow up or decay of solution. 523 00:38:05,790 --> 00:38:13,450 So now I hope we have obtained a rather detailed understanding of what goes wrong when you 524 00:38:13,605 --> 00:38:21,145 try to apply Runge Kutta's explicit Runge Kutta single step method to ODE based models 525 00:38:21,300 --> 00:38:27,660 that describe phenomena with strong decay and that this means a stability induced time 526 00:38:27,815 --> 00:38:32,815 step constraint just necessary to avoid exponential blow up. 527 00:38:32,970 --> 00:38:39,790 And now the next step will be to find methods that do not suffer from this stability induced 528 00:38:39,945 --> 00:38:40,225 time step constraint.