1 00:00:00,160 --> 00:00:04,000 Welcome back to chapter 12 of the numerical methods course. 2 00:00:04,159 --> 00:00:07,719 Well, by now we know what are stiff initial value problems. 3 00:00:07,878 --> 00:00:13,398 We also know how to tell whether a given initial value problem has this particular property 4 00:00:13,557 --> 00:00:14,717 of stiffness. 5 00:00:14,876 --> 00:00:20,976 And we have a rather thorough understanding of why stiff initial value problems give so 6 00:00:21,135 --> 00:00:25,715 much trouble to the usual explicit single step methods. 7 00:00:25,874 --> 00:00:32,034 Well, just imagine a client wants you to solve an initial value problem and you tell him 8 00:00:32,193 --> 00:00:38,153 or her, I'm sorry, your initial value problem is probably stiff and so I'm afraid I can't 9 00:00:38,312 --> 00:00:39,952 solve it efficiently. 10 00:00:40,111 --> 00:00:45,551 Even when using my very sophisticated adaptive time stepping schemes. 11 00:00:45,710 --> 00:00:50,270 Well, after you've watched this video, you will be in a much better situation. 12 00:00:50,429 --> 00:00:57,429 Then you will know what kind of single step methods can deal with stiff initial value 13 00:00:57,588 --> 00:00:58,588 problems. 14 00:00:58,747 --> 00:01:01,987 It's the class of implicit Runge Kutta single step methods. 15 00:01:02,146 --> 00:01:05,466 And I'm going to tell you the details right now. 16 00:01:05,625 --> 00:01:13,305 Remember that all the difficulties that arose from an initial value problem being stiff, 17 00:01:13,464 --> 00:01:19,824 the difficulties for numerical integration were somehow tied to the class of explicit 18 00:01:19,983 --> 00:01:21,663 single step methods. 19 00:01:21,822 --> 00:01:26,662 Most of numerical integrators that we have seen so far belong to this class, but not 20 00:01:26,821 --> 00:01:27,821 all of them. 21 00:01:27,980 --> 00:01:31,020 We know two that are not explicit. 22 00:01:31,179 --> 00:01:35,939 And the simplest that we have already seen is the implicit Euler method. 23 00:01:36,098 --> 00:01:37,358 Well, let's try it. 24 00:01:37,517 --> 00:01:43,297 Let's try it for our linear scalar model problem of the decay equation. 25 00:01:43,456 --> 00:01:50,376 In a numerical experiment, let's study how the implicit Euler method behaves and compare 26 00:01:50,535 --> 00:01:54,435 it to the behavior of the explicit Euler method. 27 00:01:54,594 --> 00:01:57,374 The behavior of the latter, well, is predictable. 28 00:01:57,533 --> 00:02:00,793 We have studied it on many occasions. 29 00:02:00,952 --> 00:02:08,012 And when the decay rate is big, applying the explicit Euler method with a not too small 30 00:02:08,171 --> 00:02:10,711 time step will give us blow up. 31 00:02:10,870 --> 00:02:16,250 As you see here in these arrow curves that I think I've already displayed earlier, for 32 00:02:16,409 --> 00:02:22,689 large time steps, we have blow up and this is how well stability induced time step constraint 33 00:02:22,848 --> 00:02:25,168 manifests themselves. 34 00:02:25,327 --> 00:02:27,527 Now let's apply the implicit Euler method. 35 00:02:27,686 --> 00:02:29,566 How to do it is straightforward. 36 00:02:29,725 --> 00:02:35,705 We know the formula of the recursion and we can plot the same arrow curves for the implicit 37 00:02:35,864 --> 00:02:36,864 Euler method. 38 00:02:37,023 --> 00:02:43,383 And to our surprise, its behavior is much better regardless of the decay rate and even 39 00:02:43,542 --> 00:02:48,882 regardless of the time step, even for very large time steps, no blow up. 40 00:02:49,041 --> 00:02:55,001 And as we reduce the time step, asymptotically we achieve first order convergence. 41 00:02:55,160 --> 00:03:00,960 So apparently the implicit Euler method is immune to instability. 42 00:03:01,119 --> 00:03:03,299 There's no time step constraint. 43 00:03:03,458 --> 00:03:09,078 And now let's try to achieve a deeper understanding of this amazing trait. 44 00:03:09,237 --> 00:03:14,397 Let's do our familiar linear model problem analysis that we have already conducted for 45 00:03:14,556 --> 00:03:17,196 the explicit Euler. 46 00:03:17,355 --> 00:03:20,595 Let's apply it to the implicit Euler method. 47 00:03:20,754 --> 00:03:23,894 So we start from the linear model problem. 48 00:03:24,053 --> 00:03:30,713 Decay equation, y dot is lambda y with negative lambda or at least negative real part so that 49 00:03:30,872 --> 00:03:34,832 the exact solution will decay to zero for large times. 50 00:03:34,991 --> 00:03:40,071 Recall the implicit Euler recursion now applied to this particular linear model problem. 51 00:03:40,230 --> 00:03:46,210 So then it's easy to rewrite, say, in an explicit way. 52 00:03:46,369 --> 00:03:47,369 Let's do this. 53 00:03:47,528 --> 00:03:50,088 So here we solve for yk plus one. 54 00:03:50,247 --> 00:03:57,487 And what we get is that yk plus one is then one over one minus h lambda yk. 55 00:03:57,646 --> 00:03:59,406 And this is the recursion. 56 00:03:59,565 --> 00:04:06,065 And this immediately yields a formula for the members of the sequence of states generated 57 00:04:06,224 --> 00:04:12,964 by the implicit Euler method, one minus one minus h lambda to the power k for the k state 58 00:04:13,123 --> 00:04:15,623 of course multiplied with the initial state. 59 00:04:15,782 --> 00:04:21,782 And now it becomes important that lambda is negative because this means now that the modulus 60 00:04:21,941 --> 00:04:28,121 of one over one minus h lambda is definitely smaller than one. 61 00:04:28,280 --> 00:04:33,620 So letting k tend to infinity, this will tend to zero. 62 00:04:33,779 --> 00:04:36,899 And it will do so for all positive h. 63 00:04:37,058 --> 00:04:39,058 So for all time step sizes. 64 00:04:39,217 --> 00:04:40,377 And now it's clear. 65 00:04:40,536 --> 00:04:41,956 There's no instability. 66 00:04:42,115 --> 00:04:46,815 It cannot happen for this simple scalar linear model problem. 67 00:04:46,974 --> 00:04:50,734 Now what about the next more complicated case? 68 00:04:50,893 --> 00:04:53,373 The case of a linear system of ODEs. 69 00:04:53,532 --> 00:04:58,132 y dot is m of y with a n by n matrix m. 70 00:04:58,291 --> 00:05:00,171 This is the initial value problem. 71 00:05:00,330 --> 00:05:05,730 We do exactly the same as for the analysis of the explicit Euler method. 72 00:05:05,889 --> 00:05:12,169 We rely on the diagonalization technique on expansion into eigenvectors of m. 73 00:05:12,328 --> 00:05:16,088 Of course under the assumption that m can be diagonalized. 74 00:05:16,247 --> 00:05:22,427 So let's assume this diagonalization and then we multiply the recursion of the implicit 75 00:05:22,586 --> 00:05:25,206 Euler method with v inverse. 76 00:05:25,365 --> 00:05:27,405 Call v inverse y. 77 00:05:27,564 --> 00:05:31,284 The set states, this is expansion into eigenvectors. 78 00:05:31,443 --> 00:05:38,243 And then we notice that the recursion now for the set case states is a decoupled recursion 79 00:05:38,402 --> 00:05:44,262 because what pops up here is the diagonal matrix D. D is diagonal. 80 00:05:44,421 --> 00:05:46,041 So we get decoupling. 81 00:05:46,200 --> 00:05:53,260 And now we have, we end up with n scalar recursions that we can solve as we did it in the case 82 00:05:53,419 --> 00:05:55,039 of the scalar model problem. 83 00:05:55,198 --> 00:06:01,078 So we immediately end up with this formula for the members of the set case sequence. 84 00:06:01,237 --> 00:06:09,717 And again we notice that if all eigenvalues of m have negative real part, all components 85 00:06:09,876 --> 00:06:16,756 of the set case state sequence will tend to zero regardless of the time step size. 86 00:06:16,915 --> 00:06:23,575 So if we have eigenvalues with negative real part throughout, then there's no instability 87 00:06:23,734 --> 00:06:25,514 regardless of the time step size. 88 00:06:25,673 --> 00:06:33,433 And this, well, irrelevance of the time step size in terms of stability is called unconditional 89 00:06:33,592 --> 00:06:34,592 stability. 90 00:06:34,751 --> 00:06:41,471 So unconditional stability in other words means there's no time step constraint that 91 00:06:41,630 --> 00:06:48,030 is due to stability issues in stark contrast to the explicit Rungicutta single step methods. 92 00:06:48,189 --> 00:06:51,669 All right, this is implicit Euler, a first order method. 93 00:06:51,828 --> 00:06:58,188 But now of course you know that higher order methods for time stepping are highly desirable. 94 00:06:58,347 --> 00:07:04,147 For explicit Rungicutta single step method, I showed you this bootstrapping trick to construct 95 00:07:04,306 --> 00:07:05,306 them. 96 00:07:05,465 --> 00:07:11,205 And now let me show you an approach that can come up with generalization of the implicit 97 00:07:11,364 --> 00:07:12,464 Euler method. 98 00:07:12,623 --> 00:07:14,603 And this is the collocation approach. 99 00:07:14,762 --> 00:07:21,222 So our goal is to generalize the implicit Euler method towards finding unconditionally 100 00:07:21,381 --> 00:07:24,821 stable single step methods of higher order. 101 00:07:24,980 --> 00:07:26,420 Well, this is the setting. 102 00:07:26,579 --> 00:07:34,539 Now let's consider a general ODE with a right hand side function f in an n dimensional state 103 00:07:34,698 --> 00:07:36,018 space. 104 00:07:36,177 --> 00:07:40,297 Let me introduce you to the already mentioned collocation approach. 105 00:07:40,456 --> 00:07:44,116 We aim for single step methods, so it's enough to specify the first step. 106 00:07:44,275 --> 00:07:51,135 So I will tell you how to get starting from the initial state, the first state in the 107 00:07:51,294 --> 00:07:55,174 sequence and then everything is again is applied again and again. 108 00:07:55,333 --> 00:08:00,493 And first step covers the time interval from t naught to t1. 109 00:08:00,652 --> 00:08:06,652 So now the fundamental idea is that we approximate the solution trajectory of the initial value 110 00:08:06,811 --> 00:08:13,211 problem with the initial state y naught by a function from a restricted function space. 111 00:08:13,370 --> 00:08:20,030 Here it's called V. So this V is referred to usually as the trial space. 112 00:08:20,189 --> 00:08:26,049 So this is the space from which we want to pick the functions that we use to approximate 113 00:08:26,208 --> 00:08:27,408 the trajectory. 114 00:08:27,567 --> 00:08:32,267 And a key property of the trial space is that it's finite dimensional. 115 00:08:32,426 --> 00:08:38,646 And here we assume that the trial space of these functions, these are vector value functions 116 00:08:38,805 --> 00:08:46,325 with n components, has dimension n times s plus 1 for a given integer s. 117 00:08:46,484 --> 00:08:47,484 And this is fixed. 118 00:08:47,643 --> 00:08:52,283 So then of course we need criteria for selecting the approximating function. 119 00:08:52,442 --> 00:08:57,542 And the selection will be based on so called collocation conditions. 120 00:08:57,701 --> 00:09:00,321 And this is the gist of the collocation approach. 121 00:09:00,480 --> 00:09:07,360 So we want to impose certain requirements on the function yh that we take from the trial 122 00:09:07,519 --> 00:09:13,839 space V. Of course we wanted to satisfy the initial condition, but we also wanted to satisfy 123 00:09:13,998 --> 00:09:14,998 the ODE. 124 00:09:15,157 --> 00:09:18,597 Well, satisfying the ODE is too much to ask. 125 00:09:18,756 --> 00:09:25,796 We wanted to satisfy the ODE in a few isolated points here called tau j. 126 00:09:25,955 --> 00:09:29,635 So this is the ODE at tau j. 127 00:09:29,794 --> 00:09:33,394 And this is to be satisfied pointwise. 128 00:09:33,553 --> 00:09:40,373 For well given tau j's in the interval, we start with tau 1 and tau 2. 129 00:09:40,532 --> 00:09:45,692 And we have s of these taus in the interval t0 to t1. 130 00:09:45,851 --> 00:09:48,271 They are given and they have a special name. 131 00:09:48,430 --> 00:09:50,710 They are called collocation points. 132 00:09:50,869 --> 00:09:56,829 So the gist of the collocation method is to demand that for a function from the trial 133 00:09:56,988 --> 00:10:03,908 space, the ODE must be satisfied pointwise in these s collocation points. 134 00:10:04,067 --> 00:10:09,607 And we hope that these collocations, collocation conditions fix a unique yh. 135 00:10:09,766 --> 00:10:15,786 And then of course the state at the end of the time step is just yh at t1. 136 00:10:15,945 --> 00:10:24,065 Well, practical collocation method, almost all of them use trial spaces of polynomials. 137 00:10:24,224 --> 00:10:29,584 And for natural reasons, because polynomials are easy to deal with. 138 00:10:29,743 --> 00:10:34,223 They are in particular easy to integrate and differentiate. 139 00:10:34,382 --> 00:10:40,702 So let's look at now at this special version of the collocation method, the polynomial 140 00:10:40,861 --> 00:10:45,541 collocation method, which I would advertise as the standard collocation method. 141 00:10:45,700 --> 00:10:48,060 So we use vector valued polynomials. 142 00:10:48,219 --> 00:10:55,179 And then for the sake of convenience, we normalize the collocation points to the interval 0,1 143 00:10:55,338 --> 00:10:56,538 by an affine mapping. 144 00:10:56,697 --> 00:11:04,817 So we define that tau j's in the interval t0 to t1 by this formula, tau t0 plus cjh, 145 00:11:04,976 --> 00:11:12,416 where now we have these parameters c1 up to cs between 0 and 1, while h is the time step 146 00:11:12,575 --> 00:11:13,855 size as usual. 147 00:11:14,014 --> 00:11:20,454 And now we need another tool, a tool from polynomial interpolation, namely Lagrange 148 00:11:20,613 --> 00:11:21,613 polynomials. 149 00:11:21,772 --> 00:11:28,932 And here I would like to ask you to flip back to chapter 5 and refresh yourself on Lagrange 150 00:11:29,091 --> 00:11:30,091 polynomials. 151 00:11:30,250 --> 00:11:32,970 Well, what are the Lagrange polynomials that we need? 152 00:11:33,129 --> 00:11:39,489 Call them as usual lj, but now they have a degree s minus 1 and we have s of them. 153 00:11:39,648 --> 00:11:45,388 And they are uniquely defined by their cardinal basis property with respect to the normalized 154 00:11:45,547 --> 00:11:50,767 collocation points given by these parameters c1 up to cs. 155 00:11:50,926 --> 00:11:56,946 So we have that lj at points c, little l, is delta lj, the Kronecker symbol. 156 00:11:57,105 --> 00:11:59,645 So this is the cardinal basis property. 157 00:11:59,804 --> 00:12:03,724 You may wonder, why do we use the degree s minus 1? 158 00:12:03,883 --> 00:12:09,163 There's a mismatch to the degree of the polynomials in the space V. Well, we use the degree s 159 00:12:09,322 --> 00:12:16,542 minus 1 because these Lagrange polynomials will be used to write down the derivative 160 00:12:16,701 --> 00:12:17,701 of yh. 161 00:12:17,860 --> 00:12:20,480 And that derivative has degree s minus 1. 162 00:12:20,639 --> 00:12:28,559 So let's write down yh dot the derivative in Lagrange polynomial representation, also 163 00:12:28,718 --> 00:12:31,058 normalized to the interval 0, 1. 164 00:12:31,217 --> 00:12:35,317 So we know this is a polynomial of degree s minus 1. 165 00:12:35,476 --> 00:12:40,636 It's in the span of the Lagrange polynomials and thanks to their cardinal basis property, 166 00:12:40,795 --> 00:12:43,355 we can just write it like this. 167 00:12:43,514 --> 00:12:48,894 The weights of the Lagrange polynomial are just nodal values of the y dot function. 168 00:12:49,053 --> 00:12:50,733 And this gives this formula. 169 00:12:50,892 --> 00:12:53,432 And now we can use the ODE. 170 00:12:53,591 --> 00:12:58,091 This is the, well, now we put in, in other words, the collocation conditions. 171 00:12:58,250 --> 00:13:00,890 So let me indicate this here. 172 00:13:01,049 --> 00:13:02,929 Collocations conditions enter here. 173 00:13:03,088 --> 00:13:08,468 So replacing this y dot value with the value of the right hand side. 174 00:13:08,627 --> 00:13:12,707 Well, you may object that this is maybe not so useful here. 175 00:13:12,866 --> 00:13:17,886 But anyway, it's not so useful because yh is unknown. 176 00:13:18,045 --> 00:13:21,705 So we don't know yh at this point. 177 00:13:21,864 --> 00:13:23,384 So we cannot evaluate f. 178 00:13:23,543 --> 00:13:25,783 So it seems we are turning in circles. 179 00:13:25,942 --> 00:13:31,022 But let's call this, these unknown quantities kj. 180 00:13:31,181 --> 00:13:32,741 This is a vector in Rn. 181 00:13:32,900 --> 00:13:36,180 And now we have a formula for yh dot. 182 00:13:36,339 --> 00:13:38,439 Let's recover yh by integration. 183 00:13:38,598 --> 00:13:39,598 Now integrate. 184 00:13:39,757 --> 00:13:47,517 So this will give us yh dot t naught plus, no, we have integrated, sorry, integrate with 185 00:13:47,676 --> 00:13:49,976 respect to psi. 186 00:13:50,135 --> 00:13:56,295 We remove the derivative, fundamental theorem of calculus, of course, then insert the initial 187 00:13:56,454 --> 00:13:57,454 value. 188 00:13:57,613 --> 00:14:03,313 And this h comes from, well, integration by substitution because there's a psi times h. 189 00:14:03,472 --> 00:14:05,452 And let's write down the formula. 190 00:14:05,611 --> 00:14:12,371 So here now we have to use this short notation kj and we integrate the Lagrange polynomial 191 00:14:12,530 --> 00:14:16,770 lj, theta, d theta. 192 00:14:16,929 --> 00:14:18,969 So this is now what we have. 193 00:14:19,128 --> 00:14:25,488 Well, and now we have a formula for, say, these unknown vectors. 194 00:14:25,647 --> 00:14:29,567 And now we can write down an equation for the kj's. 195 00:14:29,726 --> 00:14:31,406 So it's a bit tedious. 196 00:14:31,565 --> 00:14:36,165 But now we use the formula that we have obtained by integration, plug it into f. 197 00:14:36,324 --> 00:14:38,644 This gives us the k vector. 198 00:14:38,803 --> 00:14:42,703 And so what we have now is an equation for the kis. 199 00:14:42,862 --> 00:14:44,842 They occur here and here. 200 00:14:45,001 --> 00:14:53,841 Or more precisely, it's s equations for k1 up to ks, all vectors of length n. 201 00:14:54,000 --> 00:14:56,460 Now let's assume we solve this. 202 00:14:56,619 --> 00:15:00,519 Then we can write down an equation for the next state. 203 00:15:00,678 --> 00:15:03,478 Let me copy this. 204 00:15:03,637 --> 00:15:09,677 And this is the equation for the next state when I plug in 1 for psi because t0 plus h 205 00:15:09,836 --> 00:15:18,956 is t1 and yh at t1 is what we have defined as the next state, so the result of a single 206 00:15:19,115 --> 00:15:20,115 time step. 207 00:15:20,274 --> 00:15:23,894 Well, I had to go relatively fast through this collocation approach. 208 00:15:24,053 --> 00:15:27,093 But I hope you got the main ideas. 209 00:15:27,252 --> 00:15:33,352 Final approximation of the trajectory plus collocation conditions, which imply that we 210 00:15:33,511 --> 00:15:41,471 want the ODE to hold at a finite sequence of points in the current time interval. 211 00:15:41,630 --> 00:15:48,110 I've already mentioned, and the formulas are still on the left page, that we have to solve 212 00:15:48,269 --> 00:15:53,069 some equations to determine these unknown vectors kj. 213 00:15:53,228 --> 00:15:57,968 And this of course means that we face an implicit method. 214 00:15:58,127 --> 00:16:00,547 So let me call this star. 215 00:16:00,706 --> 00:16:10,466 This star, well, encodes an n times s system, usually nonlinear system of equations for 216 00:16:10,625 --> 00:16:15,065 the vectors kj, k1 up to ks in Rn. 217 00:16:15,224 --> 00:16:18,424 So this is why this is definitely implicit. 218 00:16:18,583 --> 00:16:26,383 So we have to solve nonlinear systems of equations in every time step to determine the next state. 219 00:16:26,542 --> 00:16:31,462 We can no longer simply plug in states into the right hand side vector field and get the 220 00:16:31,621 --> 00:16:33,661 next state in this simple way. 221 00:16:33,820 --> 00:16:39,060 Now let's see how such collocation methods perform in an experiment. 222 00:16:39,219 --> 00:16:45,659 We first look at a simple nonlinear scalar initial value problem for the logistic differential 223 00:16:45,818 --> 00:16:48,398 equations of resource limited growth. 224 00:16:48,557 --> 00:16:49,557 This is it. 225 00:16:49,716 --> 00:16:55,056 This is a classical example of a stiff problem because lambda is equal to 100. 226 00:16:55,215 --> 00:17:00,595 And so very hard to solve with explicit single step methods. 227 00:17:00,754 --> 00:17:04,054 And now we try a few collocation single step methods. 228 00:17:04,213 --> 00:17:11,573 Let's first use equidistant collocation points, the most natural choice if you have no other 229 00:17:11,732 --> 00:17:12,972 clue. 230 00:17:13,131 --> 00:17:14,691 Equidistant collocation points. 231 00:17:14,850 --> 00:17:21,210 And here you see the error plot in doubly logarithmic scale, error versus uniform time 232 00:17:21,369 --> 00:17:22,369 step h. 233 00:17:22,528 --> 00:17:25,928 And the error points are very nicely aligned. 234 00:17:26,087 --> 00:17:30,667 And so algebraic convergence is fairly obvious. 235 00:17:30,826 --> 00:17:34,206 And we can even read off the rate when you determine the slope. 236 00:17:34,365 --> 00:17:40,045 I can tell you that we have the rate s for even s. 237 00:17:40,204 --> 00:17:42,164 s is the number of collocation points. 238 00:17:42,323 --> 00:17:44,723 And for odd s, even s plus 1. 239 00:17:44,882 --> 00:17:51,622 And of course, as you know, this rate is a good predictor for the order of the collocation 240 00:17:51,781 --> 00:17:52,881 single step method. 241 00:17:53,040 --> 00:17:59,960 So they work rather well, at least for these equidistant collocation points. 242 00:18:00,119 --> 00:18:03,239 And the main message is we can really achieve higher order. 243 00:18:03,398 --> 00:18:11,818 I should have added, why is there hope that these methods enjoy unconditional stability? 244 00:18:11,977 --> 00:18:15,117 And I would like to ask a question to you. 245 00:18:15,276 --> 00:18:19,836 I claim that the implicit Euler method is a collocation method. 246 00:18:19,995 --> 00:18:22,515 Please find out the set of collocation points. 247 00:18:22,674 --> 00:18:25,094 So equidistant collocation points work. 248 00:18:25,253 --> 00:18:30,093 But now let's look at these single step methods from a slightly different perspective because 249 00:18:30,252 --> 00:18:34,832 this will give us another hint on how to choose collocation points. 250 00:18:34,991 --> 00:18:39,471 Let's link collocation single step methods and quadrature. 251 00:18:39,630 --> 00:18:45,190 So we already know that numerical integration is, of course, closely related to numerical 252 00:18:45,349 --> 00:18:49,349 quadrature because when you solve such an initial value problem, scalar initial value 253 00:18:49,508 --> 00:18:54,488 problem, approximately, then this amounts to computing the integral of the right hand 254 00:18:54,647 --> 00:18:55,867 side function f. 255 00:18:56,026 --> 00:18:59,186 And now, well, let's briefly switch back. 256 00:18:59,345 --> 00:19:02,185 We had this formula for the next state. 257 00:19:02,344 --> 00:19:06,144 And now let's apply this formula for an f that depends only on t. 258 00:19:06,303 --> 00:19:12,743 Then, of course, these kj vectors are given straightforwardly. 259 00:19:12,902 --> 00:19:16,502 We don't need to solve any nonlinear system. 260 00:19:16,661 --> 00:19:20,861 They're just f of t naught plus cjh. 261 00:19:21,020 --> 00:19:27,620 And so what we find is that for this really simple case, well, maybe you should try by 262 00:19:27,779 --> 00:19:28,779 yourself now. 263 00:19:28,938 --> 00:19:32,898 Maybe pause the video briefly and find the formula for y1. 264 00:19:33,057 --> 00:19:44,297 Well, I told you the kj vectors are just, say, the f evaluated t naught plus cjh, j 265 00:19:44,456 --> 00:19:46,536 from 1 to s. 266 00:19:46,695 --> 00:19:51,215 And then we integrate the Lagrange polynomials. 267 00:19:51,374 --> 00:19:54,934 And this now is a classical quadrature formula. 268 00:19:55,093 --> 00:19:59,333 This amounts to an s point quadrature formula. 269 00:19:59,492 --> 00:20:07,012 Of course, adapted to the interval here at 0h, when you choose h equals 1, then the nodes 270 00:20:07,171 --> 00:20:08,731 would be cj. 271 00:20:08,890 --> 00:20:16,850 So the cj parameters, they apparently play the role of quadrature nodes in this example. 272 00:20:17,009 --> 00:20:25,809 So we may guess that it's a good idea to choose them from suitable sets of quadrature nodes. 273 00:20:25,968 --> 00:20:30,288 And now remember the chapter on numerical quadrature. 274 00:20:30,447 --> 00:20:33,527 There we found ideal quadrature nodes. 275 00:20:33,686 --> 00:20:40,686 And these ideal quadrature nodes that gave us quadrature rules with maximal order were 276 00:20:40,845 --> 00:20:42,085 the Gauss points. 277 00:20:42,244 --> 00:20:49,484 OK, so let's try the Gauss points, of course, transformed to the interval 0, 1 as collocation 278 00:20:49,643 --> 00:20:50,963 point cj. 279 00:20:51,122 --> 00:20:55,942 And this is what I did also for the logistic ODE in a numerical experiment. 280 00:20:56,101 --> 00:21:02,441 And now we observe very nice algebraic convergence again, but with an amazing order. 281 00:21:02,600 --> 00:21:08,800 Because now it seems when we have the s point collocation single step method that we really 282 00:21:08,959 --> 00:21:11,799 observe an order of 2s. 283 00:21:11,958 --> 00:21:17,718 And this is not a coincidence because collocation single step methods have been investigated 284 00:21:17,877 --> 00:21:18,897 thoroughly. 285 00:21:19,056 --> 00:21:22,596 They admit a rather comprehensive theory. 286 00:21:22,755 --> 00:21:28,515 And the central result of that theory is the following theorem, namely that the order of 287 00:21:28,674 --> 00:21:34,514 a collocation single step method agrees with the order of the quadrature rule defined by 288 00:21:34,673 --> 00:21:36,313 the collocation points. 289 00:21:36,472 --> 00:21:37,472 They're the same. 290 00:21:37,631 --> 00:21:42,191 And we know that s point Gauss quadrature has order 2s. 291 00:21:42,350 --> 00:21:46,750 And this is inherited by the collocation single step method. 292 00:21:46,909 --> 00:21:52,649 Now with the collocation single step methods we have seen, let's call it an arbitrary order 293 00:21:52,808 --> 00:21:55,228 generalization of the implicit Euler method. 294 00:21:55,387 --> 00:22:02,067 Now it turns out that these collocation single step methods very nicely fit a slight generalization 295 00:22:02,226 --> 00:22:04,186 of Runge Kutta single step method. 296 00:22:04,345 --> 00:22:07,145 Let's recall their defining equations. 297 00:22:07,304 --> 00:22:12,184 First for the collocation single step methods based on the collocation point cj. 298 00:22:12,343 --> 00:22:16,623 These are the equations we got for the kj vectors. 299 00:22:16,782 --> 00:22:20,942 These implicit definition through a nonlinear system of equations. 300 00:22:21,101 --> 00:22:24,981 And then the equation for the next state. 301 00:22:25,140 --> 00:22:29,020 Now I've introduced coefficients and defined them on the right. 302 00:22:29,179 --> 00:22:34,419 Deliberately say of course, because now we would like to compare this with the familiar 303 00:22:34,578 --> 00:22:38,218 formulas for explicit Runge Kutta single step method. 304 00:22:38,377 --> 00:22:41,897 And now you clearly see the correspondence. 305 00:22:42,056 --> 00:22:47,936 The implicit equations for the kis are very similar to the increment equations. 306 00:22:48,095 --> 00:22:53,815 And the formula for the next state clearly reminds you of this update part of the Runge 307 00:22:53,974 --> 00:22:55,394 Kutta single step method. 308 00:22:55,553 --> 00:23:00,693 And now you see why I've named the abbreviations aij and bi. 309 00:23:00,852 --> 00:23:05,812 Because now we see that the formulas for the collocation single step method are almost 310 00:23:05,971 --> 00:23:13,651 the same as for Runge Kutta single step method, explicit ones, with one slight difference. 311 00:23:13,810 --> 00:23:17,610 And that slight difference is the summation bound. 312 00:23:17,769 --> 00:23:22,569 For explicit Runge Kutta single step methods we sum up to i minus one here in the increment 313 00:23:22,728 --> 00:23:24,968 equation and here we sum up to s. 314 00:23:25,127 --> 00:23:26,367 This is the only difference. 315 00:23:26,526 --> 00:23:32,846 And now it's very natural to generalize the definition of general Runge Kutta single step 316 00:23:33,005 --> 00:23:33,985 methods. 317 00:23:34,144 --> 00:23:40,364 Just replace i minus one with s and this gives us now the full family of Runge Kutta single 318 00:23:40,523 --> 00:23:41,523 step method. 319 00:23:41,682 --> 00:23:50,962 Now notice that the increment equations boil down to an ns times ns nonlinear systems of 320 00:23:51,121 --> 00:23:52,921 equations. 321 00:23:53,080 --> 00:23:55,720 This accounts for the implicit nature. 322 00:23:55,879 --> 00:24:02,619 Otherwise most aspects of that we have found with explicit Runge Kutta single step method 323 00:24:02,778 --> 00:24:05,178 carry over to this larger family. 324 00:24:05,337 --> 00:24:12,917 In particular again they will enjoy algebraic convergence of an integer order with an integer 325 00:24:13,076 --> 00:24:15,516 rate which also defines their order. 326 00:24:15,675 --> 00:24:17,475 There's nothing new in this respect. 327 00:24:17,634 --> 00:24:19,914 Convergence theory is almost the same. 328 00:24:20,073 --> 00:24:23,473 And what also carries over is the notation. 329 00:24:23,632 --> 00:24:29,372 Of course the general Runge Kutta single step method will usually be written down in a butcher 330 00:24:29,531 --> 00:24:30,591 scheme too. 331 00:24:30,750 --> 00:24:36,270 Now what is new now is that the butcher matrix A in the butcher scheme will now be a general 332 00:24:36,429 --> 00:24:37,989 dense matrix. 333 00:24:38,148 --> 00:24:44,488 No longer only a strictly triangular lower triangular matrix as in the explicit case. 334 00:24:44,647 --> 00:24:51,487 Now we still have to find out whether this full class of general Runge Kutta single step 335 00:24:51,646 --> 00:24:59,286 method really contains methods that succeed in tackling stiff initial value problems efficiently. 336 00:24:59,445 --> 00:25:05,485 But before I'm doing this I would like to address a few sort of implementation aspects. 337 00:25:05,644 --> 00:25:12,344 So the first is of course a crucial one namely how to solve the nonlinear increment equation. 338 00:25:12,503 --> 00:25:18,163 And in this context first I would like to introduce the so called stage form. 339 00:25:18,322 --> 00:25:21,922 Well for the simpler case of an autonomous ODE. 340 00:25:22,081 --> 00:25:23,221 What is the stage form? 341 00:25:23,380 --> 00:25:31,000 The stage form relies on introducing new auxiliary vectors usually called GI's. 342 00:25:31,159 --> 00:25:39,939 And these GI vectors are just an abbreviation for arguments to the f function in the increment 343 00:25:40,098 --> 00:25:41,098 equations. 344 00:25:41,257 --> 00:25:43,477 So what you see here this is the GI. 345 00:25:43,636 --> 00:25:49,036 So of course they are closely related to the increments because once you know the GI's 346 00:25:49,195 --> 00:25:53,815 you can recover the KI's by f of y naught plus GI. 347 00:25:53,974 --> 00:26:01,234 So these GI's are called the stages in contrast to the KI's that are called the increments 348 00:26:01,393 --> 00:26:02,593 as we have seen. 349 00:26:02,752 --> 00:26:07,472 Now we can rewrite the increment equations in a different form and then they are called 350 00:26:07,631 --> 00:26:09,431 stage equations fittingly. 351 00:26:09,590 --> 00:26:17,730 So from the increment equations that are given in the definition we obtain the stage equations 352 00:26:17,889 --> 00:26:24,489 by well we write down more or less the definition of the stage and plug in this should be s 353 00:26:24,648 --> 00:26:32,308 here and plug in instead of the increment f y naught plus g j. 354 00:26:32,467 --> 00:26:38,387 Ok now this is an ns times nf nonlinear system of equations for the stages. 355 00:26:38,546 --> 00:26:44,506 And then of course also the update equation can be rewritten in terms of the stage vectors. 356 00:26:44,665 --> 00:26:47,705 Why is this stage form interesting? 357 00:26:47,864 --> 00:26:53,104 It has to do with the following observation namely if you use a small time step which 358 00:26:53,263 --> 00:26:58,903 may be necessary for the sake of accuracy then since the time step multiplies the sum 359 00:26:59,062 --> 00:27:05,542 in the definition of the j stage the stages will probably be very small in terms of norm. 360 00:27:05,701 --> 00:27:13,661 So and as a consequence we have that zero is probably a good initial guess for an iterative 361 00:27:13,820 --> 00:27:16,980 method for solving the stage equation. 362 00:27:17,139 --> 00:27:22,019 For instance in the context of Newton's method and this is what I am going to discuss now. 363 00:27:22,178 --> 00:27:26,058 How can we use Newton's method to solve the stage equations? 364 00:27:26,217 --> 00:27:31,457 Remember the benefit of the stage equations is that the initial guess is very clear. 365 00:27:31,616 --> 00:27:37,236 So well these are the stage equations that we want to solve for the g i vectors. 366 00:27:37,395 --> 00:27:44,315 Now in order to have a compact notation we stack all the g i's vectors on top of each 367 00:27:44,474 --> 00:27:50,054 other and call this n times s vector g arrow. 368 00:27:50,213 --> 00:27:57,013 Now let's rewrite the stage equations as a zero point equation for a function capital 369 00:27:57,172 --> 00:28:01,652 F of course to set the stage for formulating the Newton method. 370 00:28:01,811 --> 00:28:10,131 Ok the stage equation so this is the left hand side and then on the right we also would 371 00:28:10,290 --> 00:28:13,450 like to rewrite this in a compact form. 372 00:28:13,609 --> 00:28:20,329 Now we introduce the Butcher matrix and take its Kronecker product with the n by n identity 373 00:28:20,488 --> 00:28:21,488 matrix. 374 00:28:21,647 --> 00:28:25,767 So this is the Kronecker product that we had in the beginning of the course. 375 00:28:25,926 --> 00:28:32,366 Kronecker product of matrices and then we multiply this again with a big vector f of 376 00:28:32,525 --> 00:28:41,285 y naught plus g 1 and this is again an n times s vector because Kronecker product is an ns 377 00:28:41,444 --> 00:28:44,984 times ns matrix and this should be zero. 378 00:28:45,143 --> 00:28:47,163 So this is the capital F. 379 00:28:47,322 --> 00:28:53,662 And now what is commonly done for solving this is not applying the full Newton method 380 00:28:53,821 --> 00:28:55,721 but the simplified Newton method. 381 00:28:55,880 --> 00:29:01,700 There was a remark in the chapter on solving nonlinear systems on the simplified Newton 382 00:29:01,859 --> 00:29:02,879 method. 383 00:29:03,038 --> 00:29:06,238 You may now pause the video and look up that remark. 384 00:29:06,397 --> 00:29:11,817 In a nutshell the simplified method, Newton method is obtained from the regular Newton 385 00:29:11,976 --> 00:29:14,316 method by freezing the Jacobian. 386 00:29:14,475 --> 00:29:21,435 So freezing the Jacobian means that you use the Jacobian at initial guess all the time 387 00:29:21,594 --> 00:29:22,794 for all the steps. 388 00:29:22,953 --> 00:29:29,753 And so this is the recursion for the simplified Newton method for say f of g vector equals 389 00:29:29,912 --> 00:29:32,212 zero with initial guess zero. 390 00:29:32,371 --> 00:29:36,231 And the initial guess zero of course we have just motivated here. 391 00:29:36,390 --> 00:29:39,770 Now let's compute that Jacobian. 392 00:29:39,929 --> 00:29:47,309 So dF zero you have to differentiate with respect to the g vector. 393 00:29:47,468 --> 00:29:48,468 Here this is easy. 394 00:29:48,627 --> 00:29:50,927 This is just the identity. 395 00:29:51,086 --> 00:29:54,746 So the Jacobian is a big identity matrix. 396 00:29:54,905 --> 00:30:00,585 Multiplication with a matrix is well can just be pulled out and then we have the Fs and 397 00:30:00,744 --> 00:30:04,024 here the components of this vector we have to differentiate. 398 00:30:04,183 --> 00:30:12,043 Then we obtain the dF at y zero but of course thanks to the initial guess zero that's all. 399 00:30:12,202 --> 00:30:19,902 So we get the Kronecker product of the Butcher matrix and the Jacobian of f at y zero. 400 00:30:20,061 --> 00:30:21,741 I admit this is fast. 401 00:30:21,900 --> 00:30:28,700 You may also pause the video now and elaborate the details say maybe for the scalar case 402 00:30:28,859 --> 00:30:32,739 to understand how this Jacobian emerges. 403 00:30:32,898 --> 00:30:36,578 This is of course an ns times ns matrix. 404 00:30:36,737 --> 00:30:44,317 And now we see that this Jacobian is equal to the identity matrix for h equals zero. 405 00:30:44,476 --> 00:30:48,136 So if h is small then it will certainly be invertible. 406 00:30:48,295 --> 00:30:53,735 So there's no problem in carrying out this simplified Newton iteration. 407 00:30:53,894 --> 00:30:57,734 All right so this is about implementation. 408 00:30:57,893 --> 00:31:04,033 But now the urgent question is are these methods suitable or at least some of them, some of 409 00:31:04,192 --> 00:31:09,212 these new implicit Runge Kutta methods, are they suitable for stiff problems? 410 00:31:09,371 --> 00:31:13,571 To understand this well we start with model problem analysis. 411 00:31:13,730 --> 00:31:19,250 And as you know from the previous section is once we understand the behavior for linear 412 00:31:19,409 --> 00:31:25,649 model problems then thanks to the linearization idea well we can also draw conclusions about 413 00:31:25,808 --> 00:31:28,608 the behavior for general nonlinear ODs. 414 00:31:28,767 --> 00:31:35,467 So linear model problem, scalar case, apply say that a general Runge Kutta single step 415 00:31:35,626 --> 00:31:36,626 method. 416 00:31:36,785 --> 00:31:41,545 And then what we will find out is that as in the case of explicit Runge Kutta single 417 00:31:41,704 --> 00:31:48,884 step method with exactly the same manipulations and arguments we find that the resulting recursion 418 00:31:49,043 --> 00:31:56,323 is defined by a multiplication with a number and that number is s of h lambda where s is 419 00:31:56,482 --> 00:32:00,442 a particular function, the stability function and also in this case. 420 00:32:00,601 --> 00:32:06,941 And the stability function can be computed as in the explicit case, almost as in the 421 00:32:07,100 --> 00:32:08,100 explicit case. 422 00:32:08,259 --> 00:32:10,319 And this is why I'm skipping the details here. 423 00:32:10,478 --> 00:32:12,298 I'm showing you the final result. 424 00:32:12,457 --> 00:32:16,257 And the final result gives concrete formulas for the stability function. 425 00:32:16,416 --> 00:32:22,956 So by either Bloch Gauss elimination you obtain the first formula, this one, or by Cramers 426 00:32:23,115 --> 00:32:28,875 rule, let me give you the keyword, Cramers rule for linear systems, you get this second 427 00:32:29,034 --> 00:32:31,614 formula involving determinants. 428 00:32:31,773 --> 00:32:36,953 And this second formula is rather illuminating. 429 00:32:37,112 --> 00:32:42,652 Remember that the arguments of the determinants are matrices that are linear functions of 430 00:32:42,811 --> 00:32:45,231 the set argument of the stability function. 431 00:32:45,390 --> 00:32:52,670 And remember that the determinant of an S by S matrix is a sum of S size products of 432 00:32:52,829 --> 00:32:53,869 its entries. 433 00:32:54,028 --> 00:32:57,828 Here of S monomials in set. 434 00:32:57,987 --> 00:33:06,047 And this means that both determinants are polynomials of degree at most S in set. 435 00:33:06,206 --> 00:33:09,746 But now both are generic polynomials. 436 00:33:09,905 --> 00:33:17,345 So here now the stability function is a quotient of polynomials and such a function is called 437 00:33:17,504 --> 00:33:21,024 a rational function, a term that you may have heard in analysis. 438 00:33:21,183 --> 00:33:23,423 I hope you've heard in analysis. 439 00:33:23,582 --> 00:33:29,462 More precisely, we have that the stability function now is a quotient of two polynomials 440 00:33:29,621 --> 00:33:37,261 of degree at most S. So P and Q call it numerator denominator polynomials of degree S. So instead 441 00:33:37,420 --> 00:33:42,340 of plain polynomials, now we get rational functions as stability functions. 442 00:33:42,499 --> 00:33:47,739 Let me give you a few examples of simple implicit Runge Kutta single step method. 443 00:33:47,898 --> 00:33:54,898 The two that we have already seen earlier, implicit Euler, 1 minus 1 over Z, essentially 444 00:33:55,057 --> 00:33:59,617 we already computed this function when we did linear model problem analysis for the 445 00:33:59,776 --> 00:34:01,576 decay equation. 446 00:34:01,735 --> 00:34:06,255 Implicit midpoint rule, well, you can just evaluate the formula and you get this rational 447 00:34:06,414 --> 00:34:07,594 function. 448 00:34:07,753 --> 00:34:13,853 Now of course, you know that this stability function tells us everything about the behavior 449 00:34:14,012 --> 00:34:18,732 of the single step method when applied to a linear ODE. 450 00:34:18,891 --> 00:34:22,811 In particular, it tells us about possible blowup. 451 00:34:22,970 --> 00:34:30,150 If lambda H is not in the region of stability, then there will be blowup when we apply the 452 00:34:30,309 --> 00:34:33,689 single step method to the scalar linear ODE. 453 00:34:33,848 --> 00:34:38,608 And the region of stability, remember, is defined as those complex numbers for which 454 00:34:38,767 --> 00:34:41,947 the modulus of the stability function is smaller than 1. 455 00:34:42,106 --> 00:34:46,726 And now let's look at the regions of stability for the two implicit Runge Kutta single step 456 00:34:46,885 --> 00:34:47,985 methods. 457 00:34:48,144 --> 00:34:49,764 Implicit Euler, implicit midpoint rule. 458 00:34:49,923 --> 00:34:55,463 I've computed them and here they are displayed, the green shaded regions. 459 00:34:55,622 --> 00:34:57,722 These are the regions of stability. 460 00:34:57,881 --> 00:35:01,341 And first of all, we observe that they are unbounded. 461 00:35:01,500 --> 00:35:04,720 So S of psi is now unbounded. 462 00:35:04,879 --> 00:35:10,899 And secondly, we observe that they contain the left half plane of the complex plane, 463 00:35:11,058 --> 00:35:12,278 this area. 464 00:35:12,437 --> 00:35:15,197 It's even exactly the left half plane. 465 00:35:15,356 --> 00:35:21,276 So now, if in the linear model problem, the real part of lambda is negative, this means 466 00:35:21,435 --> 00:35:30,435 that for any time step size, H lambda will be contained in the region of stability for 467 00:35:30,594 --> 00:35:32,254 any positive lambda. 468 00:35:32,413 --> 00:35:37,473 And this means, of course, that there will never be a blowup in the case of the decay 469 00:35:37,632 --> 00:35:38,632 equation. 470 00:35:38,791 --> 00:35:41,131 Well, this is what is unconditional stability. 471 00:35:41,290 --> 00:35:48,070 So this implicit Runge Kutta single step method have the potential for unconditional stability 472 00:35:48,229 --> 00:35:50,269 thanks to their rational stability function. 473 00:35:50,428 --> 00:35:55,368 Well, for this unconditional stability, it's essential that the left half plane belongs 474 00:35:55,527 --> 00:35:57,307 to the region of stability. 475 00:35:57,466 --> 00:36:01,966 And this property is so important that it has been given a separate name. 476 00:36:02,125 --> 00:36:03,925 It's called A stability. 477 00:36:04,084 --> 00:36:05,264 And here is the definition. 478 00:36:05,423 --> 00:36:10,923 The left half plane, abbreviated with C minus, if it's contained in the region of stability, 479 00:36:11,082 --> 00:36:15,782 then we call the single step method A stable from absolutely stable. 480 00:36:15,941 --> 00:36:22,501 Could also be called U stable, unconditionally stable, because, well, the property is that 481 00:36:22,660 --> 00:36:26,740 there's no time step constraint so that blowup is avoided. 482 00:36:26,899 --> 00:36:34,059 So now we can transfer this to general stiff problems and make the statement that these 483 00:36:34,218 --> 00:36:39,018 A stable methods, they don't have to struggle with stability induced time step constraint. 484 00:36:39,177 --> 00:36:40,397 They are none. 485 00:36:40,556 --> 00:36:44,776 So they are, I repeat this, unconditionally stable. 486 00:36:44,935 --> 00:36:49,615 Let's look at A stability and regions of stability a bit more closely. 487 00:36:49,774 --> 00:36:54,954 Let's find out what would be the ideal region of stability. 488 00:36:55,113 --> 00:37:01,413 And to understand what I mean by ideal, I'd like to ask you to remember that actually 489 00:37:01,572 --> 00:37:06,772 for a single step method, the stability function is an approximation of the exponential function, 490 00:37:06,931 --> 00:37:09,011 at least in a neighborhood of 0. 491 00:37:09,170 --> 00:37:13,890 So I would call the ideal stability function the exponential function. 492 00:37:14,049 --> 00:37:18,249 And now ideal region of stability, what can this be? 493 00:37:18,408 --> 00:37:22,708 It's the region of stability of the ideal stability function. 494 00:37:22,867 --> 00:37:29,767 So complex numbers for which the exponential function attains a modulus smaller one. 495 00:37:29,926 --> 00:37:34,126 Now question to you, what is this S psi star? 496 00:37:34,285 --> 00:37:36,405 Pause the video briefly and try to figure out. 497 00:37:36,564 --> 00:37:44,344 You know that X of a complex number now written in terms of real and imaginary part, the modulus 498 00:37:44,503 --> 00:37:46,523 is just e to the x. 499 00:37:46,682 --> 00:37:48,582 And when is this smaller than one? 500 00:37:48,741 --> 00:37:51,741 Of course, when X is negative. 501 00:37:51,900 --> 00:37:56,360 So we have that this ideal stability function is really C minus. 502 00:37:56,519 --> 00:38:00,919 The ideal stability function is, well, really ideal. 503 00:38:01,078 --> 00:38:07,278 Silly wording, but when you look at the consequences for the scalar linear model problem, so when 504 00:38:07,437 --> 00:38:11,317 the stability function of a single step method is the ideal one. 505 00:38:11,476 --> 00:38:15,596 And then for negative lambda, what can we conclude? 506 00:38:15,755 --> 00:38:17,115 And for positive lambda? 507 00:38:17,274 --> 00:38:24,914 And what will be the limit of both the exact solution and the sequence of approximate states 508 00:38:25,073 --> 00:38:30,153 as k or t tends to infinity, respectively? 509 00:38:30,312 --> 00:38:32,552 You see the pause sign, try to figure out. 510 00:38:32,711 --> 00:38:38,111 Well, this is not difficult when you remember the recursion and how it's connected to the 511 00:38:38,270 --> 00:38:39,270 stability function. 512 00:38:39,429 --> 00:38:43,869 The S lambda H is just the multiplying factor in each time step. 513 00:38:44,028 --> 00:38:50,028 Then of course we see if real part of lambda smaller than zero, thanks to the ideal stability, 514 00:38:50,187 --> 00:38:54,667 region of stability, this S H lambda has a modulus smaller than one. 515 00:38:54,826 --> 00:38:56,306 So this will tend to zero. 516 00:38:56,465 --> 00:39:00,685 And this is great because then also the solution will decay to zero. 517 00:39:00,844 --> 00:39:06,184 And for lambda positive, this modulus of H lambda will be larger equal to one. 518 00:39:06,343 --> 00:39:08,463 Let's assume it's larger one in most cases. 519 00:39:08,622 --> 00:39:12,142 So this will tend to infinity and this will also tend to infinity. 520 00:39:12,301 --> 00:39:18,061 So the numerical solution very well captures the qualitative decay and growth behavior 521 00:39:18,220 --> 00:39:21,540 of the exact solution. 522 00:39:21,699 --> 00:39:26,259 Are there implicit Rumi Kutta single step methods with ideal region of stability? 523 00:39:26,418 --> 00:39:28,418 Yes, of course. 524 00:39:28,577 --> 00:39:30,617 We see it on the right here. 525 00:39:30,776 --> 00:39:34,676 For the implicit midpoint rule, we have the ideal region of stability. 526 00:39:34,835 --> 00:39:40,535 So we already have a specimen that enjoys this remarkable property, implicit midpoint 527 00:39:40,694 --> 00:39:41,694 method. 528 00:39:41,853 --> 00:39:47,533 And this is only one of infinitely many because the implicit midpoint method is in fact the 529 00:39:47,692 --> 00:39:51,612 simplest Gauss collocation single step method. 530 00:39:51,771 --> 00:39:57,851 And now let's look at the stability functions of other Gauss collocation single step methods. 531 00:39:58,010 --> 00:40:00,050 Here I've plotted their level lines. 532 00:40:00,209 --> 00:40:02,369 Well, implicit midpoint rule. 533 00:40:02,528 --> 00:40:06,528 Well the level equal one agrees with the imaginary axis. 534 00:40:06,687 --> 00:40:09,967 But now we see the same for all others. 535 00:40:10,126 --> 00:40:14,726 Modulus of S is smaller one to the left and bigger one to the right. 536 00:40:14,885 --> 00:40:19,845 And this is a rigorous property of all Gauss collocation single step method. 537 00:40:20,004 --> 00:40:23,884 Their regions of stability are all ideal. 538 00:40:24,043 --> 00:40:29,803 In particular, they are all A stable because A stability is a consequence of an ideal region 539 00:40:29,962 --> 00:40:31,042 of stability. 540 00:40:31,201 --> 00:40:32,201 Impressive. 541 00:40:32,360 --> 00:40:33,360 Problem solved. 542 00:40:33,519 --> 00:40:36,759 With Gauss collocation methods, they can deal with everything. 543 00:40:36,918 --> 00:40:40,918 Unfortunately, this is not the full truth. 544 00:40:41,077 --> 00:40:42,517 Let's look at an example. 545 00:40:42,676 --> 00:40:49,116 So this example studies a really stiff scalar initial value problem. 546 00:40:49,275 --> 00:40:52,635 You see here, the decay rate is a million. 547 00:40:52,794 --> 00:40:55,234 This is a non autonomous case. 548 00:40:55,393 --> 00:41:04,313 So in this case, the solution will almost instantaneously approach this sine curve. 549 00:41:04,472 --> 00:41:08,912 So the sine solution is really strongly attractive. 550 00:41:09,071 --> 00:41:13,191 So this is what should be produced by a stable single step method. 551 00:41:13,350 --> 00:41:15,750 Now let's look at numerical solutions. 552 00:41:15,909 --> 00:41:22,029 The dashed green curve, let me highlight it here, is the exact solution. 553 00:41:22,188 --> 00:41:24,188 It goes here. 554 00:41:24,347 --> 00:41:25,867 This is the sine curve. 555 00:41:26,026 --> 00:41:31,306 Implicit Euler red pluses, they exactly lie on the sine curve. 556 00:41:31,465 --> 00:41:36,505 So implicit Euler, a first order method, produces a perfect solution. 557 00:41:36,664 --> 00:41:41,064 So implicit Euler is essentially on the exact solution. 558 00:41:41,223 --> 00:41:48,223 Unfortunately, the much wanted Gauss collocation single step method, they do a really poor 559 00:41:48,382 --> 00:41:49,382 job. 560 00:41:49,541 --> 00:41:50,541 You see, they are really off. 561 00:41:50,700 --> 00:41:55,220 They either oscillate or just by a factor of two off. 562 00:41:55,379 --> 00:41:56,379 How come? 563 00:41:56,538 --> 00:42:00,618 Because these methods are of order two, four, six, eight. 564 00:42:00,777 --> 00:42:03,337 In order eight cannot fix the problem. 565 00:42:03,496 --> 00:42:05,196 So what's going wrong? 566 00:42:05,355 --> 00:42:07,175 Let's look at the stability functions. 567 00:42:07,334 --> 00:42:11,614 For implicit Euler, we know that the stability function, well, it's like this. 568 00:42:11,773 --> 00:42:13,133 It's the case two zero. 569 00:42:13,292 --> 00:42:17,452 It's one over one minus Z for Z to minus infinity. 570 00:42:17,611 --> 00:42:24,371 For Gauss collocation methods, negative arguments, they stay below one or above minus one. 571 00:42:24,530 --> 00:42:25,530 Perfect. 572 00:42:25,689 --> 00:42:28,369 Ideal region of stability, of course. 573 00:42:28,528 --> 00:42:35,388 But what we see is that the modulus of the stability function tends to one as the argument 574 00:42:35,547 --> 00:42:37,207 tends to minus infinity. 575 00:42:37,366 --> 00:42:43,806 And this is not the right thing because you know that the ideal stability function is 576 00:42:43,965 --> 00:42:48,245 the exponential function which tends to zero in the same limit. 577 00:42:48,404 --> 00:42:55,204 Now, when the stability function is very close to one, remember stability function is evaluated 578 00:42:55,363 --> 00:43:01,843 for H lambda is the multiplying factor that defines the time stepping recursion. 579 00:43:02,002 --> 00:43:04,202 Then we have hardly any damping. 580 00:43:04,361 --> 00:43:07,201 And this is what we see in the experiment above. 581 00:43:07,360 --> 00:43:13,400 Because what we see is that instead of instantly dropping to this sine curve, well, the numerical 582 00:43:13,559 --> 00:43:19,839 solution just lingers around somewhere else and doesn't show a tendency of tending to 583 00:43:19,998 --> 00:43:20,998 the sine curve. 584 00:43:21,157 --> 00:43:28,597 There's too little damping in the numerical method to capture this very fast decay in 585 00:43:28,756 --> 00:43:29,756 this model. 586 00:43:29,915 --> 00:43:34,915 Implicit Euler, in contrast, does a great job because there we have a stability function, 587 00:43:35,074 --> 00:43:41,954 as you know, one over one minus Z, which clearly tends to zero for Z to infinity, minus infinity. 588 00:43:42,113 --> 00:43:48,093 And this is the right behavior because it reflects the same behavior for the exponential 589 00:43:48,252 --> 00:43:49,252 function. 590 00:43:49,411 --> 00:43:53,951 So this stability or ideal domains of stability is not enough. 591 00:43:54,110 --> 00:43:56,150 We need some stronger stability. 592 00:43:56,309 --> 00:43:58,709 And this is known as L stability. 593 00:43:58,868 --> 00:44:01,508 I can't tell you where this letter L comes from. 594 00:44:01,667 --> 00:44:04,107 This is the name used in computational math. 595 00:44:04,266 --> 00:44:05,366 What is L stability? 596 00:44:05,525 --> 00:44:09,605 It's not only A stability, which is the first requirement. 597 00:44:09,764 --> 00:44:15,144 It's also the requirement that you have this zero limit for the stability function as the 598 00:44:15,303 --> 00:44:19,483 argument tends to minus infinity or the real part tends to minus infinity. 599 00:44:19,642 --> 00:44:20,882 This is L stability. 600 00:44:21,041 --> 00:44:29,621 So you could say an L stable method is A stable and the stability function at minus infinity 601 00:44:29,780 --> 00:44:31,060 is zero. 602 00:44:31,219 --> 00:44:34,159 Sloppily speaking, I'm not writing here. 603 00:44:34,318 --> 00:44:37,218 Now let's see when a method is L stable. 604 00:44:37,377 --> 00:44:43,017 This can be clearly pinpointed by a particular property of the Butcher scheme. 605 00:44:43,176 --> 00:44:48,256 So we know the formula for the stability function from this theorem 344. 606 00:44:48,415 --> 00:44:53,135 And now let's rewrite this formula slightly by pulling out the set. 607 00:44:53,294 --> 00:44:59,854 And we have a B transpose and then we have a one over set minus A, Butcher matrix inverse, 608 00:45:00,013 --> 00:45:01,853 and then we have this vector of all ones. 609 00:45:02,012 --> 00:45:08,012 Okay, the one over set is convenient because when you plug in, say, infinity, then it's 610 00:45:08,171 --> 00:45:09,171 zero. 611 00:45:09,330 --> 00:45:15,030 So from this formula, we read off that the stability function at infinity will attain 612 00:45:15,189 --> 00:45:20,689 the value B transpose, Butcher matrix inverse multiplied with all one vector. 613 00:45:20,848 --> 00:45:23,928 And now this should be zero for L stability. 614 00:45:24,087 --> 00:45:25,327 This is what we want. 615 00:45:25,486 --> 00:45:28,206 So now when will this be zero? 616 00:45:28,365 --> 00:45:37,325 It will be zero if, say, this vector, Butcher matrix inverse times one multiplied with B 617 00:45:37,484 --> 00:45:39,764 transpose is equal to one. 618 00:45:39,923 --> 00:45:47,443 And this will be satisfied if one, say, rho of A is equal to B. 619 00:45:47,602 --> 00:45:49,282 So let me write it. 620 00:45:49,441 --> 00:45:53,521 Well, this is a notation for a single rho of A. 621 00:45:53,680 --> 00:46:00,960 If a single rho of A is equal to B, then of course B transpose A inverse means multiplying 622 00:46:01,119 --> 00:46:03,039 that rho with A inverse. 623 00:46:03,198 --> 00:46:11,718 Then as a consequence, we have that B transpose A inverse is a Cartesian coordinate vector, 624 00:46:11,877 --> 00:46:18,017 and then its inner product, a Cartesian coordinate rho vector, then its inner product with the 625 00:46:18,176 --> 00:46:20,016 all one vector will be one. 626 00:46:20,175 --> 00:46:26,435 And so this is a sufficient condition for this L stability requirement. 627 00:46:26,594 --> 00:46:34,034 And now from the Butcher scheme, we can see this algebraic condition because all we need 628 00:46:34,193 --> 00:46:38,753 to check is whether one rho of the Butcher matrix agrees with the B vector. 629 00:46:38,912 --> 00:46:41,932 Well, usually it will be the bottom row. 630 00:46:42,091 --> 00:46:47,111 So L stable implicit Runge Kutta method are usually written like this. 631 00:46:47,270 --> 00:46:53,150 So the bottom row of the Butcher matrix is equal to the weight vector B. 632 00:46:53,309 --> 00:46:58,869 Fortunately, the methods that look like this also enjoy a decent order. 633 00:46:59,028 --> 00:47:00,568 Let me give you examples. 634 00:47:00,727 --> 00:47:05,347 And the prime examples are the so called Gauss Radau methods. 635 00:47:05,506 --> 00:47:11,106 They are collocations single step methods where the collocation point Cs for s stages 636 00:47:11,265 --> 00:47:12,345 is equal to one. 637 00:47:12,504 --> 00:47:17,304 Then we end up with a Butcher scheme of the desired form. 638 00:47:17,463 --> 00:47:22,263 And when we then maximize the order, we arrive at the Gauss Radau methods. 639 00:47:22,422 --> 00:47:26,302 So these are the three simplest representatives. 640 00:47:26,461 --> 00:47:28,821 Simplest is the implicit Euler method. 641 00:47:28,980 --> 00:47:31,940 We already seen it's clearly L stable. 642 00:47:32,099 --> 00:47:37,259 And then the next one has order three and then the three stage method already has order 643 00:47:37,418 --> 00:47:38,938 five. 644 00:47:39,097 --> 00:47:47,137 And in general, this family can be constructed up to arbitrary order and the s stage representative 645 00:47:47,296 --> 00:47:49,216 will have order two s minus one. 646 00:47:49,375 --> 00:47:57,835 And now we can deal with initial value problem with extreme stiffness with a very rapid built 647 00:47:57,994 --> 00:47:58,994 in decay. 648 00:47:59,153 --> 00:48:01,133 Let's look at this finally examples. 649 00:48:01,292 --> 00:48:05,492 It's again the arrested blow up equation. 650 00:48:05,651 --> 00:48:13,011 Lambda equal 500 means that the fixed point y equal one is extremely strongly attracting. 651 00:48:13,170 --> 00:48:18,510 So all solutions very rapidly approach the state y equal one. 652 00:48:18,669 --> 00:48:23,369 So it has a very strong built in decay towards y equal one. 653 00:48:23,528 --> 00:48:28,648 We investigate single step methods with uniform time steps and these are the results. 654 00:48:28,807 --> 00:48:35,287 On the left Gauss collocation, what we see is that here we have extreme overshoots when 655 00:48:35,446 --> 00:48:41,146 the solution approaches the stable state because there's just too little damping. 656 00:48:41,305 --> 00:48:47,685 And so the solution either oscillates for a while or they are a bit off. 657 00:48:47,844 --> 00:48:54,484 In contrast, so here the Gauss Radau single step methods, no problem, only tiny overshoot, 658 00:48:54,643 --> 00:48:57,483 but not a significant problem. 659 00:48:57,642 --> 00:49:02,482 Excellent performance, good accuracy, and of course, unconditional stability because 660 00:49:02,641 --> 00:49:08,601 also for Gauss Radau single step methods, this is a theoretical result, is that they 661 00:49:08,760 --> 00:49:12,440 are all L stable, in particular, A stable. 662 00:49:12,599 --> 00:49:18,999 And with this numerical experiment and the good news that we can deal with even extremely 663 00:49:19,158 --> 00:49:24,418 stiff initial value problems using these Gauss Radau single step methods, I would like to 664 00:49:24,577 --> 00:49:25,857 close this section. 665 00:49:26,016 --> 00:49:31,556 And now we are equipped with methods to deal with almost every initial value problem. 666 00:49:31,715 --> 00:49:37,795 And in the remaining sections of this chapter, we'll see a few more special variants that 667 00:49:37,954 --> 00:49:41,734 can be used when the ODE has special structure.