 everybody, we have introduced some iterative methods for solving linear systems. We have introduced Jacobi method, Gauss-Seidel method. We have also introduced successive over relaxation method. In this lecture, we will implement the Jacobi method as a python code. Let us quickly recall the Jacobi method. In our theory class, we have derived the Jacobi method in particular for a 3 by 3 system. Let us quickly recall the formula for the Jacobi method. In the case of a 3 by 3 system, we are given a initial guess x naught. In this case, it is a three dimensional vector. Then, the Jacobi iterative sequence can be computed using this formula and this has to be done for each k is equal to 0, 1, 2 and so on. Let us recall the formula. There are two important points to be noted in this formula as far as the coding is concerned. In each of the component of the vector x k plus 1, what we are doing is that we are keeping the diagonal element on the left hand side and moving all the non-diagonal element to the right hand side. Therefore, in each component, we have two terms. In this, we do not include the diagonal term. That is the first point to note and second thing is we are dividing both sides by the diagonal element a i i. In the first equation, it is a 1 1 and second equation, it is a 2 2 and in the third equation, it is a 3 3. Therefore, we have to make sure that all these diagonal elements are non-zero. These are the important points that we have to keep in mind when we are developing the code for Jacobi method. To have a clear idea, let us put this formula in the general case when our system has n linear equations. In this case, our initial guess will be an n dimensional vector, which is generally chosen as per our wish and then for each iteration, you have to compute now n components of the vector x k plus 1 and each component is computed using this formula. Again, you note that the summation on the right hand side should exclude the diagonal term. That is, the summation should run from j equal to 1 to n, where you should not include the term j equal to i here. That is one thing and second thing is when you compute the components, you have to make sure that the diagonal elements are non-zero. Before going to the code, let us see how many loops are involved in this computation. Well, the first loop in our code should run for the iterations and the index is denoted by k. Therefore, our first loop is for the iteration and for each k, we have to compute the vector x k plus 1. Remember, x k plus 1 has n components. Therefore, our second loop should run for the components of the vector x k plus 1. That is your second loop and there is one more loop in each component of the vector x k plus 1. We have to find the value of this sum. Therefore, the third loop should run for j is equal to 1 to n and that is your third loop. Remember, in this third loop, you have to exclude one term. That is the diagonal term, which means you have to always check whether j is equal to i or not. If j is not equal to i, then you have to include the term into your summation. If j is equal to i, then you have to exclude it. We can see that all these loops will be nested loops and therefore, the code will be little confusing. Mathematically, Jacobi method is a very simple and easy to understand method. However, the computer code will be little confusing. Therefore, we will carefully try to understand this. Let us go into the code now. The first part of our code is to include some modules into our program. Here, we are including two modules. One is numpy, which is mainly used to declare the arrays. The second one is the math plot lib, which is used to plot some graphs. Mainly, in our program at the end, we will plot the L2 norm of the error involved in each iteration. For that purpose, I am also importing math plot lib module. The next part of our code is input. Here, we are taking the left hand side coefficient matrix in the variable called a. We are taking the right hand side vector in the variable b. Finally, the initial guess is taken in the variable x. You can see that a is a two dimensional array, whereas b and x are one dimensional array. With this, you can see that we are trying to solve the linear system x1 plus 0 x2 plus x3 is equal to 0 minus x1 plus x2 plus 0 x3 is equal to 0 minus x1 plus x2 plus 0 x3 equal to 0. Finally, x1 plus 2 x2 minus 3 x3 is equal to 0. The initial guess x0 is taken to be 1, 1 and 1. This variable corresponds to the number of iterations that we wish to run. I have taken it as 12. If you want to have more iterations, you can increase the value of capital N. With all these inputs, let us go to the main part of the program. At the first stage, we are trying to see what is the length of the array b. We are taking that into the variable small n. This is to sense what is the dimension of our system. In this particular input, we have taken a 3 cross 3 system. Therefore, we could have directly written 3 here, but I want to make our code general so that you go to enter a matrix A, which is say 4 cross 4 matrix or 5 cross 5 matrix and correspondingly the vectors b and x, then that should be automatically taken into consideration in our main program. That is why I do not want to directly feed the value 3 in this variable n. Rather, I want to find it from the length of the array b, which will precisely be the dimension of our system. Then I am declaring a variable called indicator. We will see why we are using this variable. Then the next part of our code is to check whether all the diagonal elements are non-zero. This is one of the main points, which we have noted from our previous slide that we have to check that all the diagonal elements are non-zero. If at least one diagonal element is zero, then we do not want to go ahead with the Jacobi method. What we are doing is we are simply printing that the Jacobi method fails and then we are putting the value in the variable indicator as zero. Remember, we have initialized the variable indicator as 1. Now, we have put the value zero into the indicator. That says that we should not compute the Jacobi iteration. Let us see how we are sensing that. Let us get into the main part of the Jacobi method. This main part is executed only when the value of the indicator is 1. If the value in the indicator is zero, then the Python control will not enter into the if statement. Since, from our previous for loop, you can see that if any one of the diagonal elements is zero, then the indicator will become zero. Therefore, the control will not enter into this if statement. If all the diagonal elements are non-zero, then this part of the statement is never executed. Therefore, the indicator will hold the value 1, which was initialized previously. Therefore, the Python will enter into this if statement. That is the idea. Remember for the system that we have given, you can see that all the diagonal elements are non-zero. Therefore, we will get into this if loop and we will start computing the iterations. As we have seen, the first loop of our method is going to be for the iteration and it is denoted by k mathematically. The k runs from zero up to capital N minus 1. Remember, capital N is the variable we have used to store how many iterations we want to compute in our method. If you recall, we have taken N is equal to 12 in the input. Therefore, this loop will run for k is equal to 0, 1, 2 up to 11. Then, the first step is to compute the value of this loop is to allocate memory for x. If you recall, we already allocated memory for x 0 and we have stored the vector 1, 1, 1 in that memory. Now, what we are doing is we are allocating another memory for this for the variable x 1, which is the vector from the first iteration of the Jacobi method. When we create the memory for this variable, we are initializing it as 0, 0, 0 and then we will try to compute each of this coordinates using the Jacobi formula. That is the idea. If you recall, what is the Jacobi formula? The Jacobi formula is given like this for the coordinates of the vector x k plus 1. Let us see how to compute that. The next loop is for the components of the vector x k plus 1 and that is our second loop. Let us write the formula step by step and see how to implement it here. If you recall, x k plus 1 i is equal to 1 by a i i and then you have the bracket in that you have b i minus some summation. In that, we are first taking the value of b i into x k plus 1 i here. Remember, when we initialized this vector x k plus 1, remember when k equal to 0, it is simply x 1 and we have initialized it as 0, 0, 0 and then we came into the second loop, which runs through the vector x k plus 1 for each component of the vector x 1. In that, we are first storing the value of b i into x i of 1 and then we are getting into our third loop. What is our third loop? Third loop is for the summation, which is involved on the right hand side of this formula. If you recall, that is sigma j equal to 1 to n and you have to exclude the diagonal element and then it is a i j x j of k, the previous iterations value. So, let us see how to implement this summation. For that, we are running this innermost loop with the index as j. The first line is to check whether j is not equal to i and this summation will happen only if j is not equal to i. Then, what we do is we will take the previous value of a i x k plus 1 i and with that we will subtract a i j into x j. This is something like initially we had x i k plus 1 is equal to b i. That was the first value, which is assigned into x i k plus 1 and then when j equal to 0, what you are doing is you are subtracting a i 0 x 0 of k, the previous iteration. You see you have k here, that is because you are taking the value from the previous iteration and this is like the first component of the ith row and then when it goes to the second time in the loop j becomes 1. Therefore, when it comes for the second time of course, it will check the if condition. If this if condition is satisfied then x i k plus 1 is holding now the value given by this expression. From that value now you are again subtracting a i 1 into x 1 of k and like that this loop will go on till j touches n minus 1. When you come out of this loop you will have the value of this expression for the corresponding ith coordinate. Finally, what you have to do once you finished computing this expression inside the bracket then you have to divide that value by a i i. That is what we are doing now in this step after executing this loop you are coming out of this loop and then taking whatever value you obtained that is precisely the value of this expression inside the bracket and then you are dividing it by a i i. You see here you are dividing it by a i i that completes the calculation of one component of the vector x k plus 1 that is x k plus 1 i and this has to go for each i equal to 0 1 2 up to n minus 1. Once that is done then you get the vector x k plus 1. Once x k plus 1 is done then you will go back to the outermost loop you will increment k by 1 more and then you will initialize the memory for the new x which is the vector of the next iteration and then you come to compute the value of the components of that vector like that the loop will go on. How many times well the number of times that we have specified in the variable capital N. Once the iterations are calculated for n number of times then python will come out of this entire loop and print this statement which says that approximate solution after that n number of times remember we have given the value of capital N as 12. Therefore, it has to print approximate solution after 12 iterations equal to it will print the vector x capital N that is whatever the iteration it has calculated at the nth stage will be printed as the output. Next we are also interested to have a feeling on how the error looks like for that what we are doing here is that we are calculating the l 2 norm of the error involved in the k th iteration that is we are computing e k and then taking the l 2 norm of it. First let us see what is the error involved in the k th iteration by definition the error involved in the k th iteration is equal to the exact value that is the exact solution of our system minus the vector computed at the k th iteration right. Recall that we have taken the right hand side vector of our system as the 0 vector therefore, the exact solution is 0 because our matrix A is invertible and hence the error involved in the k th iteration is simply minus x k only. Now, let us see what is the formula for the l 2 norm of any vector in particular what is the l 2 norm of e k 2 that is nothing back square root of e first component square plus e second component square and so on till the n th component of the vector e k square right. So, that is the formula for the l 2 norm and you can observe that we have computed l 2 norm for each k is equal to 0 to n capital n as a first step we are initializing the variable e r r l 2 this is going to have the l 2 norm of the error of all the iterations. We are initializing this variable with 0 remember it is created as an array of dimension capital n plus 1 and then for each k th iteration from 0 to n therefore, you can see that this for loop will run for n plus 1 times for each time it will compute the l 2 norm of the error involved in the k th iteration of the Jacobi method and then once it computes the inner sum then it takes the square root and that will finally give you the l 2 norm of the error at the k th iteration. I am also printing the l 2 norm of the last iteration here and finally, we are plotting the graph of e r r l 2 as a function of k that is done in this part of our code I will leave it to you to see how this is done there is nothing to explain here. Let us see the output of this code the code is executed and let us see the output recall that as the first stage we are printing the approximate solution after 12 iterations. So, that is what we expect to be printed as the first output of our code you can see that it is printed here approximate solution after 12 iteration is given by this vector you can see that it is still far away from the 0 vector recall the exact solution of the system that we have taken as a input has the exact solution as the 0 vector. Here you can see that after 12 iterations it has only given this as the approximate solution which is rather a bad approximation and also you can see what is the l 2 error involved in this obviously it is around 0.8. Let us see the plot of the l 2 norm of the error and that is given here the iteration started somewhere here and the error kept on decreasing well it came up to 0.8 around and then it stopped because we asked the code to compute only 12 iterations. Maybe if you go on further in the iteration you may get better and better approximation to our solution let us also see that let us do say 50 iterations and see how the approximation is improved. Let us run the code with capital N is equal to 50 and see what is the improvement in our result. Now you can see that the l 2 error is almost 0.09 also let us see the plot of the l 2 norm of the error you can see that we were somewhere here and then further the error decreased and going closer and closer to 0. Let us go further in the iteration and see what is happening. Let us say go up to 100 iterations and see what is the result now. You can see that the approximation is quite good now the l 2 norm of the error is around 0.005 you can also see the solution the solution is now pretty close to the 0 vector. So, you can see here the computed solution is pretty close to the exact solution and we can also see that the error is considerably reduced and going very close to 0. Now if you go on with the iteration it will go more and more close to 0. Well let us also see what happens if one of the diagonal elements is 0. Say for instance let us take a 2 2 as 0 and see what is happening let us execute the program. Now you can see that the code hit the error message it just flashed the error message that Jacobi iteration failed and it did not get into the main loop that is because a 2 2 is 0 therefore it satisfied this if condition. So, Python entered into this if statement and printed this command and put indicator as 0. Since indicator is 0 this statement is false and therefore Python could not enter into the main part of the program that is Python could not further compute the iterations of the Jacobi method and therefore it could not also print the solution and so on. It just entered the program by printing the error message only and that is why we only see the error message as the output here. Well there are 2 main disadvantages in this program. The first disadvantage is the way we have managed the memory. You can see that we are allocating a fresh memory for x at every iteration that is we have first initialized x 0 and then when you come to compute x 1 you have allocated a fresh memory for this x 1 and then stored the computed values of the components of x 1 in that memory and when you go to the second iteration again you created a fresh memory for x 2 and store the values of the components of x 2 in this memory and so on. Every time you go for the fresh iteration you are creating a memory for that. This can be quite costly if your system is very large and you want to do many iterations then you need lot of memory while running this code but if you carefully observe once you compute x 1 you no longer need the values stored in x 0 right because x 2 is computed by only keeping x 1 values we never use x 0 values right. In general to compute x k plus 1 we only need x k we do not need x k minus 1 x k minus 2 and so on right. Therefore remembering all these computed values is unnecessary as long as the Jacobi method computation is concerned unless you want to remember all these vectors from the method you really do not need to remember all these vectors. There is one more main disadvantage in our code if you see we are asking the code to run the iterations for some pre assigned number of times. Now we really do not know after computing the iterations this many times whether we have achieved the accuracy that we want right. So, for instance at the first stage we have given only 12 iterations and we have not achieved a good accuracy in this particular example we know the exact solution. Therefore, we can see that we have not achieved a good accuracy and so we increase the number of iterations from 12 to 50 then we saw that it is good but still not enough. So, we increased from 50 to 100 like that you can go on but for that you need to know the exact solution right. If you do not know the exact solution then blindly you have to specify some big number it may happen that after computing so many iterations you might have not still achieved the required accuracy or you might have achieved the accuracy much before than what you specified here. So, in that way going with giving number of iterations and computing the approximation is rather a blind approach we need a good stopping criteria that is we want some condition that will tell us that we have achieved the required accuracy therefore, we can stop the iteration at this stage. Let us see how to design a good stopping criteria for that let us see what are all the ways that we can measure the error involved in our computed solution. Let x star denotes the computed solution from any method remember this is not only restricted to Jacobi method what we are discussing here as a stopping criteria can be used for any method. So, we are assuming x star to be the computed solution from some method then we know that we know that the error involved in x star when compared to x is defined as x minus x star and it is denoted by E. Now, remember you may also go for finding the L 2 norm or L infinity norm or L 1 norm of this error and then check whether this is less than say for instance 10 to the power of minus 5 or any tolerance parameter that you specify say for instance you can give some epsilon value to your code and every time once you finish the iteration you can check for the condition that the error involved in that iteration whether it is less than that epsilon or not right, but unfortunately this cannot be implemented practically because to check this condition you need to know the exact solution right. There is another way of sensing this error that is using the residual vector residual vector is defined as or is equal to b minus a into x star. You can see that if x star is very close to x then a x star will be very close to a x which is actually equal to b right. Therefore, when x star is very close to the exact solution then r will be very close to 0 just like how x if x star is very close to x then e will be very close to 0. Similarly, even r will go closer and closer to 0 if x star goes closer and closer to the exact solution. You can also see an interesting fact that r can be written as a x because simply a x is equal to b minus a x star and since a is a linear transformation you can also write this as a of x minus x star and therefore, this can be written as a e equal to r. An interesting observation here is that the error e can be computed as a solution of this system and therefore, you can compute e without knowing the exact solution. Remember if you want to use this formula to compute e you need to know the exact solution whereas, by computing e from this system you just do not need to know the exact solution you can compute e from all the available information. Remember r is fully known to you because x star is known to you b is known to you and a is also known to you. Therefore, r can be computed explicitly and then it can be plugged in as the right hand side vector and you can solve this system to get e without knowing the exact solution of your original system. Coming back to the stopping criteria as I told ultimately we would like to go for checking this norm to be less than epsilon where epsilon is some preassigned positive quantity which is the tolerance of your error, but this is not practically applicable because you just do not know this x or if you want to go for this condition then you may have to use this system solve this system to get e that may be equivalent to solving your original system itself. An alternate idea is to compute the residual vector which can be computed directly from your approximate solution and then check for this condition. Say for instance in our code we could have computed this with a l 2 norm and we could have checked whether that quantity is less than some number say 10 to the power of minus 5 or whatever the tolerance number that you give. You can check this condition after computing each iteration that is for each k once you finish the iteration for that particular k you can compute the residual error and then find the l 2 norm of that residual error and check this condition. If that condition is satisfied you come out of the loop otherwise you go for the next iteration in that way you can get the solution as accurate as you want and you can stop the iteration as soon as you achieve your accuracy. So, this is one popular way to stop your iteration computationally. However, it has a danger that if you are working with a system where a is almost a singular matrix then you may have a very small residual error, but still you may be far away from the solution. Therefore, one has to be very careful in using this stopping criteria. However, if your a is a nice matrix then this is a good stopping criteria with this note let us finish this lecture. Thank you for your attention.