1 00:00:00,140 --> 00:00:05,160 Welcome to the last part of chapter 11 of the course on numerical methods. 2 00:00:05,320 --> 00:00:10,120 We are studying numerical methods for the approximate solution of initial value problems 3 00:00:10,280 --> 00:00:15,000 for ordinary differential equations with a particular focus on single step methods. 4 00:00:15,160 --> 00:00:18,160 We have seen the important class of Runge Kutta single step methods. 5 00:00:18,320 --> 00:00:24,760 I have told you about theoretical properties in particular order of convergence. 6 00:00:24,920 --> 00:00:32,080 And now in this last part we are going to discuss a very important aspect of the practical 7 00:00:32,240 --> 00:00:34,320 use of single step methods. 8 00:00:34,480 --> 00:00:39,960 Because so far we have always taken for granted that a temporal mesh was available. 9 00:00:40,120 --> 00:00:46,580 However, in practice when you are given a single step method, you have no idea what 10 00:00:46,740 --> 00:00:50,100 could be a good mesh to be used for time stepping. 11 00:00:50,260 --> 00:00:53,080 You may opt for the simplest choice of uniform time step. 12 00:00:53,240 --> 00:00:56,920 But even then fixing that uniform time step is not clear. 13 00:00:57,080 --> 00:01:01,080 What time step will yield results with sufficient accuracy? 14 00:01:01,240 --> 00:01:08,020 All these issues are now addressed in the following video unit and they are essential 15 00:01:08,180 --> 00:01:13,400 for using single step methods in application contexts. 16 00:01:13,560 --> 00:01:17,600 The key word now is time step adaptation. 17 00:01:17,760 --> 00:01:22,860 In other words, we want an algorithm that is smart enough so that it can choose the 18 00:01:23,020 --> 00:01:28,720 temporal mesh, the time steps, well according to requirements specified by the user more 19 00:01:28,880 --> 00:01:34,100 precisely so that a user specified error tolerance is met. 20 00:01:34,260 --> 00:01:39,560 Why is it essential to have such algorithms for the solution of initial value problems? 21 00:01:39,720 --> 00:01:41,320 Let's start with a few examples. 22 00:01:41,480 --> 00:01:49,480 The first example is from computational reaction kinetics and studies a particular fairly complex 23 00:01:49,640 --> 00:01:50,640 chemical reaction. 24 00:01:50,800 --> 00:01:54,800 A famous one known as the Oregonator reaction. 25 00:01:54,960 --> 00:01:56,760 Well I'm not going into detail. 26 00:01:56,920 --> 00:02:03,240 Here you see the reaction scheme, the chemical compounds and from the reaction we can derive 27 00:02:03,400 --> 00:02:07,120 an ODE model describing the reaction kinetics. 28 00:02:07,280 --> 00:02:13,640 That is describing the temporal evolution of the concentrations of the involved compounds. 29 00:02:13,800 --> 00:02:17,320 This is the model with a particular reaction parameters. 30 00:02:17,480 --> 00:02:23,800 These reaction parameters can be fitted based on measurements and then when we solve these 31 00:02:23,960 --> 00:02:28,080 reaction equations, right now I'm not going to tell you how this is done, but let's assume 32 00:02:28,240 --> 00:02:29,720 that we can solve them. 33 00:02:29,880 --> 00:02:30,880 We see the following. 34 00:02:31,040 --> 00:02:38,560 So these are plots of the time dependent concentrations of two of the involved compounds and it's 35 00:02:38,720 --> 00:02:45,120 conspicuous that these concentrations vary heavily in time and in a highly non uniform 36 00:02:45,280 --> 00:02:45,780 way. 37 00:02:45,940 --> 00:02:51,660 Note its logarithmic scale, so the right plot shows that one of the concentrations almost 38 00:02:51,820 --> 00:02:57,360 stays constant for quite some time and then suddenly increases massively. 39 00:02:57,520 --> 00:03:04,260 There are these typical spikes and on the left you see a sort of exponential decrease 40 00:03:04,420 --> 00:03:10,820 of the concentrations over some time and then something violent happens and there's a sudden 41 00:03:10,980 --> 00:03:11,980 upshot. 42 00:03:12,140 --> 00:03:18,600 It goes without saying that to capture such a behavior one should use highly non uniform 43 00:03:18,760 --> 00:03:25,720 temporal meshes with mesh points clustering where things happen, where concentrations 44 00:03:25,880 --> 00:03:28,680 change rapidly over a short time. 45 00:03:28,840 --> 00:03:34,040 This is a typical example where adaptive step size control is inevitable. 46 00:03:34,200 --> 00:03:38,800 Another example that we have already seen is the explosion ODE. 47 00:03:38,960 --> 00:03:46,600 Very simple scalar ODE, y dot is y squared and some course units ago I showed you that 48 00:03:46,760 --> 00:03:50,720 this involves a blow up of the solutions. 49 00:03:50,880 --> 00:03:55,320 The solutions will only exist for some finite time and then they escape to infinity. 50 00:03:55,480 --> 00:04:02,400 So well, now imagine you want to do time stepping for the explosion ODE on an equidistant mesh 51 00:04:02,560 --> 00:04:06,160 using say the explicit Euler method. 52 00:04:06,320 --> 00:04:12,240 Well think about it, pause the video and answer the question what will happen if you use uniform 53 00:04:12,400 --> 00:04:15,520 H explicit Euler time stepping. 54 00:04:15,680 --> 00:04:22,560 Well it's clear when you apply explicit Euler to this ODE it will just go on forever. 55 00:04:22,720 --> 00:04:24,140 You will never see any blow up. 56 00:04:24,300 --> 00:04:30,440 It will of course produce states whose values are larger and larger but will never stop. 57 00:04:30,600 --> 00:04:33,980 So it fails to capture a really essential feature of the solution. 58 00:04:34,140 --> 00:04:40,680 Now let me give you a glimpse of what the algorithms to be presented in this section 59 00:04:40,840 --> 00:04:41,840 can accomplish. 60 00:04:42,000 --> 00:04:47,520 These are the states produced by this algorithm applied to the explosion ODE and what you 61 00:04:47,680 --> 00:04:53,760 see is this adaptive time stepping algorithm really manages to capture the blow up. 62 00:04:53,920 --> 00:05:00,600 So what it does is it reduces the time step more and more and more as it approaches the 63 00:05:00,760 --> 00:05:02,180 blow up time. 64 00:05:02,340 --> 00:05:07,880 And then well the computed states really seem to tend to infinity as functions of time. 65 00:05:08,040 --> 00:05:12,720 And this is exactly the behavior of the exact solution that is now reflected also by the 66 00:05:12,880 --> 00:05:14,180 numerical solution. 67 00:05:14,340 --> 00:05:15,340 Great job isn't it? 68 00:05:15,500 --> 00:05:18,540 So alright now let's see how this can be done. 69 00:05:18,700 --> 00:05:22,680 Let's see how we can achieve such an adaptive step size control. 70 00:05:22,840 --> 00:05:28,220 I would like to focus on methods that are local in time and these are the methods that 71 00:05:28,380 --> 00:05:30,920 are most widely used and most popular. 72 00:05:31,080 --> 00:05:39,440 So the goal is to find an algorithm that during time stepping builds the temporal mesh. 73 00:05:39,600 --> 00:05:44,660 Let's look at the guiding principles for the construction of such algorithms. 74 00:05:44,820 --> 00:05:52,900 In fact we have already seen this in the context of adaptive quadrature section 7.6 if we want 75 00:05:53,060 --> 00:05:54,200 to turn pages. 76 00:05:54,360 --> 00:05:58,380 The same ideas, I should say almost the same ideas were applied then. 77 00:05:58,540 --> 00:06:00,380 Now what are these ideas? 78 00:06:00,540 --> 00:06:06,880 The objective is twofold or I shouldn't call it, it's not a really realistic goal, let's 79 00:06:07,040 --> 00:06:08,400 call it a dream. 80 00:06:08,560 --> 00:06:14,000 We want to find a temporal mesh with as few as possible time points so that on that temporal 81 00:06:14,160 --> 00:06:19,120 mesh the error measure of interest is below a tolerance. 82 00:06:19,280 --> 00:06:23,560 And this tolerance should be specified by the user. 83 00:06:23,720 --> 00:06:28,900 So the algorithm should in principle be able to meet any prescribed tolerance. 84 00:06:29,060 --> 00:06:32,340 Of course taking into account the limits of machine precision. 85 00:06:32,500 --> 00:06:36,600 Why do we want the mesh to have as few points as possible? 86 00:06:36,760 --> 00:06:39,440 Of course we want our method to be efficient. 87 00:06:39,600 --> 00:06:45,040 The asymptotic theory about convergence of single step method of course tells us sending 88 00:06:45,200 --> 00:06:51,680 the mesh with 10 to 0, sending the number of mesh points to infinity, we can bring the 89 00:06:51,840 --> 00:06:54,280 error down below any threshold. 90 00:06:54,440 --> 00:06:58,400 But of course this could be at formidable, unacceptable cost. 91 00:06:58,560 --> 00:07:05,000 So mesh should be as coarse as possible and the solution should have a guaranteed quality. 92 00:07:05,160 --> 00:07:07,940 This is called the method should be reliable. 93 00:07:08,100 --> 00:07:11,440 Now what is the strategy for this? 94 00:07:11,600 --> 00:07:15,720 First of all the focus here is on methods that are local in time. 95 00:07:15,880 --> 00:07:16,880 So what does this mean? 96 00:07:17,040 --> 00:07:23,720 This means that when we are at a particular time step say tk, then the next point in time 97 00:07:23,880 --> 00:07:32,480 so tk plus 1 is chosen only based on time local information which is tk of course. 98 00:07:32,640 --> 00:07:34,720 Then of course we can use the state. 99 00:07:34,880 --> 00:07:42,760 This is what we know, current state and some time step size that we used previously. 100 00:07:42,920 --> 00:07:44,540 Why this locality in time? 101 00:07:44,700 --> 00:07:48,960 What we want to avoid is that the method tells us okay you're off, start again. 102 00:07:49,120 --> 00:07:54,640 Of course sometimes this would be maybe appropriate but it would not be satisfactory because it 103 00:07:54,800 --> 00:07:55,980 would be too costly. 104 00:07:56,140 --> 00:08:00,140 So this is cheap, usually with fixed cost per step. 105 00:08:00,300 --> 00:08:03,760 So locality in time is a constraint on our design. 106 00:08:03,920 --> 00:08:06,400 Then under this constraint what can we do? 107 00:08:06,560 --> 00:08:11,720 We have seen that the total error in a single step method well arises from the accumulation 108 00:08:11,880 --> 00:08:13,620 of one step errors. 109 00:08:13,780 --> 00:08:16,800 When we insist on locality in time we have little choice. 110 00:08:16,960 --> 00:08:19,600 All we can tackle are one step errors. 111 00:08:19,760 --> 00:08:21,800 Remember what is a one step error? 112 00:08:21,960 --> 00:08:28,040 It's the discrete evolution applied to the current state minus exact evolution applied 113 00:08:28,200 --> 00:08:29,560 to the current state. 114 00:08:29,720 --> 00:08:32,720 Here I'm writing it for autonomous ODE. 115 00:08:32,880 --> 00:08:36,840 So the policy is to control the one step error. 116 00:08:37,000 --> 00:08:43,040 So make the current time step small enough so this one step error is below a certain 117 00:08:43,200 --> 00:08:44,200 threshold. 118 00:08:44,360 --> 00:08:49,440 And then also predicting the next time step so that the one step error in the next step 119 00:08:49,600 --> 00:08:52,600 will probably be below our error threshold. 120 00:08:52,760 --> 00:08:59,520 The trouble is this one step error is not accessible because there is the exact evolution 121 00:08:59,680 --> 00:09:01,760 operator involved. 122 00:09:01,920 --> 00:09:06,280 So all we can do is to estimate this one step error. 123 00:09:06,440 --> 00:09:08,760 So what we need is another tool. 124 00:09:08,920 --> 00:09:12,720 We need a locality in time one step error estimator. 125 00:09:12,880 --> 00:09:18,240 Well before I'm going to discuss this let's briefly recall the one step error. 126 00:09:18,400 --> 00:09:21,560 Remember this was the picture for explicit Euler. 127 00:09:21,720 --> 00:09:27,920 It's of the one step where we start from the approximate state and compute the exact solution 128 00:09:28,080 --> 00:09:30,920 of a local initial value problem. 129 00:09:31,080 --> 00:09:40,160 Let's say y dot is f of y, so autonomous, and y of tk, sorry, y of tk is yk, so the 130 00:09:40,320 --> 00:09:41,640 approximate state. 131 00:09:41,800 --> 00:09:47,400 And this gives us the one step error if you indicate it by this orange arrow. 132 00:09:47,560 --> 00:09:49,440 Now how to estimate this? 133 00:09:49,600 --> 00:09:54,200 Again we rely on a strategy already employed in the context of numerical quadrature. 134 00:09:54,360 --> 00:10:00,860 It's the old idea of comparing two methods from which we can expect vastly different 135 00:10:01,020 --> 00:10:02,640 levels of accuracy. 136 00:10:02,800 --> 00:10:07,760 Here in the single step method context, methods of different order. 137 00:10:07,920 --> 00:10:15,200 So we assume that we've got two single step methods described by two discrete evolution 138 00:10:15,360 --> 00:10:18,640 operators which have different orders. 139 00:10:18,800 --> 00:10:22,760 And then we apply a heuristic consideration. 140 00:10:22,920 --> 00:10:28,780 We hope that the higher order method is much more accurate than the lower order method. 141 00:10:28,940 --> 00:10:36,440 So the higher order, well, approximation of the one step can be used as a Bonafide replacement 142 00:10:36,600 --> 00:10:39,620 for the action of the exact evolution operator. 143 00:10:39,780 --> 00:10:45,320 So what we do in the definition of the one step error, we replace the exact evolution 144 00:10:45,480 --> 00:10:47,840 by the higher order discrete evolution. 145 00:10:48,000 --> 00:10:51,180 There's no guarantee that this is reasonable. 146 00:10:51,340 --> 00:10:52,740 It's heuristics. 147 00:10:52,900 --> 00:10:54,520 It works in most cases. 148 00:10:54,680 --> 00:11:02,780 The benefit is that now the resulting difference between higher order and low order discrete 149 00:11:02,940 --> 00:11:06,640 evolutions is a difference of states that can be computed. 150 00:11:06,800 --> 00:11:10,080 It's available in the algorithm. 151 00:11:10,240 --> 00:11:14,440 And this is how we obtain this error estimator. 152 00:11:14,600 --> 00:11:19,160 This difference will be used as a replacement for the unknown one step error. 153 00:11:19,320 --> 00:11:21,960 All right, now we have an error estimate. 154 00:11:22,120 --> 00:11:26,900 We have forged a tool that we need to pursue our strategy. 155 00:11:27,060 --> 00:11:32,560 Now let's see how to embed this estimator in a larger algorithm. 156 00:11:32,720 --> 00:11:36,100 Let's discuss local time step adaptation. 157 00:11:36,260 --> 00:11:41,280 Let's assume that we've got an estimate for the one step error at our disposal using this 158 00:11:41,440 --> 00:11:45,220 comparison of two methods with different order. 159 00:11:45,380 --> 00:11:47,560 The approach is the following. 160 00:11:47,720 --> 00:11:50,060 So we compare the estimate. 161 00:11:50,220 --> 00:11:55,660 So this is this ESTK with the tolerances specified by the user. 162 00:11:55,820 --> 00:11:58,200 And now there are two tolerances. 163 00:11:58,360 --> 00:12:04,800 And we have seen this in the termination criteria for iterative methods and in several other 164 00:12:04,960 --> 00:12:05,960 contexts. 165 00:12:06,120 --> 00:12:11,460 We specify a relative tolerance and which is multiplied with the norm of the current 166 00:12:11,620 --> 00:12:13,500 state to get a tolerance. 167 00:12:13,660 --> 00:12:19,620 This expresses the desire to have a control of the relative error, which is the error 168 00:12:19,780 --> 00:12:20,840 that really matters. 169 00:12:21,000 --> 00:12:24,440 But now of course it could be that the current state is actually zero. 170 00:12:24,600 --> 00:12:29,340 Then our algorithm might not terminate. 171 00:12:29,500 --> 00:12:33,800 So we need a second tolerance, the absolute tolerance. 172 00:12:33,960 --> 00:12:38,480 So this is called relative and this is the absolute tolerance. 173 00:12:38,640 --> 00:12:45,240 The absolute tolerance is a safeguard for the case of encountering a zero state. 174 00:12:45,400 --> 00:12:48,800 These tolerances are supplied by the user. 175 00:12:48,960 --> 00:12:54,040 So you give them to the argue to the algorithm as control parameters. 176 00:12:54,200 --> 00:12:56,680 OK, so this is the key decision. 177 00:12:56,840 --> 00:13:00,280 We compare the estimate of the one step error with the tolerance. 178 00:13:00,440 --> 00:13:05,840 Now if it is smaller, say here, yes, this inequality is satisfied. 179 00:13:06,000 --> 00:13:09,280 Then what we do, we accept the current step. 180 00:13:09,440 --> 00:13:15,640 Remember the estimate is computed based on some current time step size. 181 00:13:15,800 --> 00:13:18,880 So we use this time step size and advance. 182 00:13:19,040 --> 00:13:20,880 So this is what we do. 183 00:13:21,040 --> 00:13:22,160 Advance to next steps. 184 00:13:22,320 --> 00:13:25,020 Just carry out one step of the single step method. 185 00:13:25,180 --> 00:13:32,100 And as next state we of course use, well, I should say of course the state provided 186 00:13:32,260 --> 00:13:35,620 by the higher order method and its discrete evolution. 187 00:13:35,780 --> 00:13:41,000 So OK, we have accepted the step because in the current step, the one step error probably 188 00:13:41,160 --> 00:13:43,260 was below the specified tolerance. 189 00:13:43,420 --> 00:13:51,200 Now we show some conservative optimism because now we also slightly increase the current 190 00:13:51,360 --> 00:13:52,880 length of the time step. 191 00:13:53,040 --> 00:13:56,440 So with some factors, say maybe 1.1. 192 00:13:56,600 --> 00:13:58,800 So this is a conservative increase. 193 00:13:58,960 --> 00:14:01,100 Why do we do this? 194 00:14:01,260 --> 00:14:07,160 Because if we accept the current step, it's not unlikely that we would have even accepted 195 00:14:07,320 --> 00:14:08,580 a slightly larger step. 196 00:14:08,740 --> 00:14:14,400 And for the sake of efficiency, we should use steps as large as possible as long as 197 00:14:14,560 --> 00:14:16,860 we still meet the prescribed tolerance. 198 00:14:17,020 --> 00:14:21,360 So this is why this increase of the time step is done for the sake of efficiency. 199 00:14:21,520 --> 00:14:27,320 Now if the good case does not occur, but we find that our estimate is larger than the 200 00:14:27,480 --> 00:14:34,360 tolerance, then we reject the current step and what we do now, we repeat with half the 201 00:14:34,520 --> 00:14:35,520 step size. 202 00:14:35,680 --> 00:14:36,680 So here we are pessimistic. 203 00:14:36,840 --> 00:14:42,280 We believe, well, if the tolerance was not met, then better safe than sorry. 204 00:14:42,440 --> 00:14:45,520 Let's try a significantly reduced step size. 205 00:14:45,680 --> 00:14:48,880 And so here we carry on and here we repeat. 206 00:14:49,040 --> 00:14:55,480 This means we stay at the current point in time, Tk, and try it again with half the step 207 00:14:55,640 --> 00:14:56,760 size. 208 00:14:56,920 --> 00:15:01,000 Let's also look at the actual implementation of these ideas. 209 00:15:01,160 --> 00:15:04,840 Okay, so this is the C++ function. 210 00:15:05,000 --> 00:15:12,380 As arguments, it takes two discrete evolution operators connected with a low order and high 211 00:15:12,540 --> 00:15:13,960 order method. 212 00:15:14,120 --> 00:15:15,860 Usually the order will just differ by one. 213 00:15:16,020 --> 00:15:20,680 And then we pass the two tolerances, relative and absolute. 214 00:15:20,840 --> 00:15:27,600 hmin is meant to stop the algorithm when it tries to reduce the step size indefinitely. 215 00:15:27,760 --> 00:15:32,320 So then of course, we have to jump out of the loop at one point. 216 00:15:32,480 --> 00:15:37,360 All right, so now let's look at the main time stepping loop. 217 00:15:37,520 --> 00:15:44,960 Well, first, we evaluate the discrete evolutions, compute the next states for both methods, 218 00:15:45,120 --> 00:15:49,900 then compute the estimate, the EST and then check. 219 00:15:50,060 --> 00:15:52,380 And this is exactly what I've just discussed. 220 00:15:52,540 --> 00:15:55,860 If the test is passed, then we advance. 221 00:15:56,020 --> 00:15:58,600 So we go to the next time. 222 00:15:58,760 --> 00:16:00,040 So this is advanced. 223 00:16:00,200 --> 00:16:06,200 If we fail the test, then step size is half, we do not advance, and this means repeat. 224 00:16:06,360 --> 00:16:10,880 So well, this advance is done here where we store the next state. 225 00:16:11,040 --> 00:16:19,180 As I told you, hmin is a safeguard against sending step size to zero so that in this 226 00:16:19,340 --> 00:16:21,200 case, the loop will terminate. 227 00:16:21,360 --> 00:16:28,840 If this is activated, this test, then we know that the algorithm has ground to a halt. 228 00:16:29,000 --> 00:16:30,320 When can this happen? 229 00:16:30,480 --> 00:16:35,480 Well, remember the blow up example where the time step has to become smaller and smaller 230 00:16:35,640 --> 00:16:40,160 and smaller to capture this blow up of the states to infinity. 231 00:16:40,320 --> 00:16:46,720 And this finally will abort the main loop until the user now we are close enough to 232 00:16:46,880 --> 00:16:47,880 the singularity. 233 00:16:48,040 --> 00:16:54,280 So I hope this has made clear the main ideas underlying adaptive time stepping. 234 00:16:54,440 --> 00:16:57,060 And now let's look at its performance. 235 00:16:57,220 --> 00:17:03,960 in a few examples, because simple examples are in which we solve scalar ordinary differential 236 00:17:04,120 --> 00:17:07,080 equations and in which we know the analytic solution. 237 00:17:07,240 --> 00:17:12,640 So in the first example, we solve this initial value problem, use explicit Euler as the low 238 00:17:12,800 --> 00:17:18,160 order method and the explicit trapezoidal rule as the higher order method. 239 00:17:18,320 --> 00:17:22,440 So be a method of order one and order two in combination. 240 00:17:22,600 --> 00:17:25,860 And this is, these are plots of the results on the left. 241 00:17:26,020 --> 00:17:32,080 We see the states produced by the adaptive algorithm compared to the exact solution drawn 242 00:17:32,240 --> 00:17:33,240 in green. 243 00:17:33,400 --> 00:17:39,120 And what we see very nicely is that where the solution changes rapidly, so at this internal 244 00:17:39,280 --> 00:17:42,300 layer here, the mesh becomes very fine. 245 00:17:42,460 --> 00:17:49,000 So the pink pluses represent the mesh points and the algorithm does an excellent job about 246 00:17:49,160 --> 00:17:53,680 resolving this steep section of the solution. 247 00:17:53,840 --> 00:17:56,880 And then afterwards, the mesh becomes coarser again. 248 00:17:57,040 --> 00:18:02,360 Only slowly so because of our conservative increase in the step size. 249 00:18:02,520 --> 00:18:05,340 On the right, I've plotted the true error. 250 00:18:05,500 --> 00:18:11,840 So this is the difference of the approximate state and the true state for the time points 251 00:18:12,000 --> 00:18:15,400 and also the estimated error in pink. 252 00:18:15,560 --> 00:18:21,320 What we observe here is that they have little to do with each other at this steep part of 253 00:18:21,480 --> 00:18:22,480 the solution. 254 00:18:22,640 --> 00:18:26,640 So the estimated error and the true error are poorly correlated. 255 00:18:26,800 --> 00:18:31,240 But nevertheless, the adaptive time stepping based on the estimated error does a really 256 00:18:31,400 --> 00:18:35,640 good job about resolving the features of the solution. 257 00:18:35,800 --> 00:18:38,000 This is the magic of this strategy. 258 00:18:38,160 --> 00:18:44,440 There's a very profound explanation for this using control theory, but I'm not going to 259 00:18:44,600 --> 00:18:45,900 discuss this in this course. 260 00:18:46,060 --> 00:18:51,240 So though the estimated error is far off, the true error is still useful for time step 261 00:18:51,400 --> 00:18:51,920 control. 262 00:18:52,080 --> 00:18:52,880 Amazing property. 263 00:18:53,040 --> 00:18:59,920 This experiment discusses to what extent the adaptive strategy really has a benefit. 264 00:19:00,080 --> 00:19:01,340 Same scale OD. 265 00:19:01,500 --> 00:19:05,200 Now we start with initial value zero. 266 00:19:05,360 --> 00:19:07,720 So we are right in the steep part of the solution. 267 00:19:07,880 --> 00:19:12,060 I ran the adaptive algorithm with many different tolerances. 268 00:19:12,220 --> 00:19:14,720 And here you see the states. 269 00:19:14,880 --> 00:19:17,640 The plot doesn't convey so much information. 270 00:19:17,800 --> 00:19:20,660 On the left, on the right, it's much clearer. 271 00:19:20,820 --> 00:19:26,080 What is plotted here is the number of time steps of the adaptive algorithm versus the 272 00:19:26,240 --> 00:19:29,080 error, the maximal error in the time points. 273 00:19:29,240 --> 00:19:31,080 These are the pink stars. 274 00:19:31,240 --> 00:19:37,800 The red pluses give the same information for a uniform temporal mesh. 275 00:19:37,960 --> 00:19:39,440 So equidistant time steps. 276 00:19:39,600 --> 00:19:43,000 And what we see is, let's look at a particular error. 277 00:19:43,160 --> 00:19:49,440 So here this error level means with uniform time step size, we need several hundred time 278 00:19:49,600 --> 00:19:53,480 steps with adaptive time step size clearly below 100. 279 00:19:53,640 --> 00:20:00,540 So there's a big benefit, a big gain in efficiency by using adaptive time stepping in this example. 280 00:20:00,700 --> 00:20:03,520 Now let's look at another example. 281 00:20:03,680 --> 00:20:10,600 Same OD, almost the same initial value problem, but now the initial value is negative. 282 00:20:10,760 --> 00:20:15,920 So the solution has this S shape, as you see on the left. 283 00:20:16,080 --> 00:20:17,920 This is this logistic curve. 284 00:20:18,080 --> 00:20:20,760 Then I varied the relative tolerance. 285 00:20:20,920 --> 00:20:22,120 Now let's see what we get. 286 00:20:22,280 --> 00:20:29,900 Oops, this is a big disappointment because now to achieve a certain level of error, well, 287 00:20:30,060 --> 00:20:33,420 for which the uniform time step has no problem. 288 00:20:33,580 --> 00:20:36,560 The adaptive method doesn't even go so low. 289 00:20:36,720 --> 00:20:40,460 So adaptive time stepping is much worse than uniform time step. 290 00:20:40,620 --> 00:20:42,400 So what's going on here? 291 00:20:42,560 --> 00:20:43,640 The following. 292 00:20:43,800 --> 00:20:48,800 The solution of this initial value problem displays a sensitive dependence on the initial 293 00:20:48,960 --> 00:20:50,600 data. 294 00:20:50,760 --> 00:20:56,880 When you change the initial value only very slightly, the position of this step of the 295 00:20:57,040 --> 00:21:00,320 steep section changes significantly. 296 00:21:00,480 --> 00:21:04,820 And so the adaptive algorithm thinks, okay, in the beginning little happens. 297 00:21:04,980 --> 00:21:08,160 Solution is almost constant, so I can afford to do large time steps. 298 00:21:08,320 --> 00:21:12,060 Large time steps means moderately large local errors. 299 00:21:12,220 --> 00:21:17,960 But now if you make some moderately large local errors in the beginning, say here, then 300 00:21:18,120 --> 00:21:21,120 this will massively shift the position of the step. 301 00:21:21,280 --> 00:21:24,240 So later you're invariably off. 302 00:21:24,400 --> 00:21:28,320 So the adaptive algorithm will use large time steps here. 303 00:21:28,480 --> 00:21:32,560 And then this will mean later that here you have huge errors. 304 00:21:32,720 --> 00:21:40,160 Well, fortunately the picture is not so bleak because this adaptive computation is still 305 00:21:40,320 --> 00:21:44,960 stable in the sense of backward analysis. 306 00:21:45,120 --> 00:21:46,120 What does this mean? 307 00:21:46,280 --> 00:21:51,280 That the apparently bad solution produced by the adaptive algorithm is still a good 308 00:21:51,440 --> 00:21:55,400 solution for slightly changed initial data. 309 00:21:55,560 --> 00:22:01,140 And now remember that input data for numerical algorithms invariably come with errors. 310 00:22:01,300 --> 00:22:05,820 So you can control them up to some relative error only. 311 00:22:05,980 --> 00:22:12,020 And so if you take this into account, even the solution with the significantly shifted 312 00:22:12,180 --> 00:22:13,960 step is okay. 313 00:22:14,120 --> 00:22:20,640 It's the exact solution to data initial values that had been changed only very slightly within 314 00:22:20,800 --> 00:22:25,600 the level of tolerance for the input data. 315 00:22:25,760 --> 00:22:27,640 So the solution is still acceptable. 316 00:22:27,800 --> 00:22:34,560 It's far off in terms of this error measure, but it's okay from the perspective of backward 317 00:22:34,720 --> 00:22:35,720 error analysis. 318 00:22:35,880 --> 00:22:39,760 It's still, as you may call it, a typical solution or acceptable solution. 319 00:22:39,920 --> 00:22:44,960 So this is a sort of rehabilitation of the adaptive algorithm in this setting. 320 00:22:45,120 --> 00:22:51,800 All right, now we have a viable adaptive algorithm, but we can do even better because so far we 321 00:22:51,960 --> 00:22:53,600 waste a lot of information. 322 00:22:53,760 --> 00:22:59,760 We just accept and reject steps based on our estimate, but we don't use the actual value 323 00:22:59,920 --> 00:23:01,260 of the estimate. 324 00:23:01,420 --> 00:23:03,400 So now we want to aim for more. 325 00:23:03,560 --> 00:23:06,520 We don't only want to say accept, reject. 326 00:23:06,680 --> 00:23:13,680 We also want to propose a perfect time step, size, or predict. 327 00:23:13,840 --> 00:23:19,900 And this, of course, should be based on information we have, say based on our estimator or more 328 00:23:20,060 --> 00:23:26,920 precisely on the ratio of the estimated one step error and the tolerance. 329 00:23:27,080 --> 00:23:33,720 So let's recall the formula for the estimator and let's assume that the higher order method 330 00:23:33,880 --> 00:23:37,840 has order p plus one and the lower order method has order p. 331 00:23:38,000 --> 00:23:41,680 So there's a gap of one in the orders. 332 00:23:41,840 --> 00:23:44,560 Now let's look at the one step errors. 333 00:23:44,720 --> 00:23:47,940 So one step error for the lower order method. 334 00:23:48,100 --> 00:23:50,400 Now h is the current time step. 335 00:23:50,560 --> 00:23:54,200 We know that this is an order p method. 336 00:23:54,360 --> 00:23:56,880 So then remember the definition of the order. 337 00:23:57,040 --> 00:24:04,240 So the one step error is proportional or is bounded by h to the p plus one and some remainder. 338 00:24:04,400 --> 00:24:07,480 And this will be h to the p plus two in order. 339 00:24:07,640 --> 00:24:09,560 Of course, everything for h to zero. 340 00:24:09,720 --> 00:24:16,360 Well for the higher order method, we would start here with some constants h to the p 341 00:24:16,520 --> 00:24:21,560 plus two, but we put everything in just a remainder term, p plus two. 342 00:24:21,720 --> 00:24:25,080 Okay, now the estimator is the difference. 343 00:24:25,240 --> 00:24:32,560 Now we can subtract the two equations and we find that the estimator is the h p plus 344 00:24:32,720 --> 00:24:38,040 one and then there's an order h to the p plus two remainder. 345 00:24:38,200 --> 00:24:44,280 But now this is what we sweep under the rug by defining this equality with a dot as an 346 00:24:44,440 --> 00:24:49,600 equality up to terms of order h to the p plus two. 347 00:24:49,760 --> 00:24:51,840 The estimator is known, S k. 348 00:24:52,000 --> 00:24:54,520 This is a computationally available quantity. 349 00:24:54,680 --> 00:25:00,320 And so from this knowledge, we can now get an estimate for the constant. 350 00:25:00,480 --> 00:25:06,840 So we can solve for the constant and this is E s t k h to the p plus one. 351 00:25:07,000 --> 00:25:08,520 Let's label this. 352 00:25:08,680 --> 00:25:14,240 And now we want to use the estimator for predicting the optimal time step. 353 00:25:14,400 --> 00:25:19,320 Optimal in sense of being as large as possible for the sake of efficiency. 354 00:25:19,480 --> 00:25:20,940 Let's call this h star. 355 00:25:21,100 --> 00:25:28,640 And h star should be so big that in fact if you use that h star, then our estimate would 356 00:25:28,800 --> 00:25:30,760 just be the same as the tolerance. 357 00:25:30,920 --> 00:25:34,480 So this is the acceptance limit for the time step. 358 00:25:34,640 --> 00:25:44,000 Okay, but now we can combine this with a formula for the constant and then we find so that 359 00:25:44,160 --> 00:25:54,040 let me compute that c h star p plus one is tall and we replace c with a value that we've 360 00:25:54,200 --> 00:26:01,600 just obtained as h p plus one where h is the current time step. 361 00:26:01,760 --> 00:26:10,600 And so we can solve for h star and obtain that h star is now just h times this quotient 362 00:26:10,760 --> 00:26:17,040 and this can be computed to well the p plus first square root of the quotient. 363 00:26:17,200 --> 00:26:21,640 And now this is the formula for proposing a better time step. 364 00:26:21,800 --> 00:26:28,000 And this h star is now used when we repeat the step instead of halving the time step, 365 00:26:28,160 --> 00:26:29,440 we now use h star. 366 00:26:29,600 --> 00:26:37,080 And if we accept the step, then h star will be our tentative step size for the next step. 367 00:26:37,240 --> 00:26:39,760 So h star serves two purposes. 368 00:26:39,920 --> 00:26:40,960 Great. 369 00:26:41,120 --> 00:26:44,480 I also would like to show you the implementation of this method here. 370 00:26:44,640 --> 00:26:49,940 Well very similar to the previous algorithm but with a new twist that now we have this 371 00:26:50,100 --> 00:26:52,220 formula for the h star. 372 00:26:52,380 --> 00:26:57,340 And well this is now very elegant because when we accept the step then we just move 373 00:26:57,500 --> 00:27:04,080 on we move to the next step and store the current state as an acceptable approximation. 374 00:27:04,240 --> 00:27:07,640 And if we reject the step then we just do nothing. 375 00:27:07,800 --> 00:27:13,280 In both cases we compute the h star and use it either for repetition or for the next step. 376 00:27:13,440 --> 00:27:19,840 And this is what is done in state of the art numerical libraries providing numerical integrators. 377 00:27:20,000 --> 00:27:23,380 Now let's look at how this is really done in the libraries. 378 00:27:23,540 --> 00:27:30,040 Because the numerical libraries they mainly rely on so called embedded Runge Kutta methods. 379 00:27:30,200 --> 00:27:31,520 What is this? 380 00:27:31,680 --> 00:27:34,800 Recall the definition of an explicit Runge Kutta single step method. 381 00:27:34,960 --> 00:27:35,960 We have seen this. 382 00:27:36,120 --> 00:27:42,040 We need sets of coefficients and while executing an implicit Runge Kutta single step method 383 00:27:42,200 --> 00:27:47,720 we have to compute these so called increments in every time step by repeatedly evaluating 384 00:27:47,880 --> 00:27:49,120 the right hand side function. 385 00:27:49,280 --> 00:27:54,360 Well now time step control is based on using two methods with different orders. 386 00:27:54,520 --> 00:28:00,180 We can use two unrelated explicit Runge Kutta methods but then for each we have to compute 387 00:28:00,340 --> 00:28:01,540 all these increments. 388 00:28:01,700 --> 00:28:06,960 And now already back in the 60s computational mathematicians made an observation. 389 00:28:07,120 --> 00:28:14,160 It's possible to use the same set of say i, i, i, j and c, i coefficients to obtain two 390 00:28:14,320 --> 00:28:17,760 methods that just varied in the b, i coefficients. 391 00:28:17,920 --> 00:28:23,520 So you have two Runge Kutta single step methods that rely on the same increments and the only 392 00:28:23,680 --> 00:28:29,680 difference is that then in the final step these increments are combined in a different 393 00:28:29,840 --> 00:28:31,340 way. 394 00:28:31,500 --> 00:28:37,600 They have the same set of a and c coefficients and just differ in the b coefficients. 395 00:28:37,760 --> 00:28:40,640 These methods are called embedded Runge Kutta methods. 396 00:28:40,800 --> 00:28:44,360 So the notation is that you extend the Butcher scheme. 397 00:28:44,520 --> 00:28:51,060 So the Butcher matrix is the same for both methods and it's just these b, i coefficients. 398 00:28:51,220 --> 00:28:59,160 So here b, i and well these are the b, i coefficients and the convention is that the top row of b, 399 00:28:59,320 --> 00:29:04,920 i coefficients describes the discrete evolution of the low order method and this is the method 400 00:29:05,080 --> 00:29:10,640 of order p plus one if the low order method has order p, the one that I marked with a 401 00:29:10,800 --> 00:29:11,820 tilde. 402 00:29:11,980 --> 00:29:15,280 This is one particular example called the Merson embedded method. 403 00:29:15,440 --> 00:29:19,480 Here it's p equals four, p plus one equals five. 404 00:29:19,640 --> 00:29:25,880 And these methods are used in numerical libraries, these embedded Runge Kutta methods. 405 00:29:26,040 --> 00:29:30,320 They are the workhorses for numerical integration nowadays. 406 00:29:30,480 --> 00:29:35,040 So when you use them, typically you have to provide the right hand side function as a 407 00:29:35,200 --> 00:29:41,440 functor then you have to provide an initial value, of course initial time in the of course 408 00:29:41,600 --> 00:29:47,880 in the general non autonomous case and then the two tolerances, absolute and relative 409 00:29:48,040 --> 00:29:54,080 tolerance, a warning here being able to specify tolerances does not mean that you have the 410 00:29:54,240 --> 00:30:00,880 power to impose a strict error threshold and this is the address in the following remark. 411 00:30:01,040 --> 00:30:07,400 The warning is that even though you seem to supply tolerances, this is just a sort of 412 00:30:07,560 --> 00:30:10,400 control for the adaptive strategy. 413 00:30:10,560 --> 00:30:17,140 It does not mean that actually the discretization error stays below these tolerances in some 414 00:30:17,300 --> 00:30:23,240 sense and this is clear because we use these tolerances to control the one step error and 415 00:30:23,400 --> 00:30:30,760 this one step error as we have seen may only be loosely related to the full discretization 416 00:30:30,920 --> 00:30:32,280 error. 417 00:30:32,440 --> 00:30:39,640 Be careful, don't have the illusion that providing tolerances permits you to find control of 418 00:30:39,800 --> 00:30:45,500 the accuracy of the final solution and this is in particular true for initial value problems 419 00:30:45,660 --> 00:30:47,440 where you have sensitive dependence. 420 00:30:47,600 --> 00:30:53,080 We have seen this in an example and with this remark I would like to close this final section 421 00:30:53,240 --> 00:31:00,280 of chapter 11 on now the practical schemes used for numerical integration and provided 422 00:31:00,440 --> 00:31:01,720 by numerical libraries. 423 00:31:01,880 --> 00:31:06,760 I would also like to encourage you to look at the review questions just to make sure 424 00:31:06,920 --> 00:31:08,000 that you got the main points.