CHAPTER 08.05: Higher Order and Coupled Ordinary Differential Equations:

Heun’s Method


In this example, we're going to solve a higher-order ordinary differential equation using Heun’s method. Heun's method is a kind of a Runge-Kutta second-order method, and we will see that how we can apply Heun’s method formulas to solve a second-order ordinary differential equation.

 

So, this is the ordinary differential equation which is given to us, is second order in nature. We have initial conditions given to us. We are asked to find the value at y at 0.5. The step size is 0.25. It simply implies that if I start at 0, the first step takes me to 0.25, next step takes me to 0.5. Hence, making it clear that I need to take two steps in order to be able to find the value of y at 0.5. So, the first thing which I need to do is, because all methods such as Euler's method, Runge-Kutta second-order method or Runge-Kutta fourth-order method they're all based on solving first-order ordinary differential equations. So, we do need to reduce the second-order differential equation to first-order differential equation. So, how we're going to do it is we're going to, let's suppose, say dy by dx equals to z. So, that is the substitution we're going to make. So, we're introducing another dependent variable called z, here. And now when I substitute this into the differential equation I'm going to get 2 dz by dx, because I know that hey if I take dz by dx of this, that'll be nothing but the second derivative of y with respect to x. So, plus 3z, plus 5y is equal to 11 e to the power minus x. So, what you have basically seen is that how I have been able to reduce the second-order differential equation to a first-order differential equation now, and what I'm going to do is, I'm going to rewrite this differential equation, so I'm rewriting this differential equation in this form so that it is in the form of the slope being equal to a function of the independent and the dependent variables. So, what we mean by that is as follows: we have our first differential equation, which is dy by dx equal to z. So, the substitution which we made dy by dx equal to z becomes our first differential equation and, let's suppose, we call it a function f1, and it'll be a function of three variables: x is the independent variable, and y is the dependent variable, then z is the dependent variable, which you introduced into the picture. And corresponding to this you will have an initial condition of y sub 0 equal to 7. Because the initial condition is required on the dependent variable which is there in the differential equation.

 

The second differential equation would be dz by dx equal to this quantity here, and this will be some function f2, let's suppose. And again, it is a function of three variables. One independent variable x and two dependent variables y and z. And we have to correspond it to an initial condition, and we call it z sub 0 equal to 13. And where did I get this from? What is z? Z is dy by dx. And we were already given the initial condition that the dy by dx at 0 is 13, and since dy by dx is z, we can write the initial condition like this. So, it's very important that the first thing before you solve any problem by any of the methods which you've learned so far is to rewrite your differential equation in the form, such as given here, so that you can use the methods which you've learned, numerical methods that you've learned, to be solve the problem. So, since by these two differential equations, how will I be able to find the values of y's and z's as I go along from one step to another, we do need to write down the equations. It's a good idea to always do that. So, y sub i plus one will be equal to y sub i, plus one half the slope of y at x i, and this will be the value of the slope at x sub i plus one. And so you're basically taking the average of the two slopes at the point where you are and the point where you want to go, and then you multiply by h in order to be able to find out what the future value of y is. Since z sub i plus 1 is equal to z sub i plus one half the slope of z at x of i, plus one-half k2 of z times h. So, that will be the value of that's how we're going to be able to find the values of y and z. However, what are k1y? k1y is the slope of y, which is given by the function f1, at this particular point. What is k2 of y? Is also f sub 1 but, at the point where we are going in the one step, will be x of i plus h, y sub i plus k1y times h, because you're trying to find the value of y at this point forward by h. And that will be given by this, and then you have to also look for hey what will the value of z at that point and again using the slope of z, and we are able to do that. The same thing we'll have for k1z will be the function f2. And keep in mind that the arguments now are the same for the values of k1z and k1y and same will be the case for k2z, because the only difference is that you are using, you are finding the slope of z now rather than y. But it has to be found at the same point. So, once we have these formulas in place, now it is simply being substituting at the proper spot to be able to find the values of the slopes, and then be able to find the values of the corresponding values of y and z at a point ahead. Now, keep in mind that in order to be able to apply these formulas, one has to find k1y first, and then k1z. You cannot simply go to k2y as the second quantity to be calculated because you need k1z here in order to get k2y. So, this will become the next one, and then this will be the next one to be solved.

 

So, you have to think about that all the formulas might be given in a certain sequence, it doesn't necessarily mean that the same sequence can be used in order to apply your algorithm. So, again, going back to the thing, you will first calculate k1y, then you will calculate k1z, and then you will calculate k2y, and k2z. So, the first, the first step: i is equal to 0. Let's write down what the corresponding values which are known to us. We know that x naught equals 0, y naught is 7, and z naught is equal to 13. And these come from the initial conditions and the points at which those initial conditions are given to us. So, k1y will be equal to the value of the function x1 at these arguments, and since we know these arguments, and since f1 is just z it will be just 13. Let's calculate k1z now: k1z will be now the, the slope of z, which is given by the function f2, and at the same arguments. And it will be 0 7 13 and you will get here will be 11 e to the power minus 0, minus 5 times 7, minus 3 times 13 divided by 2. So, it is simply applying the expression which we have for f2, and substituting the values of xy and z in there. And this value which we'll get is minus 31.5. Once we have done that, let's start calculating our k2y and k2z. So, k2y will be equal to, sorry, we now need to calculate uh k2y, that's right. We're going to calculate k2y. Will be equal to f sub 1, now at a point ahead, and that will be given by this, and this will be equal to f1 and we substitute the values of x naught h, and then, of course, we just calculated k1 of y, which is 13, so we substitute that here, and we just subs, calculated the value of k1z which is minus 31.5. And, so what we are basically doing by doing this is, that we are calculating what the arguments of the function need to be, f1 need to be, in order to be able to calculate the values of the slope of y at a point ahead which is at 0.25. And this one is just z, so just 5.125. Now, let's calculate k2 of z, the k2 of z will be equal to the value of the function x2, at x naught plus h, y naught plus k1y h, z naught plus k1z h and again keep in mind that now we don't need to do the substitutions because they are the same arguments as for the function f1, for k2 of y. So, we can just simply write it hey we need to calculate it at these argument values there. And now we're going to go to the f2 function, which is the slope of z at 0.25, and that will be given 11 e to the power minus 0.25, minus 5 times 10.25, minus 3 times 5.125, divided by 2. And that turns out to be minus 29.03.

 

So, now what do we have is, we have all the values of the slopes which we need, we needed four different slopes, two slopes of y, two slopes of z respect to x, but at two different points, and that's why the total number is four. And now we can calculate the value of y1 and z1 based on these slopes. So, how do we go about doing that? We have y sub-one equal to y naught, plus one-half k1y, plus one half k2y times h, and this one turns out to be the following numbers. So, this is the slope of y at 0, this is the value of the slope of y respect to x at x equal to, x equal to 0.25. So, we take the average of those two and use that to be able to calculate our value of y1, which is the approximate value of y at 0.25. We do need to calculate z1, you might say hey let's go and jump to finding y2, we need z1 because that's going to form the basis of calculating y2, and you will see that in a little bit. So, z1 is z naught, plus one-half k1z, plus one half k2z times h. So, z naught is 13, and k1z is minus 31.5. k2z, which we already calculated, turns out to be that, so you get that. And once you do this you get 5.434. And this is the approximate value of z at 0.25. So, this is the procedure which you have to take and now, in order to be, do this for one more step, uh we'll be, finally be able to calculate what we started to wanting to calculate, which was y sub 2. So, how do we go about doing that? We now, again, say hey i is equal to 1 now, my x1 is equal to 0.25, my y1 is 9.266, and my z1 is equal to 5.434. These numbers I got from the previous step. Now I'll calculate my k1y, which is the value of the function f1 at these arguments, and this will turn out to be just substituting these arguments which I just calculated. And from there I get 5.434 because k1 is nothing but just the value of z. Now, k1z, again keep in mind that I need to calculate k1z, k1z is nothing but the value of the function f2 at the same points that we where we calculated our k1 of y. So, this one here now turns out to be equal to, you just do the substitution like you did before, and you will get minus 27.03. You'll see that you don't need this number because we are just stopping at y2, but if you're going to one more step you would need k1z, but I am calculating these numbers just for purposes of, purposes of showing that hey how do we do from step to step. 

 

So, now we have k2 of y which will be f1 x1 plus h, y1 plus k1y h, z1 plus k1z h, and if we do that we're going to get the following… And the value which you get for the arguments is 0.5 yeah of course because that is the next step, the value will be 0.5, the value of the slope of…Not, the value of y which you're going to use at the point ahead turns out to be 10.62, and the value of the slope, of the slope which is the slope of z, now slope of y, I should say. I shouldn't call it slope of slope, this is the value of the slope of y at um at the point ahead, so that's what we get. And once we have that, we can calculate this value and turns out to be minus 1.324 itself because it is just z. And then let's calculate k2 of z now, k2 of z is the is the value of the function f2, had the same arguments as we had for k2 of y, so we don't need to necessarily go through the same calculations because we'll get the same thing. So, it turns out to be this, and so this value here turns out to be minus 21.23. And again, this value is not needed also for what we are trying to calculate, but I’m showing you to you for purposes of completion.

 

So, again now we have calculate the values of the four slopes, the slope of y with respect to x at two different points, the slope of z with respect to x, at two different points, and we're going to use this information to be able to calculate our y2 and z2. And since we don't need to calculate z2 I'm just going to show you how to calculate y2. So, in this case, your y2 turns out to be y1, plus one half the slope of y at the point where you are, plus one half the slope what you want to be. And this one turns out to be 9.266 plus this quantity here. And this turns out to be 9.780. That's the value of y2. So, that is the approximate value of y at 0.5. So, it is important for you to think about that hey, of course, this is the answer which you are looking for, but that hey what are these intermediate values which we just got for the different slopes? So, for example, the value of 5.434 was the value of the slope of y with respect to x at 0.25, however, this value is the value of some approximate value of the slope of y with respect to x at x equal 0.5. So, it is the average of those two slopes which allows you to find a better estimate than Euler's method, let's suppose, of the value of y at the point you are in of interest to you. So, that's the end of this example and I'm hoping that you will be able to see for yourself how accurate this answer is, based on knowing the exact solution, because an exact solution is available for this kind of problem. And that's the end of this segment.