 we are discussing numerical methods for ordinary differential equations. So far we have discussed numerical methods for initial value problems. In this lecture we will discuss some numerical methods for 2 point boundary value problems. Let us start our discussion with linear problems consider the linear second order boundary value problem given by minus y double dash plus p of x y dash plus q of x y equal to f of x, where p q and f are given continuous functions and the problem is posed on an interval a b. We are also given 2 conditions now, one is at the point x equal to a at this point the function y takes the value y a and at the point x equal to b the function y takes the value y b and we want the solution y in such a way that it satisfies these 2 conditions at the boundaries of the interval and at all the interior points it has to satisfy this second order linear ordinary differential equation. Since the conditions are specified at 2 different points the problem is said to be a boundary value problem as we did in the initial value problem the first step is to get a discretization or the partition for the given interval a to b. Let us partition the interval into subintervals with the end points as x naught x 1 up to x n plus 1 where each x j is given by x naught plus j into h where h equal to b minus a by n plus 1. We will always assume that a unique solution for the given boundary value problem exists and it is sufficiently smooth. The first step for us to devise a numerical method for this problem is to use the Taylor's theorem to approximate the derivative y dash. You can see that we are using the central difference formula to approximate y dash therefore, it is given by y of x j plus 1 minus y of x j minus 1 divided by 2 h and you have the reminder term now given at some point xi j where xi j lies between x j minus 1 to x j plus 1. Similarly, you can get an approximation for y double dash also we will take again the central difference approximation for y double dash. If you recall in one of our previous lectures on finite difference formulas we have derived this formula using method of undetermined coefficients right. We have also derived this truncation error there and therefore, this formula is familiar to us. Now, what we will do is we will replace y double dash by this central difference formula and y dash by the central difference formula for y dash and thereby the error that we are committing is of order 2 remember we have already divided by h and h square for y dash and y double dash respectively. Therefore, when you replace the central difference formulas in the equation you get a numerical method which is of order 2 now right. So, let us do that when you replace y double dash and y dash by their corresponding central difference formulas we get this equation. This is the central difference method for the second order ODE given here where we have replaced the central difference formula for y double dash in the first term and we have replaced the central difference formula for y dash in the second term and this holds for all j equal to 1 2 up to n. Here we have used the notation y j to represent the approximate value of the solution at y of x j. As we told in the previous slide the order of the method is 2 it is a second order method remember we also have the boundary conditions at y of a remember y of a is nothing but y of x naught right that we have indicated by y naught. Similarly, y of b is equal to y of x n plus 1 and that we use the notation y n plus 1 right. Therefore, we have this finite difference method together with the boundary conditions y naught is equal to y a and y n plus 1 equal to y b remember we are not making any approximation at j equal to 0 and j equal to n plus 1 they are directly taken from the given boundary conditions. If there is no rounding error these are exactly represented right. Let us rewrite this expression in a different way let us collect all the terms of y j minus 1 and keep them together. Similarly, all the terms of y j are gathered and kept as the second term and finally, all the terms of y j plus 1 are grouped and kept as the third term and of course, we have the right hand right. So, we just rearranged the terms and we considered the finite difference method in this form. Now, let us see how this equation looks like when you put j equal to 1 remember this is the finite difference method that we have after rearrangement and this holds for j equal to 1 2 up to n right. For j equal to 1 you can see that the first term is this with j equal to 1 here that is what we are writing here and then y of j minus 1 becomes y naught right. Similarly, the second term is the same here only thing is now we have put j equal to 1 here and y j is now becoming y 1 and similarly, the last term is the coefficient with y 2 remember p q and f are given functions in our problem. Therefore, all this terms are known to us and also f of x 1 is known to us. Now, look at the first term the first term is appearing with y naught which is also known to us from the boundary condition right. Therefore, the first term is fully known to us and so, you can simply push the first term to the right hand side and write the equation for j equal to 1 as this into y 1 remember y 1 is unknown here. This is known and y 2 is unknown equal to this full term is known to us ok thanks to the boundary condition because of that y naught is also known to us. Similarly, let us take j equal to capital N and then you write this equation you can see that the first term is given by this expression into y n minus 1 which is unknown plus this expression with j equal to n into y n again this is unknown and the last term is this expression with j equal to n and now, we have y n plus 1 again you can observe that y n plus 1 is the right side boundary condition and therefore, this is also known to us. So, you can just take it to the right hand side and write the equation for j equal to capital N as the coefficient with y n minus 1 which is unknown and then the corresponding coefficient with y n which is also unknown equal to the full known quantity where y capital N plus 1 is now known to us from the right side boundary condition. Therefore, we have a set of equations like this where the first equation has two terms on the left hand side because we have pushed the first term to the right hand side because of the boundary condition. Similarly, the last equation has only two terms on the left hand side and the right hand side has two terms where this term coming from the boundary condition all the interior points that is j equal to 2 to n minus 1 you have three terms for y j minus 1 y j and y j plus 1. So, for you have two terms y j minus and each j you have the diagonal term the lower diagonal term and the upper diagonal term and all the coefficients are known to us. You can see that the unknowns are y 1 y 2 up to y capital N and all these coefficients are known you can observe that this forms a linear system which can be written as a y equal to b where y is the left hand the approximate solution for our boundary value problem given by y 1 y 2 up to y capital N and what is this coefficient matrix A well the coefficient matrix A is coming from these terms let us use a notation a j to denote this term that is the lower diagonal term that is the lower diagonal term. Let us use the notation b for the diagonal coefficient and c for the upper diagonal coefficient in that way you can see that the first equation has only two terms because you have pushed to the boundary term on the right hand side similarly the last equation has only two terms where the again the boundary term from the right side boundary is pushed to the right hand side all the interior equations have three terms the lower diagonal term diagonal term and the upper diagonal term in that way we got a tridiagonal system if you recall we have discussed an algorithm to solve a tridiagonal system it is Thomas algorithm you can use Thomas algorithm to obtain the solution which is the approximate solution of your given linear boundary value problem you can immediately see that you can develop a python code to obtain the approximate solution of the given linear boundary value problem once you are given y a y b which are the boundary conditions p q and f which are the coefficient functions you can generate the elements of the matrix a which are a j b j and c j you can send this information into the system subroutine for Thomas algorithm and you can get the unknown vector y which is the approximate solution of the linear boundary value problem I hope you can code this method just to have a clarity let us take a very simple boundary value problem and compute the solution manually and see how it works let us take minus y double dash plus y equal to minus x this is our equation we are given the homogeneous boundary condition y of 0 equal to y of 1 equal to 0 let us take h equal to 1 by 4 you can see from the boundary condition that we are interested in solving this problem in the interval 0 to 1 and we have taken the step size as 1 by 4 by this you can see that we will have a 3 cross 3 tridiagonal system here p of x equal to 0 q of x equal to 1 and f of x equal to minus x you can compare this equation with the given general equation of our problem from there you can see this information once you have this information you can go back to the definition of a j b j and c j and you can compute their values for each j right since p and q are constants you can see that a b and c are constants that is they do not change for different j's right only f will change and they are given by this values once you have all this information you can immediately write the tridiagonal system a y equal to b and it is given by this in fact you can use the Gaussian elimination method also to solve this system you can see that the solution is given like this so with this we have completed the numerical method the finite difference method for linear boundary value problem let us move on to non-linear boundary value problem consider the non-linear boundary value problem y double dash equal to f of x comma y comma y dash in general f can be a non-linear function of y and y dash in which case the boundary value problem will be a non-linear boundary value problem and it is posed on the interval a comma b with the conditions that y of a equal to y a and y of b equal to y b y of a and y of b are given to us y a and y b are some real numbers right therefore we have a boundary value problem again now it is a non-linear boundary value problem you can also apply the finite difference method that we have discussed in the previous section but in that case you may land up with a non-linear system of equations for which you may have to use the Newton's method that we discussed in the non-linear chapter but here we will use another interesting method called shooting method before getting into the method let us quickly see the existence and uniqueness theorem of this non-linear boundary value problem we will only state the theorem but the proof of this theorem is not the subject of numerical analysis let f be a continuous function defined on the domain x y z where x belongs to the closed and bounded interval a b and y and z are in R we assume that the partial derivative of f with respect to y and z be continuous on the domain D and further if f y is positive for all x y z in D and f z is bounded in the domain D then a unique solution for the given boundary value problem exists therefore when we devise the numerical method we have to make sure that all these conditions are satisfied by the right hand side function that we are considering otherwise what happens is you may be devising a numerical method for which there may not be any solution in order to avoid such a situation it is more safe for us to work with those functions f that satisfies all this hypothesis anyway this is just a mathematical remark let us go to the method that we are interested in it is called the shooting method main idea of the shooting method is to first choose some eta you may choose it arbitrarily and consider the initial value problem y double dash equal to f of x comma y comma y dash posed on the interval a b remember that this is the same as the equation that you have in your original boundary value problem and now you have initial conditions one initial condition is the same as it is given in your original problem that is y of a equal to y a therefore you are having the left boundary condition from your boundary value problem fixed as the condition at x equal to a in your initial value problem also now you have one more condition because you have second order equation right therefore you need two conditions in order to have a unique solution for that reason you need one more condition in our original problem we have specified that additional condition at the point b therefore it became a boundary value problem now we are fixing that additional condition at the point x equal to a itself but now we are specifying the condition at y dash and we take y dash of a equal to eta remember eta is something that we arbitrarily choose and then fix it here now if you can solve this initial value problem and get the solution then we will denote that solution as y of x comma eta where eta is the parameter and x is the independent variable in your problem an interesting observation here is to see that if you chose eta in such a way that y of b comma eta is equal to y of b then that y will also be the solution of your boundary value problem why because your y already satisfies your equation and also it satisfies the left side condition y of a equal to y a right now if you have chosen your eta such that y of b with that eta is equal to y b then you are completely done with your boundary value problem solution also right but we really do not know what is that eta which gives you the solution y such that y b comma eta equal to y b we really do not know this therefore in general this will define a non-linear equation whose solution is precisely the eta for which you have y b comma eta equal to y b so the idea is to choose a eta solve this initial value problem get the solution y and then plug in that y into this equation and solve this non-linear problem to get the eta remember this is the equation with variable as eta therefore its root will precisely be the value of eta at which y b eta equal to y b right so how are we going to achieve this well you can use any non-linear iterative method that we have introduced in one of our previous chapters to solve this non-linear equation in our case we will use the secant method but before showing you how to set up the secant method let us first worry about how to solve this initial value problem and get this solution because unless you get the solution you cannot go to set up the secant method to solve this non-linear problem right therefore let us first go to solve this initial value problem well you can solve this initial value problem in two ways one is you can replace y double dash by the central difference formula and replace y dash by again its corresponding central difference formula but that may give you a system of non-linear equations therefore another nice way to get the solution of the second order initial value problem is to convert the second order equation to a system of first order equation how can you do that if you would have done a basic course in ODE you would have learnt that any higher order equation can be converted to a system of first order equations here you have second order equation therefore when you converted to a first order equation you will get two first order equations the idea is very simple even if you have not done a course on ODE it is not very difficult for you to understand this what you do is first you define a function z which is equal to y dash ok once you have this then what is z dash z dash is nothing but y double dash right and y double dash is from your original equation is given by f of x comma y comma y dash that is what we are writing here z dash which is actually equal to y double dash equal to f of x comma y comma instead of y dash we will put this z here in that way you have cleverly eliminated y double dash and got this system of two equations that involves only the first order derivative of the unknown variable now instead of posing the second order equation on the interval a b now we will pose this system of first order equations on the interval a b with the same initial condition y of a equal to y a and instead of y dash of a now we will put z of a because that is the notation we are using here right therefore it is z of a equal to eta now you can see that you have two first order equations you can use any numerical method that we have developed for first order initial value problem now you see we have a system of first order equation with initial condition recall we have developed many methods to approximate solution of a first order initial value problem you can use any of that methods like forward Euler method or rangikutta method or any other multi step implicit or explicit methods to approximate solution of this system of first order equations only thing is you have to apply that method twice one for the first equation and another for the second equation to illustrate that we will consider the simplest possible method that is the forward Euler's method and see how to implement the forward Euler method for this system of first order equations there is nothing difficult here you just have to apply the forward Euler method individually to both this equation choose h equal to b minus a by n and consider for j equal to 0 1 2 up to n minus 1 the Euler forward method for the first equation which is given by y j plus 1 equal to y j plus h into z j plus h into z j and similarly you have the forward Euler method for the second equation given by z j plus 1 equal to z j plus h into f of x j y j and z j now you see starting from j equal to 0 you can go up to j equal to n minus 1 at j equal to n minus 1 you have y capital N that is precisely the approximation for y b eta for a given eta so that is approximately equal to y n that you obtained using the Euler method to obtain that you also need to find z j plus 1 for every j because it is a coupled system so once you have y n let us denote it by y n eta because we are choosing a eta and then setting up this Euler method and computing y n therefore your y n will surely depend on the choice of eta that you have. Taken at the beginning of this problem therefore we will use the notation y n comma eta now once for given eta you know how to get y n eta now we want that eta for which y n eta equal to y b right so that is what ultimately we want to have but unfortunately we do not know that value of eta now how to get that well we will try to capture that eta through an iteration we will use secant method to solve this non-linear equation by replacing y b comma eta which is the exact solution of the system now by the approximate solution y n comma eta right now we have gathered all the information to set up the secant method let us now write the algorithm in a systematic way and call it as the shooting method remember we are given a non-linear boundary value problem what is that problem that problem is here once you are given the non-linear boundary value problem you first have to set up a initial value problem right for that you will choose an eta now remember we want to use secant method in secant method you need two initial guesses right therefore we will choose eta naught and eta 1 two initial guesses for the secant method for eta naught we will first solve the initial value problem that is this initial value problem with eta equal to eta naught how we are doing it we are first converting it to a system of first order ODE's and then using the forward Euler method to get y n comma eta naught similarly you give eta equal to eta 1 correspondingly you set up the initial value problem for the first order system use forward Euler method to get y n comma eta 1 right so therefore once you choose eta naught and eta 1 you can get y n eta naught and y n eta 1 right these are the approximations of the solution of your initial value problem with eta naught and eta 1 as parameters once you have these two values you are now ready to set up the secant method remember you have to apply the secant method for the non-linear equation y of b comma eta minus y b equal to 0 right only thing is instead of exact solution coming from your initial value problem now you will plug in this approximate solutions ok thereby you are supposed to set up this as the iterative formula coming from the secant method applied to this equation you should go back to our non-linear equations chapter recall the formula for secant method for a given non-linear equation and come back and see this is the formula here f of eta is given by this ok now what you will do is instead of using the exact value now you will use the approximate value that is y n eta computed using the Euler method similarly you can also apply Runge-Kutta method or any other method trapezoidal method or any predictor character method anything you can use to approximate this solutions here to have simplicity we have used forward Euler method and similarly you also have y n eta naught and you can plug in these values into the secant method formula and get eta 2 now once you have eta 2 again you set up the initial value problem that is you go back to the step 2 and then you set up the initial value problem here now with eta 2 and again you do this process get y n 2 once you have y n eta 2 you can use these two values to get eta 3 from the secant method and like that you can keep on iterating eta's and if this sequence converges you will eventually get that eta for which y b comma eta is equal to y b so that is the idea of the shooting method let us illustrate the shooting method with this simple non-linear ODE we have y double dash equal to minus y dash square it is a non-linear ODE and we are given the boundary conditions as y of 0 equal to 0 and y of 1 equal to 1 let us take h equal to 0.5 because we are going to do the computation manually therefore it is better to take some big h so that there is no much computation involved well to start with we have to choose two initial guesses eta 0 and eta 1 you can choose them arbitrarily here I have chosen them as eta 0 equal to 1 and eta 1 equal to 1.5 the first step is to take eta 0 and use the numerical method to solve the initial value problem well we decided to take Euler method that is forward Euler method therefore we will apply the forward Euler method to the initial value problem with initial condition as y naught equal to 0 that is given here and y dash of 0 that is y dash of 0 that is now z naught because in our initial value problem we have taken z equal to y dash that is why z naught equal to y dash of 0 and that we will take as eta naught the parameter that we have chosen. So, once you have this initial value problem you will go to apply the Euler method component wise that is for the equation y dash equal to z that gives you y 1 and then you apply the Euler method to the second component that is z dash equal to f of x comma y comma z right and that gives you z 1 for the first step. Now once you have the first step you have to go for the second step that is j equal to 0 gives you this remember y 1 is the approximate solution for y of x naught plus h right that is y of x naught is 0 therefore 0 plus x naught h is 0.5 therefore it is 0.5 this y 1 is the approximation of y of 0.5 similarly you need to have y 2 which is an approximation for y of 1 and that is the right hand side boundary therefore with this h you just have to go 2 steps in the Euler method y 2 is given by 0.75 you can check that and that is precisely what you want to have as the value for the right side boundary with eta naught as the parameter right well you may have to also find z 2 but that is not required because for our shooting method we only want this value once you have this value you do not need to find z 2 but you need to find z 1 because that is plugged in the expression for y 2 therefore you need to find z 1 once you have z 1 you can compute y 2 once you have y 2 your purpose is achieved for this particular eta naught therefore you need not compute z 2. Now let us take eta 1 remember we have already computed y 2 of eta naught corresponding to eta naught now we have to find y 2 of eta 1 corresponding to the parameter eta 1 right for that again we will apply the Euler method with the initial condition as y naught equal to 0 and z naught equal to eta 1 you can again find y 1 z 1 plug in y 1 and z 1 into y 2 and get y 2 as this and denote it by y 2 comma eta 1 remember you already got y 2 comma eta naught now you got y 2 comma eta 1 therefore you are now ready to apply the secant method to get eta 2 and that is given by this formula and when you plug in all these values into eta 2 remember we have y 2 comma eta naught minus y b right what is y b y b is 1 that is why we have here y 2 comma eta 1 minus 1 and similarly y 2 comma eta naught minus 1 here right and that reduces to this value and that is your eta 2 once you have eta 2 again you go back to the Euler method with eta 2 and compute y 1 z 1 and then y 2 corresponding to eta 2 once you have y 2 comma eta 2 you can come to compute eta 3 which is equal to eta 2 minus y 2 comma eta 2 minus 1 into eta 2 minus eta naught divided by y let me write here y 2 comma eta 2 minus 1 minus y 2 comma eta 1 minus 1 ok so that will give you eta 3 again once you have eta 3 you will then go to the Euler method to compute y 2 comma eta 3 like that the iteration will keep on going it may be little confusing but you have to carefully understand it once you understand it is very clear how the iteration goes first set up the initial value problem and then solve that initial value problem using Euler method or any other method for setting up the initial value problem you need to choose the eta first time eta naught and eta 1 and once you have the corresponding y values on the boundary then you will come back to the secant method and get the next iteration for eta like that the iterative sequence will go and if the secant method iteration converges then this eta which comes as the limit of the sequence will be the eta for which you have y n comma that eta is equal to y b ok so let us see how this solution looks like for eta 2 you have y 1 is equal to this value z 1 is equal to this value and once you have y 1 and z 1 you will plug in to y 2 to get y 2 comma eta 2 once you have y 2 comma eta 2 you can go to compute eta 3 and the iteration goes on like this let us see how the graph of the solution looks like here the red line represents the exact solution of the given boundary value problem and the blue line gives the approximate solution for each eta computed using the Euler method the first line corresponds to eta naught and for that y 2 comma eta naught is roughly 0.75 and this is the solution computed for eta 1 remember you have y naught y 1 and y 2 which we denote by y 2 comma eta 1 right similarly for the eta 2 you have this graph and this is the solution obtained using the shooting method you can see that as you go on with eta in the secant method your y 2 eta is going more and more closer to the exact boundary condition this is y b for you remember you have taken y b as 1 right so it is trying to approach this condition of course the solution otherwise is not so good because we have taken h is equal to 0.5 if you take h to be something more smaller say 0.05 and all the approximation will be better but we cannot do it manually you can develop a python code and compute the solution for smaller values of h here I will show you the solution computed using the shooting method with h is equal to 0.05 you can see that eta naught we have taken as 1 as usual and eta 1 is taken as 1.5 for eta naught the solution is this and you have y 2 comma eta naught is something approximately 0.69 and this is the graph corresponding to eta 1 and that is giving us the boundary condition as y 2 comma eta 1 which is approximately 0.92 and this blue line corresponds to eta 2 and that gives us the boundary condition as y 2 comma eta 2 is almost 0.99 you can see if you take eta 3 eta 4 and all it will try to approach the exact boundary point which is y b is equal to 1 right. So, this is how the shooting method works it may be little confusing but if you carefully understand it is not very difficult for you to code it given that we have already learnt the coding for Euler method in the last class you can now combine the Euler method and the secant method and I hope you can develop a python code for the shooting method also with this note let us end this lecture. Thank you for your attention.