====== 3. Gaussian elimination algorithm ====== Gaussian elimination algorithm allows us to solve systems of linear equations. Let us have a system of linear equations: $$ \begin{pmatrix} a_{11}&a_{12}&\dots &a_{1n}\\ a_{21}&a_{22}&\dots &a_{2n}\\ \vdots & \vdots &\ddots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mn} \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_m \end{pmatrix}= \begin{pmatrix} b_1\\ b_2\\ \vdots\\ b_m \end{pmatrix} $$ and its augmented matrix: $$ \left(\begin{array}{rrrr|r} a_{11}&a_{12}&\dots &a_{1n}&b_1\\ a_{21}&a_{22}&\dots &a_{2n}&b_2\\ \vdots & \vdots &\ddots&\vdots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mn}&b_m \end{array}\right). $$ The following operations with a matrix are called **elementary row operations**: - Switching two rows. - Multiplying a row by a nonzero constant. - Adding to some row a nonzero constant times another row. The matrix $\begin{pmatrix}-2g & -2h & -2i\\a+d & b+e & c+f\\a & b & c\end{pmatrix}$ has been obtained from the matrix $\begin{pmatrix}a & b & c\\d & e & f\\g & h & i\end{pmatrix}$ by elementary row operations.\\ Name the elementary row operations used. ++Answer | $\begin{pmatrix}a & b & c\\d & e & f\\g & h & i\end{pmatrix}\overset{R_2+R_1}\longrightarrow\begin{pmatrix}a & b & c\\a+d & b+e & c+f\\g & h & i\end{pmatrix}\overset{-2\cdot R_3}\longrightarrow\begin{pmatrix}a & b & c\\a+d & b+e & c+f\\-2g & -2h & -2i\end{pmatrix}\overset{R_1\leftrightarrow R_3}\longrightarrow\begin{pmatrix}-2g & -2h & -2i\\a+d & b+e & c+f\\a & b & c\end{pmatrix}$ ++ Write an IPython function for **(a)** adding to a constant $l$ times row $i$ a constant $k$ times row $j$ in a given matrix, return the new matrix; **(b)** switching rows $i$ and $j$ in a given matrix, return the new matrix. \\ from numpy import array def add_row_multiples(M, i, j, l, k): # M is matrix, i is row multiplied by l and add k times row j if (i > len(M) or j > len(M)): # check do the rows i and j exist in the matrix M return 'The row(s) do not exist' else: R = array(M) # the result matrix obtains the value of the initial matrix R[i,:] = l*M[i,:] + k*M[j,:] # row i in R obtains the value l times row i added k times row j return R from numpy import array def swap_rows(M, i, j): # M is matrix, rows i and j are swapped if (i > len(M) or j > len(M)): # check whether the rows i and j exist in the matrix M return 'The rows do not exist' else: R = array(M) # the result matrix obtains the value of the initial matrix R[[i,j],:]=R[[j,i],:] # matrix rows i and j are swapped return R \\ In terms of the augmented matrix, the first elementary row operation means that we can reorder equations within the system.\\ The second elementary row operation means that we can multiply an equation with some constant, i.e. $ \begin{align*} a_{i1}x_1+\cdots+a_{in}x_n=\beta\qquad\Longleftrightarrow\qquad \gamma a_{i1}x_1+\cdots+\gamma a_{in}x_n=\gamma\beta_i\enspace. \end{align*} $ The third elementary row operation means adding two equations: $ \begin{align*} \begin{aligned} a_{k1}x_1+\cdots+a_{kn}x_n=\beta_k\\ a_{l1}x_1+\cdots+a_{ln}x_n=\beta_l\\ \end{aligned} \qquad\Longleftrightarrow\qquad (a_{k1}+\gamma a_{l1})x_1+\cdots+(a_{kn}+\gamma a_{ln})x_n=\beta_k+\beta_l.\\ \end{align*} $ **Theorem.** The elementary row operations deliver an augmented matrix for a system of equations which has the same solution set as the original system. For example, let ;#; $ \begin{align*} \left(\begin{array}{rrr|r} 1 & 3 & 6 & 25\\ 2 & 7 & 14 & 58\\ 0 & 2 & 5 & 19\\ \end{array}\right) \quad\overset{R_2-2\cdot R_1}\longrightarrow\quad \left(\begin{array}{rrr|r} 1 & 3 & 6 & 25\\ 0 & 1 & 2 & 8\\ 0 & 2 & 5 & 19\\ \end{array}\right)\enspace, \end{align*} $ ;#; then ;#; $ \begin{align*} \left\{ \begin{aligned} 1x + 3y + 6z &= 25\\ 2x + 7y + 14z &= 58\\ 2y + 5z &= 19 \end{aligned} \right. \end{align*}\quad $ and $ \quad\begin{align*} \left\{ \begin{aligned} 1x + 3y + 6z &= 25\\ 1y + 2z &= 8\\ 2y + 5z &= 19 \end{aligned} \right. \end{align*} $ ;#; have one and the same solution ($x=1,\ y=6,\ z=3$). ===== 3.1 Simple Gaussian elimination algorithm ===== Let $m=n$ (i.e. the number of equations equals to the number of unknown) and let the system of linear equations have precisely one solution, that is, there is only one set of values $x_1=c_1$, $x_2=c_2$, ..., $x_m=c_m$ for which all the equation in the system are true. Then, the idea of the Gaussian elimination algorithm is $$ \left(\begin{array}{rrrr|r} a_{11}&a_{12}&\dots &a_{1m}&b_1\\ a_{21}&a_{22}&\dots &a_{2m}&b_2\\ \vdots & \vdots &\ddots&\vdots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mm}&b_m \end{array}\right)\qquad\overset{\text{elementary row operations}}{\longrightarrow \cdots\longrightarrow}\qquad \left(\begin{array}{rrrr|r} 1&0&\dots &0&c_1\\ 0&1&\dots &0&c_2\\ \vdots & \vdots &\ddots&\vdots&\vdots\\ 0&0&\dots &1&c_m \end{array}\right), $$ i.e. the augmented matrix of a system of linear equations is converted into such ** diagonal form** by using the elementary row operations. Note that the diagonal form corresponds to the following system of linear equations ;#; $ \begin{align*} \left\{ \begin{aligned} x_1&=c_1\\ x_2&=c_1\\ \cdots&\\ x_n&=c_m \end{aligned}\right. \end{align*} $ ;#; and thus, the last column in the diagonal form obtained from the initial augmented matrix by elementary row operations, i.e. $\begin{pmatrix}c_1&\cdots &c_n\end{pmatrix}^T$, is the solution of the original system (see [[03_gaussian_elimination#t_gauss|Theorem]]). Let us now use Gaussian elimination algorithm for solving the system ;#; $ \begin{align*} \left\{ \begin{aligned} x_2-10x_3&=0\\ 3x_1-x_2-5x_3&=9\\ -2x_1+x_2+3x_3&=-6 \end{aligned}\right. \end{align*}. $ ;#; The augmented matrix of this system is $$ \left(\begin{array}{rrr|r} 0&1&-10&0\\ 3&-1&-5&9\\ -2&1&3&-6 \end{array}\right). $$ **STEP 1. (a)** Switch the first row of the augmented matrix with rows below till the first element in the first column is nonzero: $$ \left(\begin{array}{rrr|r} 0&1&-10&0\\ 3&-1&-5&9\\ -2&1&3&-6 \end{array}\right)\quad \overset{\text{switch}\ R_1\ \text{and}\ R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 3&-1&-5&9\\ 0&1&-10&0\\ -2&1&3&-6 \end{array}\right). $$ **(b)** By using the first element of the first column we zero out the rest of the entries in the first column: $$ \left(\begin{array}{rrr|r} \color{blue}3&-1&-5&9\\ \color{red}0&1&-10&0\\ \color{red}{-2}&1&3&-6 \end{array}\right) \quad\overset{3R_3+2R_1}{\longrightarrow}\quad \left(\begin{array}{rrr|r} \color{blue}3&-1&-5&9\\ \color{red}0&1&-10&0\\ \color{red}0&1&-1&0 \end{array}\right) $$ Verify the step in IPython by using the functions you created for elementary row operations (e.g., [[03_gaussian_elimination#el_row_op|swap_rows]] and [[03_gaussian_elimination#el_row_op|add_row_multiples]]).\\ from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j in matrix M if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R M = array([[0,1,-10,0],[3,-1,-5,9],[-2,1,3,-6]]) # initial matrix i = 1 R = array(M) # the result matrix while (R[0,0] == 0 and i < len(R)): # swap rows till the element a_11 is nonzero R = swap_rows(R,0,i) i = i + 1 print R # matrix with swapped rows for j in range(1, len(R)): # loop for making elements a_i1 zero from the row two if R[j,0] <> 0: l = R[0,0] k = R[j,0] if cmp(l,0) <> cmp(k,0): # choosing constants for row multiplication R = add_row_multiples(R,j,0,abs(l),abs(k)) else: R = add_row_multiples(R, j, 0, l, -k) print R # matrix where 1st column is zeroed out below 1st row \\ Solve the task in IPython for $\mathbb{Z}_{11}$. **STEP 2. (a)** Now switch the second row with rows below till the second element in the second column is nonzero, i.e. apply STEP 1 (a) for the submatrix obtained by ignoring the first row and the first column: $$ \left(\begin{array}{rrr|r} \color{lightgray}3&\color{lightgray}{-1}&\color{lightgray}{-5}&\color{lightgray}9\\ \color{lightgray}0&1&-10&0\\ \color{lightgray}0&1&-1&0 \end{array}\right). $$ In our case the second element of the second row is already nonzero, thus no switching is needed. **(b)** By using the second element of the second column we zero out the rest of the entries in the second column: $$ \left(\begin{array}{rrr|r} 3&\color{red}{-1}&-5&9\\ 0&\color{blue}1&-10&0\\ 0&\color{red}1&-1&0 \end{array}\right) \quad\overset{R_1+R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 3&\color{red}{0}&-15&9\\ 0&\color{blue}1&-10&0\\ 0&\color{red}1&-1&0 \end{array}\right) \quad\overset{R_3-R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 3&\color{red}0&-15&9\\ 0&\color{blue}1&-10&0\\ 0&\color{red}0&9&0 \end{array}\right). $$ Verify the step in IPython (like in step 1).\\ from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j in matrix M if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R M = array([[3,-1,-5,9],[0,1,-10,0],[0,1,-1,0]]) # initial matrix i = 2 R = array(M) # the result matrix while (R[1,1] == 0 and i < len(R)): # swap rows till the element a_22 is nonzero R = swap_rows(R,1,i) i = i + 1 print R # matrix with swapped rows for j in range(len(R)): # loop for making elements a_i2, except i=2, zero if (R[j,1] <> 0 and j <> 1): l = R[1,1] k = R[j,1] if cmp(l,0) <> cmp(k,0): # choosing constants for row multiplication R = add_row_multiples(R, j, 1, abs(l), abs(k)) else: R = add_row_multiples(R, j, 1, l, -k) print R # matrix printed if row i is changed \\ Solve the task in IPython for $\mathbb{Z}_{11}$. **STEP 3.** The third element in the third row is nonzero (if it were zero, the system would have more then one solution, see [[03_gaussian_elimination#gen_gaussian|General Gaussian elimination algorithm]]) and can be used to zero out the rest of the elements in the third column: $$ \left(\begin{array}{rrr|r} 3&0&\color{red}{-15}&9\\ 0&1&\color{red}{-10}&0\\ 0&0&\color{blue}{9}&0 \end{array}\right) \quad\overset{9R_1+15R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 27&0&\color{red}{0}&81\\ 0&1&\color{red}{-10}&0\\ 0&0&\color{blue}{9}&0 \end{array}\right) \quad\overset{9R_2+10R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 27&0&\color{red}0&81\\ 0&9&\color{red}0&0\\ 0&0&\color{blue}{9}&0 \end{array}\right). $$ Verify the step in IPython (like in step 2).\\ from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j in matrix M if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R M = array([[3,0,-15,9],[0,1,-10,0],[0,0,9,0]]) # initial matrix i = 3 R = array(M) # the result matrix while (R[2,2] == 0 and i < len(R)): # swap rows till the element a_33 is nonzero R = swap_rows(R,1,i) i = i + 1 print R # matrix with swapped rows for j in range(len(R)): # loop for making elements a_i3, except i=3, zero if (R[j,2] <> 0 and j <> 2): l = R[2,2] k = R[j,2] if cmp(l,0) <> cmp(k,0): # choosing constants for row multiplication R = add_row_multiples(R, j, 2, abs(l), abs(k)) else: R = add_row_multiples(R, j, 2, l, -k) print R # matrix printed for every i if a_i3 is changed \\ Solve the task in IPython for $\mathbb{Z}_{11}$. **THE RESULT.** It remains to multiply each row by the inverse element of the value of the leading entry: $$ \left(\begin{array}{rrr|r} \color{red}{27}&0&0&81\\ 0&\color{red}{9}&0&0\\ 0&0&\color{red}{9}&0 \end{array}\right) \quad\overset{\frac{1}{27}R_1,\ \frac{1}{9}R_2,\ \frac{1}{9}R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} \color{blue}{1}&0&0&3\\ 0&\color{blue}{1}&0&0\\ 0&0&\color{blue}{1}&0 \end{array}\right), $$ and thus, $x_1=3,\ x_2=0,\ x_3=0$. In general, the $i^{th}$ step of the simple Gaussian elimination algorithm can be described as follows: **(a)** if $a_{ii}=0$ and $i>m$ switch the $i^{th}$ row with a row $j$, where $j>i$ till $a_{ii}$ is nonzero; **(b)** use $a_{ii}$ to zero out all of the other elements in the $i^{th}$ column. Write a function in IPython for performing the $i^{th}$ step of the algorithm.\\ from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j in matrix M if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R ############################################################################################################# def gaussian_step(M, i): if (i < 1 or i >len(M)): return 'No such step in this matrix' R = array(M) # the result matrix s = i while (R[i-1,i-1] == 0 and s < len(R)): # swap rows till the element a_ii is nonzero R = swap_rows(R, i-1, s) s = s + 1 print R # print matrix with swapped rows if R[i-1,i-1] == 0: print 'Error: the system does not have one solution' return numpy.zeros(shape=(len(R),len(R.T))) # if we continue, we have row multiplication with zero for j in range(len(R)): # loop for making elements in row i, except a_ii, zero if (R[j,i-1] <> 0 and j <> (i - 1)): l = R[i-1,i-1] k = R[j,i-1] if cmp(l,0) <> cmp(k,0): # choosing constants for row multiplication R = add_row_multiples(R, j, i-1, abs(l), abs(k)) else : R = add_row_multiples(R, j, i-1, l, -k) print R # matrix printed every time it changes return R \\ Write a function in IPython for solving a system of linear equation with the simple Gaussian elimination algorithm (mind the assumptions for the algorithm; use the function created in the previous exercise (e.g., gaussian_step)).\\ from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j in matrix M if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R def gaussian_step(M, i): # function for step i in gaussian elimination algorithm if (i < 1 or i >len(M)): return 'No such step in this matrix' R = array(M) s = i while (R[i-1,i-1] == 0 and s < len(R)): R = swap_rows(R, i-1, s) s = s + 1 print R if R[i-1,i-1] == 0: print 'Error: the system does not have one solution' return numpy.zeros(shape=(len(R),len(R.T))) for j in range(len(R)): if (R[j,i-1] <> 0 and j <> (i - 1)): l = R[i-1,i-1] k = R[j,i-1] if cmp(l,0) <> cmp(k,0): R = add_row_multiples(R, j, i-1, abs(l), abs(k)) else : R = add_row_multiples(R, j, i-1, l, -k) print R return R ############################################################################################################ def simple_gaussian(M): R = array(M) # the result matrix for i in range(len(R)): R = gaussian_step(R,i+1) if (R == numpy.zeros(shape=(len(R),len(R.T)))).all() : print 'Error: not solvable with this method' return T = R.astype(numpy.float) # needed when initial matrix is int, but solution is float for j in range(len(R)): T[j,:] = (float(1)/R[j,j])*R[j,:] # multiplying matrix rows so that on diagonal we get 1 print T return T[:,len(R)] # returns last column of the matrix, i.e. the solution \\ Find a solution to the equation $A\mathbf{x}=\mathbf{b}$ when $$\begin{align*}A=\begin{pmatrix}1 & 1 & 1\\1 & 2 & 4\\1 & 3 & 9\end{pmatrix}\end{align*}$$ and $\mathbf{b}=(1\ 1\ 1)^T$ over $\mathbb{Z}_2$, $\mathbb{Z}_3$ and $\mathbb{Z}_5$. ===== 3.2 General Gaussian elimination algorithm ===== Let us now use the simple Gaussian elimination algorithm for the system ;#; $ \begin{align*} \left\{ \begin{aligned} 2x_1-3x_2+5x_3&=1\\ 4x_1-6x_2+2x_3&=2\\ 2x_1-3x_2-11x_3&=1\\ \end{aligned}\right.\enspace. \end{align*} $ ;#; The first step is $$ \left(\begin{array}{rrr|r} 2&-3&5&1\\ 4&-6&2&2\\ 2&-3&-11&1 \end{array}\right) \quad\overset{R_2-2R_1}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&0&-8&0\\ 2&-3&-11&1 \end{array}\right) \quad\overset{R_3-R_2}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&0&-8&0\\ 0&0&-16&0 \end{array}\right). $$ The second step can not be performed, since $a_{22}$ is zero and it can not be changed by switching second row with any row below: $$ \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&\color{red}0&-8&0\\ 0&\color{red}0&-16&0 \end{array}\right). $$ So lets skip the second step. The third step is again standard: $$ \left(\begin{array}{rrr|r} 2&-3&5&1\\ 0&0&-8&0\\ 0&0&-16&0 \end{array}\right) \quad\overset{16R_1+5R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 32&-48&0&16\\ 0&0&-8&0\\ 0&0&-16&0 \end{array}\right) \quad\overset{2R_2-R_3}{\longrightarrow}\quad \left(\begin{array}{rrr|r} 32&-48&0&16\\ 0&0&0&0\\ 0&0&-16&0 \end{array}\right). $$ Note, that the equation $0x_1+0x_2+0x_3=0$ holds always and hence, does not change the solution of the system. Thus, we may disregard the row and thus, have obtain the system ;#; $ \begin{align*} \left\{ \begin{aligned} 32x_1-48x_2&=16\\ -16x_3&=0\\ \end{aligned}\right.\enspace, \end{align*} $ ;#; which is equivalent to ;#; $ \begin{align*} \left\{ \begin{aligned} x_1&=\frac{16}{32}+\frac{48}{32}x_2\\ x_3&=0\\ \end{aligned}\right.\enspace. \end{align*} $ ;#; Note, that $x_3$ has a concrete value, $x_1$ has a criteria for evaluating it, and $x_2$ has no boundaries for valuations ($x_2$ is a **free variable** in this system). Clearly, there isn't just one set of values $x_1=c_1,\ x_2=c_2,\ x_3=c_3$ which satisfy the system of equations, that is, the system has no unique solution. In such a case, the solution is described as a set $$\left\{\left(\frac{16}{32}+\frac{48}{32}c_2,\ c_2,\ 0\right): c_2\in\mathbb{R}\right\}.$$ By giving $c_2$ a concrete value, for example, let us take $c_2=0$, we obtain one solution of the system: $(\frac{16}{32},0,0)$, but it is not the only solution. Actually, the system has infinitely many solutions, since we can choose infinitely many different values for $c_2$. Use the IPython function you created for simple Gaussian elimination algoritm (e.g., [[03_gaussian_elimination#simple_gaussian_func|simple_gaussian]]) for solving the system ;#; $ \begin{align*} \left\{ \begin{aligned} 2x_1-3x_2+5x_3&=1\\ 4x_1-6x_2+2x_3&=2\\ 2x_1-3x_2-11x_3&=1\\ \end{aligned}\right.\enspace. \end{align*} $ ;#; Analyse the result: what should be changed in order to obtain the final augmented matrix from above.\\ from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j in matrix M if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R def changed_gaussian_step(M, i): # function for step i in simple gaussian elimination if (i < 1 or i >len(M)): return 'No such step in this matrix' R = array(M) if any(R[range(i-1,len(R)),i-1] <> 0): # here is a change s = i while (R[i-1,i-1] == 0 and s < len(R)): R = swap_rows(R, i-1, s) s = s + 1 print R for j in range(len(R)): if (R[j,i-1] <> 0 and j <> (i - 1)): l = R[i-1,i-1] k = R[j,i-1] if cmp(l,0) <> cmp(k,0): R = add_row_multiples(R, j, i-1, abs(l), abs(k)) else : R = add_row_multiples(R, j, i-1, l, -k) print R else : print R return R ########################################################################################################### # it is possible to get the final augmented matrix with the changes def changed_simple_gaussian(M): R = array(M) for i in range(len(R)): R = changed_gaussian_step(R,i+1) # here is a change if (R == numpy.zeros(shape=(len(R),len(R.T)))).all() : print 'Error: not solvable with this method' return T = R.astype(numpy.float) for j in range(len(T)): if R[j,j] <> 0: T[j,:] = (float(1)/R[j,j])*R[j,:] # here is a change print T return T M = array([[2,-3,5,1],[4,-6,2,2],[2,-3,-11,1]]) changed_simple_gaussian(M) Now apply the changed simple Gaussian elimination to the system ;#; $ \begin{align*} \left\{ \begin{aligned} 2x_1-3x_2-x_3+5x_4&=1\\ 4x_1-6x_2-2x_3+2x_4&=2\\ 2x_1-3x_2-x_3-11x_4&=1\\ \end{aligned}\right.\enspace. \end{align*} $ ;#; Do you get the correct answer? ++Answer | The final augmented matrix is $\left(\begin{array}{rrrr|r}1&-\frac{3}{2}&-\frac{1}{2}&\frac{5}{2}&\frac{1}{2}\\0&0&0&1&0\\0&0&0&0&0\end{array}\right)$ ++ \\ As you can see, the simple Gaussian elimination algorithm does not work in such a case. Next we will generalize the simple Gaussian elimination algorithm in order to be able to solve any (solvable) system. So, let us have a system of linear equation without any restrictions to $m$ or $n$ or number of solutions $$ \begin{pmatrix} a_{11}&a_{12}&\dots &a_{1n}\\ a_{21}&a_{22}&\dots &a_{2n}\\ \vdots & \vdots &\ddots&\vdots\\ a_{m1}&a_{m2}&\dots &a_{mn} \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ \vdots\\ x_m \end{pmatrix}= \begin{pmatrix} b_1\\ b_2\\ \vdots\\ b_m \end{pmatrix}. $$ The goal of the general Gaussian elimination algorithm is to give to the augmented matrix of the linear system the following form\\ \\ {{ :linear_algebra:gaussian_uus.png |}} \\ by using elementary row operations and leaving out rows consisting only of zeros ($r\leq m$). In the obtained matrix, if $r=n$ and $a_{ii}=1$ for every $i=1,\dots,n$, then **the system has a unique solution**. The solution is the last column of the augmented matrix, see [[03_gaussian_elimination#simple_gaussian|simple Gaussian elimination algorithm]]. If $r Which of the following augmented matrices are the results of general Gaussian elimination algorithm. If some matrix is the result of general Gaussian elimination algorithm, then decide weather it has unique solution or not. |$\left(\begin{array}{rrrrr|r}1&0&0&5&8&0\\0&0&1&2&7&0\\0&0&0&0&0&1\end{array}\right)$ | $\left(\begin{array}{rrrrr|r}1&0&6&5&8&2\\0&0&2&2&7&3\end{array}\right)$ | |++Answer | This one is. However, the system has no solution. ++ | ++Answer | This one is not. The entry $a_{23}=2$ but in case of general Gaussian elimination algorithm there should be $1$. By multiplying the second row by $\frac{1}{2}$ we obtain the desired final form of Gaussian elimination. The system has more then one solution. ++ | |$\left(\begin{array}{rrr|r}1&0&0&3\\0&1&0&-2\\0&0&1&10\\0&0&0&0\end{array}\right)$ | $\left(\begin{array}{rrrrr|r}1&0&0&0&0&3\\0&0&9&1&0&-3\\0&1&-2&0&0&2\\0&0&0&0&1&5\end{array}\right)$ | | ++Answer | No, since the last row consists of only zeros. If we drop the last row, then the augmented matrix is in the needed form. The system has precisely one solution. ++ | ++Answer | No. The 3rd and 2nd row should be switched and instead of $9$ there should be $1$ and instead of $-2$ there should be $0$. ++ | Note that this time we do not move within the augmented matrix along the diagonal elements a_ii. However, the simple Gaussian elimination is a special case of the general Gaussian elimination: just assume, that all the blocks of arbitrary numbers are with 0 width. The general Gaussian elimination step:\\ row $i$ and column $j$ **if** there exists row $k$, $k\leq j$ such that $a_{kj}\neq 0$ then switch rows till $a_{ij}$ is nonzero and by using the new $a_{ij}$ zero out all the rest of the elements in column $j$; next step row $i+1$ and column $j+1$ \\ \\ **else** next step row $i$ and column $j+1$. Write a function in IPython for the general Gaussian elimination step. from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = array(M) R[i,:] = l*M[i,:] + k*M[j,:] return R ########################################################################################################### def general_gaussian_step(M, i, j): # i is row index and j is column index: i,j=0,1,... print M if (i < 0 or i > len(M)-1 or j < 0 or j > len(M.T)-2): return 'No such step in this matrix' R = array(M) # the result matrix if any(R[range(i,len(R)),j] <> 0): # if nonzero elements in column j starting from row i s = i + 1 while (R[i,j] == 0 and s < len(R)): # switched rows till the element a_ij is nonzero R = swap_rows(R, i, s) s = s + 1 print R for t in range(len(R)): # zeros the elements below and above if (R[t,j] <> 0 and t <> i): l = R[i,j] k = R[t,j] if cmp(l,0) <> cmp(k,0): R = add_row_multiples(R, t, i, abs(l), abs(k)) else : R = add_row_multiples(R, t, i, l, -k) print R return R \\ Write a function in IPython for solving a system of linear equations (no restriction to the system here) with the Gaussian elimination algorithm. from numpy import array def swap_rows(M, i, j): # function for swapping rows i and j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = M R[[i,j],:] = R[[j,i],:] return R def add_row_multiples(M, i, j, l, k): # function for adding l times row i and k times row j if (i > len(M) or j > len(M)): return 'The rows do not exist' else: R = M R[i,:] = l*M[i,:] + k*M[j,:] return R def general_gaussian_step(M, i, j): # function for general gaussian step, row i and col j if (i < 0 or i > len(M)-1 or j < 0 or j > len(M.T)-2): return 'No such step in this matrix' R = array(M) if any(R[range(i,len(R)),j] <> 0): s = i + 1 while (R[i,j] == 0 and s < len(R)): R = swap_rows(R, i, s) s = s + 1 print R for t in range(len(R)): if (R[t,j] <> 0 and t <> i): l = R[i,j] k = R[t,j] if cmp(l,0) <> cmp(k,0): R = add_row_multiples(R, t, i, abs(l), abs(k)) else : R = add_row_multiples(R, t, i, l, -k) print R return R ############################################################################################################# def general_gaussian(M): R = array(M) # the result matrix i = 0 # row index j = 0 # column index while (j < len(M.T)-1 and i < len(R)): R = general_gaussian_step(R,i,j) # general gaussian step for row i and column j if R[i,j] <> 0: i = i+1 # criteria for changing row index j = j+1 # in every step changes column index rowind = array(range(len(R))) # vector of row indexes in R R = R[(rowind < i) | (R[:,len(R.T)-1] <> 0),:] # leaving out rows consisting only of zeros print R if (all(R[len(R)-1,range(len(R.T)-1)] == 0)): # enough to check the last row for contradiction print 'The system has no solution' return R = R.astype(numpy.float) t = 0 # row index s = 0 # column index while (s < len(R.T) and t < len(R)): if R[t,s] <> 0: # if the leading element is not zero R[t,:] = (float(1)/R[t,s])*R[t,:] # the leading el. one by mult. the row with the inverse print R t = t+1 s = s+1 if len(R) == len(R.T)-1: # precisely one solution, return the last column print R[:,len(R.T)-1] return else : print 'The system has more then one solution' S = [0]*(len(R.T)-1) # we take the free variables equal to 0 for m in range(len(R)): for n in range(len(R.T)): if R[m,n] == 1: #if leading el., the variable equals to constant term S[n] = R[m,len(R.T)-1] break print 'One of the solutions is', S return \\ Find all solutions of the system $A\mathbf{x}=\mathbf{b}$ where $$A=\begin{pmatrix}1&1&0\\0&0&1\\0&0&2\end{pmatrix}$$ and $\mathbf{b}=\begin{pmatrix}1&2&4\end{pmatrix}^T$ over $\mathbb{Z}_3$. [[linear_algebra:04_invertible_matrices]]