1 00:00:00,200 --> 00:00:05,900 The current chapter 11 of the numerical methods course is all about numerical integration. 2 00:00:06,060 --> 00:00:11,280 And as you know, this means the search for numerical methods for the approximate solution 3 00:00:11,440 --> 00:00:15,780 of initial value problems for ordinary differential equations. 4 00:00:15,940 --> 00:00:23,560 In the second section of that chapter, I am going to introduce you to a few simple specimens 5 00:00:23,720 --> 00:00:30,980 of the most important class of numerical methods, the so called single step methods. 6 00:00:31,140 --> 00:00:42,320 The methods covered in this section all have a very nice and intuitive, both analytic and 7 00:00:42,480 --> 00:00:44,560 geometric justification. 8 00:00:44,720 --> 00:00:52,200 They are also the oldest methods for the numerical solution of initial value problems. 9 00:00:52,360 --> 00:00:57,840 While dating back to the 18th century when people still had to do these computations 10 00:00:58,000 --> 00:00:59,580 with pen and paper. 11 00:00:59,740 --> 00:01:05,760 So what we want to tackle is of course an initial value problem here written in standard 12 00:01:05,920 --> 00:01:12,920 form with a right hand side function f and specifying an initial value, initial state 13 00:01:13,080 --> 00:01:14,540 I should say. 14 00:01:14,700 --> 00:01:17,760 Why not at initial time t0. 15 00:01:17,920 --> 00:01:23,400 Of course I take for granted that an appropriate state space has already been defined. 16 00:01:23,560 --> 00:01:31,520 I would like to point out a particular practical condition that is often faced in application. 17 00:01:31,680 --> 00:01:36,520 It's the fact that that right hand side function f, which of course we assume to be at least 18 00:01:36,680 --> 00:01:41,640 Lipschitz continuous so that we have existence and uniqueness of solutions of this initial 19 00:01:41,800 --> 00:01:43,800 value problem. 20 00:01:43,960 --> 00:01:50,800 We may face the fact that f is given only in procedural form. 21 00:01:50,960 --> 00:01:55,760 So say as a functor object with the following definition. 22 00:01:55,920 --> 00:02:01,320 And this functor object might be opaque or very complicated so that all we care of access 23 00:02:01,480 --> 00:02:06,240 to is simple point evaluations of f. 24 00:02:06,400 --> 00:02:11,480 So when designing numerical methods we have to keep in mind this restriction. 25 00:02:11,640 --> 00:02:13,680 So now then what is the goal? 26 00:02:13,840 --> 00:02:16,640 We want to find an approximate solution. 27 00:02:16,800 --> 00:02:18,780 Generically the solution is a function of time. 28 00:02:18,940 --> 00:02:22,800 But to begin with we set ourselves a more modest goal. 29 00:02:22,960 --> 00:02:30,920 We are content with the approximation of y of t, my notation for the solution, on isolated 30 00:02:31,080 --> 00:02:32,800 points in time. 31 00:02:32,960 --> 00:02:36,560 Say on a so called time grid or temporal mesh. 32 00:02:36,720 --> 00:02:41,240 Yes let me also give term time grid which is defined as follows. 33 00:02:41,400 --> 00:02:47,400 Just a sequence of isolated points in time starting from initial time up to some well 34 00:02:47,560 --> 00:02:50,800 capital T also called final time. 35 00:02:50,960 --> 00:02:56,760 Well remember as a side remark that for an autonomous ODE we always choose zero as initial 36 00:02:56,920 --> 00:03:01,500 time because we can just exploit the time invariance in this case. 37 00:03:01,660 --> 00:03:07,840 The first method in these among these polygonal approximation methods is definitely the oldest 38 00:03:08,000 --> 00:03:13,580 and it's attributed to Leonhard Euler and therefore it has been called the Euler method. 39 00:03:13,740 --> 00:03:20,340 So again we start from the ODE and then rely on a very natural approximation because what 40 00:03:20,500 --> 00:03:25,940 we do is we replace the derivative with respect to time with a difference quotient. 41 00:03:26,100 --> 00:03:32,420 So y dot is replaced with a difference quotient but with one additional constraint. 42 00:03:32,580 --> 00:03:41,840 We want to use a difference quotient, write this, that is anchored on the mesh that involves 43 00:03:42,000 --> 00:03:48,680 values of say the y function only in mesh points and therefore here I write tk plus 44 00:03:48,840 --> 00:03:51,280 one y of tk. 45 00:03:51,440 --> 00:03:57,680 So this is the difference quotient of width hk where hk now is an abbreviation for the 46 00:03:57,840 --> 00:04:00,400 distance of tk plus one and tk. 47 00:04:00,560 --> 00:04:07,000 I do not assume uniform distances between the mesh points by the way and now this difference 48 00:04:07,160 --> 00:04:14,440 quotient approximates y dot and supposed to approximate y dot in tk. 49 00:04:14,600 --> 00:04:19,240 Of course we don't know the function y but now we have a relationship after replacing 50 00:04:19,400 --> 00:04:27,920 the wiggle with an equality we have a relationship between values of y in the mesh points and 51 00:04:28,080 --> 00:04:35,760 now when we call these values or their approximations yk so it should be approximation of y of tk 52 00:04:35,920 --> 00:04:39,880 then we end up with the following recursion. 53 00:04:40,040 --> 00:04:46,360 The difference of yk and yk yk plus one and yk just hk f tk yk. 54 00:04:46,520 --> 00:04:51,880 So this is what we get out of the difference quotient approximation and now this recursion 55 00:04:52,040 --> 00:05:00,720 this generates a sequence of approximate values yk say well m plus one of course y zero is 56 00:05:00,880 --> 00:05:07,600 just given by the initial state where I would like to emphasize this again yk supposed to 57 00:05:07,760 --> 00:05:14,000 be an approximation of the solution at time point tk and this sequence is generated in 58 00:05:14,160 --> 00:05:21,180 a straightforward way starting from y zero y one y two we just plug this state into the 59 00:05:21,340 --> 00:05:27,080 recursion and get the next one and all it takes to compute this sequence is evaluations 60 00:05:27,240 --> 00:05:33,720 of the right hand side functions and this is why the method is called explicit. 61 00:05:33,880 --> 00:05:42,200 There's nothing to be solved it's straightforward repeated evaluations that give us the sequence 62 00:05:42,360 --> 00:05:48,560 of approximations of the solution in the mesh points and this is the explicit Euler method 63 00:05:48,720 --> 00:05:53,600 as already announced it has a nice and simple geometric interpretation which I would like 64 00:05:53,760 --> 00:06:00,700 to give for the scalar case so capital N equal one one dimensional state space then when 65 00:06:00,860 --> 00:06:07,240 you're in a point state y naught at time t naught you just follow the tangent to the 66 00:06:07,400 --> 00:06:14,680 solution trajectory through that point evaluate the tangent at the next time t1 and this gives 67 00:06:14,840 --> 00:06:23,080 you the next approximation because well the ODE tells us that F provides directly provides 68 00:06:23,240 --> 00:06:29,960 the slope of the tangent to the solution trajectory in well well I should say the graph of the 69 00:06:30,120 --> 00:06:31,460 solution trajectory. 70 00:06:31,620 --> 00:06:40,400 This is the simple geometric interpretation which probably made Euler discover this method 71 00:06:40,560 --> 00:06:45,440 or propose this method so well this interpretation also works in other contexts so let's look 72 00:06:45,600 --> 00:06:52,040 at a non autonomous scalar ODE it's called the Riccati differential equation y dot is 73 00:06:52,200 --> 00:06:59,140 t square plus y square what I've drawn here is the tangent field so for a point in time 74 00:06:59,300 --> 00:07:05,680 and a state which corresponds to a point in in the coordinate system the blue arrow gives 75 00:07:05,840 --> 00:07:12,000 you the direction of and say not only the direction but also the say the sloppily speaking 76 00:07:12,160 --> 00:07:17,600 the length of the tangent and then what is done is we apply the explicit Euler method 77 00:07:17,760 --> 00:07:23,360 and by following the tensions we move on from state to state and the resulting piecewise 78 00:07:23,520 --> 00:07:28,960 linear continuous curve is also known as the Euler polygon by the way drawing this polygon 79 00:07:29,120 --> 00:07:36,840 gives us a possibility to recover a function of time a piecewise linear approximation of 80 00:07:37,000 --> 00:07:44,800 the solution trajectory so then we go beyond merely computing approximate states in times 81 00:07:44,960 --> 00:07:51,000 that are provided by a temporal mesh so we can also do this in higher dimensions well 82 00:07:51,160 --> 00:07:58,480 you may remember from the previous course unit the Lotka Volterra ODE this was a predator 83 00:07:58,640 --> 00:08:06,960 prey model post on the state space of a positive population densities two dimensional so if 84 00:08:07,120 --> 00:08:15,400 this is n equal to and now we can do something similar the pink arrows represent the velocity 85 00:08:15,560 --> 00:08:21,500 field defined by the right hand side of the ODE and we just follow the pink arrow the 86 00:08:21,660 --> 00:08:27,640 direction and the length of the arrow gives us the width of the step and then we construct 87 00:08:27,800 --> 00:08:33,920 a sequence of states also here they can be connected by straight lines and in this way 88 00:08:34,080 --> 00:08:42,040 we obtain trajectories and again what I've drawn in red also in this case is the Euler 89 00:08:42,200 --> 00:08:49,040 polygon you meanwhile you may get bored because this is all repetition of the same idea again 90 00:08:49,200 --> 00:08:54,360 and again but I just would like to demonstrate that this construction principle works well 91 00:08:54,520 --> 00:08:58,960 in various settings for scalar case even the vectorial case where then we get an Euler 92 00:08:59,120 --> 00:09:05,780 polygon in state space all right simple method straightforward to implement with a nice geometric 93 00:09:05,940 --> 00:09:13,560 interpretation but there was something arbitrary about this method and this is revealed when 94 00:09:13,720 --> 00:09:18,780 looking at derivation of the explicit Euler method a bit more closely so well what we 95 00:09:18,940 --> 00:09:24,400 did we replaced the white dot with a difference quotient on the temporal mesh but now you 96 00:09:24,560 --> 00:09:31,560 may wonder why do we use the Tk y of Tk for the right hand side the difference quotient 97 00:09:31,720 --> 00:09:37,960 is essentially symmetric it does not put a special weight on one of its end point so 98 00:09:38,120 --> 00:09:42,760 why not using the other end point here well when you do this then essentially what you 99 00:09:42,920 --> 00:09:50,480 use now is a backward difference quotient because it's anchored at Tk plus one and it 100 00:09:50,640 --> 00:09:54,800 goes backward in time whereas for the explicit Euler method we employed a forward difference 101 00:09:54,960 --> 00:10:02,420 quotient also this gives us a recursion that reads say yk plus one yk again an approximation 102 00:10:02,580 --> 00:10:12,400 of the solution at Tk minus yk is hk f of Tk plus one but now yk plus one well essentially 103 00:10:12,560 --> 00:10:18,080 I'm running out of highlighting options it's just this k plus one here and here that makes 104 00:10:18,240 --> 00:10:24,440 the difference well but now you may object is no longer nicely defines a sequence because 105 00:10:24,600 --> 00:10:32,520 now now we need to solve an equation to compute this yk plus one and this can even be a big 106 00:10:32,680 --> 00:10:38,600 system think of cases where the dimension capital N is big and this is why now this 107 00:10:38,760 --> 00:10:45,180 method is called implicit so the next member of the sequence is defined only implicitly 108 00:10:45,340 --> 00:10:52,800 by the recursion equation but nevertheless it's a recipe to well derive a recursion and 109 00:10:52,960 --> 00:10:58,360 it can be stated like this and this is now the so called implicit Euler method well you 110 00:10:58,520 --> 00:11:05,520 may not be entirely happy when I'm telling you that this defines a recursion and sequence 111 00:11:05,680 --> 00:11:13,520 of approximate state because this system of equation may not have a solution so is is 112 00:11:13,680 --> 00:11:19,800 it possible to solve the implicit equation is the implicit Euler method feasible now 113 00:11:19,960 --> 00:11:26,560 well let me give you some intuition let's think of the case hk is equal to zero well 114 00:11:26,720 --> 00:11:32,480 this of course seems to be pointless because then of course we have a solution namely that 115 00:11:32,640 --> 00:11:40,560 yk plus one is just yk so the state is reproduced but now I can invoke what I would call a continuity 116 00:11:40,720 --> 00:11:48,640 argument in the sense that when we have of course a solution for hk equals zero we can 117 00:11:48,800 --> 00:11:54,680 expect that for small hk small time step there will still be a solution so the intuition 118 00:11:54,840 --> 00:12:01,920 is whenever solution for h equals zero then for small hk a solution of the recursive of 119 00:12:02,080 --> 00:12:08,560 the implicit Euler equation should still exist and there's a rigorous underpinning for this 120 00:12:08,720 --> 00:12:13,080 and the corresponding theorem is called the implicit function theorem and this rigorous 121 00:12:13,240 --> 00:12:20,120 justification employs the following arguments so we start from the defining equation of 122 00:12:20,280 --> 00:12:26,800 the implicit Euler method rewrite it as a problem for finding this zeros of a function 123 00:12:26,960 --> 00:12:33,680 with which is parameterized by h and now we observe that of course well that function 124 00:12:33,840 --> 00:12:41,680 for h equals zero is zero then it no for h equals zero it has a very clear solution because 125 00:12:41,840 --> 00:12:47,760 this is just yk so for h equals zero we have a solution of course and we also find that 126 00:12:47,920 --> 00:12:55,040 the Jacobian of g for h equals zero is just the identity is a regular matrix and moreover 127 00:12:55,200 --> 00:13:02,360 g is of course assuming nice f at least a continuous function and these are the assumptions 128 00:13:02,520 --> 00:13:11,240 of the implicit function theorem which tells us that for sufficiently small h or hk the 129 00:13:11,400 --> 00:13:15,640 defining equation for the implicit Euler method has a unique solution should refer to the 130 00:13:15,800 --> 00:13:21,520 number of course the theorem does not tell you what is sufficiently small this will depend 131 00:13:21,680 --> 00:13:27,180 on well the fine properties of the right hand side function f but at least in the neighborhood 132 00:13:27,340 --> 00:13:32,720 of zero we have a solution for time steps in the neighborhood of zero and this justifies 133 00:13:32,880 --> 00:13:40,280 my saying that the implicit Euler method defines a sequence of approximate states there's always 134 00:13:40,440 --> 00:13:46,640 the assumption of sufficiently small time step now we've seen backward difference quotient 135 00:13:46,800 --> 00:13:51,840 earlier forward different difference quotient giving rise to the Euler methods but this 136 00:13:52,000 --> 00:13:57,160 is not the only possibility the same idea I mean the idea of replacing the derivative 137 00:13:57,320 --> 00:14:02,900 with respect to time with the difference quotient based on m on the mesh the same idea can be 138 00:14:03,060 --> 00:14:09,520 applied in in other versions for instance starting from the ODE we can now well view 139 00:14:09,680 --> 00:14:15,520 this difference quotient as a symmetric difference quotient and as such it approximates the y 140 00:14:15,680 --> 00:14:23,280 dot halfway between times tk and tk plus one and of course also halfway between the states 141 00:14:23,440 --> 00:14:31,920 well and this gives rise to another implicitly defined method when you define the state and 142 00:14:32,080 --> 00:14:40,920 yk plus one by the solution of this equation and well here we have no well yes just copy 143 00:14:41,080 --> 00:14:46,760 this and highlight it so this is the recursion of the so called implicit midpoint method 144 00:14:46,920 --> 00:14:51,640 and the same arguments as for the implicit Euler method confirm that for sufficiently 145 00:14:51,800 --> 00:14:58,040 small time step here again a unique solution of this equation will exist and so the method 146 00:14:58,200 --> 00:15:04,160 is well defined there's also a geometric interpretation here of this implicit midpoint method let's 147 00:15:04,320 --> 00:15:10,080 take a brief look in the end so you start from state y naught at time t naught we look 148 00:15:10,240 --> 00:15:14,880 at the first step but the same principle of course applies to all other steps so it's 149 00:15:15,040 --> 00:15:21,080 sufficient to define the first step and then you don't follow the tangent like in the explicit 150 00:15:21,240 --> 00:15:30,680 Euler method but you follow a line whose slope is given by the tangent direction direction 151 00:15:30,840 --> 00:15:37,600 at a time point between t naught and t1 and the state for which you compute the tangent 152 00:15:37,760 --> 00:15:44,280 direction is of course well the mean of the initial state and the state that you actually 153 00:15:44,440 --> 00:15:50,280 want to compute and this expresses the implicit nature so we define the state that gives us 154 00:15:50,440 --> 00:15:56,160 the tangent direction depending on the state that we actually want to compute well and 155 00:15:56,320 --> 00:16:01,460 this is another way of saying that we have to solve an equation to obtain the next state 156 00:16:01,620 --> 00:16:07,660 so and with this geometric view of the implicit midpoint method I would like to close this 157 00:16:07,820 --> 00:16:07,840 course unit.