# 381 - A17 - Linear Systems with Gaussian Elimination

## 247 days ago by Professor381

John Travis

Mississippi College

Solving Linear Systems

Linear systems of equations arise from many applications...often with a huge number of equations and a corresponding huge number of unknowns. In general, a system of m linear equations in the n unknowns $x_1, x_2, ... , x_n$ is of the form

$a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + ... + a_{1n}x_n = b_1$

$a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + ... + a_{2n}x_n = b_2$

$a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + ... + a_{3n}x_n = b_3$

...

$a_{m1}x_1 + a_{m2}x_2 + a_{m3}x_3 + ... + a_{mn}x_n = b_m$

and by collecting the given coefficients $a_{kj}$ into a mxn matrix A

$A = \begin{bmatrix} a_{11} & a_{12} & ... & a_{1n}\\ a_{21} & a_{22} & ... & a_{2n}\\ ... & ... & ... & ... \\ a_{m1} & a_{m2} & ... & a_{mn} \end{bmatrix}$

and the given numbers $b_k$ into a mx1 vector b gives the matrix equation

$Ax=b$.

Since solving a linear equation is a finite process, there should be no "trucation error" in general when devising algorithms to determine a solution.  So, the major source of error will arise up due to the use of floating point arithmetic...that is "roundoff error".

Note, in this worksheet we will be running in Matlab/Octave mode so that all calculations are done using floating point arithmetic and so that any resulting errors in the answers generated will be due to the use of floating point arithmetic.  We will also concernn ourselves with only solving a square linear system written in matrix form as $Ax = b$  where $A$ is a nxn matrix that is given and $b$ is a similar nx1 vector that has also been provided.  Our goal will be to solve for the unknown vector $x$.

# Assume A is a "diagonal" matrix A = zeros(5,5); # to create a 5x5 matrix consisting of all zeros b = [12; 11; 10; 9; 8]; A(1,1) = 2; # and now adjust so that the diagonal entries are nice values A(2,2) = -3; A(3,3) = 4; A(4,4) = -5; A(5,5) = 6; format short A b
 A = 2 0 0 0 0 0 -3 0 0 0 0 0 4 0 0 0 0 0 -5 0 0 0 0 0 6 b = 12 11 10 9 8 A = 2 0 0 0 0 0 -3 0 0 0 0 0 4 0 0 0 0 0 -5 0 0 0 0 0 6 b = 12 11 10 9 8

Notice, for this diagonal matrix the first row can be written out as the linear equation

$2x_1 + 0x_2 + 0x_3 + 0x_4 + 0x_5 = 12$ or $2x_1 = 12$

and since there is only one variable to deal with then solving for $x_1 = 6$ is trivial. The second row can be written out as the linear equation

$0x_1 - 3x_2 + 0x_3 + 0x_4 + 0x_5 = 11$ and thus $x_2 = -11/3$.

Write out the actual equations for the other rows.

Therefore, for diagonal matrices it is easy to develop an algorithm for generating the solution of $Ax=b$...

clear A = zeros(5,5); b = [12; 11; 10; 9; 8]; A(1,1) = 2; A(2,2) = -3; A(3,3) = 4; A(4,4) = -5; A(5,5) = 6; n = size(A)(1); for k=1:n if A(k,k)!=0 x(k) = b(k)/A(k,k); else disp("No Solution") end end x = x' b-A*x # to check the answer, this must be zeros
 x = 6 -3.66667 2.5 -1.8 1.33333 ans = 0 0 0 0 0 x = 6 -3.66667 2.5 -1.8 1.33333 ans = 0 0 0 0 0

Go ahead and create your own diagonal system of size 8x8 and solve using the code above.

# Put your code here. Perhaps copying and pasting and then editing the code above.

Ok. So, diagonal systems are trivial to solve.  But what about more complex systems?

Consider a "lower triangular" system which might be something of the form

$2x_1 = 4$

$-x_1 + 3x_2 = 9$

$5x_1 - 4x_2 + 6x_3 = 0$

This gives a matrix form $Ax=b$ with

$A = \begin{bmatrix}2 & 0 & 0 \\-1 & 3 & 0 \\5 & -4 & 6\end{bmatrix}, x = \begin{bmatrix}x_1 \\x_2 \\x_3\end{bmatrix}, b=\begin{bmatrix}4 \\9 \\0\end{bmatrix}$

Notice that when written out, the first equation has only one variable, the second has only two, and the third has only three.  For larger lower triangular matrices this pattern would continue.  This leads us to succesively solve such a system by getting the value for $x_1$ and then using it to reduce the number of variables left in all the other equations by one.  Indeed, the second equation would now only have one unknown variable.  Find it and broadcast to the remaining equations and repeat giving only one new variable found for each new equation considered. One can easily follow this patter to show that the final solution to this lower triangular system is $x_1 = 2, x_2 = 11/3, x_3 = 7/9$.

This is known as "forward substitution".

For the following matrix, write out the actual equations proscribed using variables $x_1, x_2, x_3, x_4$.

# Forward Substitution for a Lower Triangular System clear A = [2 0 0 0; -1 3 0 0; 3 4 -2 0; 0 1 9 5] b = [8; -2; 7; 4] n = size(b)(1); # Now forward substitution x(1) = b(1)/A(1,1); # the first equation for k = 2:n # each successive equation temp = 0; for j = 1:k-1 # accumulate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end x = x' # And to check our result, let's see if Ax actually equals b by computing b-Ax which should be zero. b - A*x
 A = 2 0 0 0 -1 3 0 0 3 4 -2 0 0 1 9 5 b = 8 -2 7 4 x = 4 0.666667 3.83333 -6.23333 ans = 0 0 0 0 A = 2 0 0 0 -1 3 0 0 3 4 -2 0 0 1 9 5 b = 8 -2 7 4 x = 4 0.666667 3.83333 -6.23333 ans = 0 0 0 0

Go ahead and create your own lower triangular system of size 6x6 and solve by modifying the code above.

# Your 6x6 lower triangular system and the modified code

Similarly, one can have a system with an "upper triangular matrix" that only has non-zero entries on and above the main diagonal.  We can solve such a system now by starting with the last equation and working our way back up through the other equations.  The last equation will only have one variable, the second to last will only have two, etc.  Doing this approach is called "backward solving" and is very similar to forward solving.

# Backward Substitution for an Upper Triangular System clear A = [2 -1 3 0; 0 3 4 1; 0 0 -2 9; 0 0 0 5] b = [8; -2; 7; 4] n = size(b)(1); # Now backward substitution x(n) = b(n)/A(n,n); # the last equation for k = n-1:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end x = x' # And to check our result, let's see if Ax actually equals b by computing b-Ax which should be zero. b - A*x
 A = 2 -1 3 0 0 3 4 1 0 0 -2 9 0 0 0 5 b = 8 -2 7 4 x = 3.31667 -1.06667 0.1 0.8 ans = 0 -2.22045e-16 0 0 A = 2 -1 3 0 0 3 4 1 0 0 -2 9 0 0 0 5 b = 8 -2 7 4 x = 3.31667 -1.06667 0.1 0.8 ans = 0 -2.22045e-16 0 0

Go ahead and create your own upper triangular system of size 7x7 and solve by modifying the code above.

# Your 7x7 upper triangular system and the modified code

Next, let's consider a general nxn system which doesn not have any special structure (such as diagonal or triangular).  For these, we use a process known as Gaussian Elimination that utilizes the following general linear algebra properties:

1. Any equation can be multiplied by a non-zero number and the solution to the resulting system does not change.
2. Any two equations can be swapped and the solution of the resulting system does not change.
3. Any multiple of one equation can be added to any other equation and then replace the second equation by this sum and the solution of the resulting system does not change.

Using these repeatedly in a manner that aims for creating zero entries below the diagonal is the essence of Gaussian Elimination.  Once this is done, one will have a new upper triangular system that has the same solution as the original system.  Finish by backsolving.

First, consider a situation in which there is no need to swap any rows...

# Gaussian Elimination (without Pivoting) clear A = [2 -1 3 0; 4 3 4 1; -1 1 -2 -3; 5 0 0 4] b = [6; 9; -12; 37] n = size(b)(1); # Let's save the originals... Aorig = A; borig = b; # Make the (1,1) position equal to 1 by dividing that row. pivot = A(1,1); A(1,:) = A(1,:)/pivot; b(1) = b(1)/pivot; disp("After getting the first pivot:") A b # Make zeros below the (1,1) position by multiplying row 1 by the negative of the existing values # in the kth row below that pivot and adding that to the existing k-th row and replace for k=2:n multiplier = -A(k,1); A(k,:) = A(k,:) + multiplier*A(1,:); # kth row, all columns b(k) = b(k) + multiplier*b(1); end disp("After zeroing out the rest of column one:") A b # Make the (2,2) position equal to 1 by dividing that row. pivot = A(2,2); A(2,:) = A(2,:)/pivot; b(2) = b(2)/pivot; disp("After getting the second pivot:") A b # Make zeros below the (2,2) position by multiplying row 2 by the negative of the existing values # in the kth row below that pivot and adding that to the existing k-th row and replace for k=3:n multiplier = -A(k,2); A(k,:) = A(k,:) + multiplier*A(2,:); # kth row, all columns b(k) = b(k) + multiplier*b(2); end disp("After zeroing out the rest of column two:") A b # Make the (3,3) position equal to 1 by dividing that row. pivot = A(3,3); A(3,:) = A(3,:)/pivot; b(3) = b(3)/pivot; disp("After getting the third pivot:") A b # Make zeros below the (3,3) position by multiplying row 3 by the negative of the existing values # in the kth row below that pivot and adding that to the existing k-th row and replace k=4 # this is the final "loop" with only one iteration multiplier = -A(k,3); A(k,:) = A(k,:) + multiplier*A(3,:); # kth row, all columns b(k) = b(k) + multiplier*b(3); disp("After zeroing out the rest of column three:") format short A b # Now, we are left with an upper triangular system which can be solved as above using backward substitution x(n) = b(n)/A(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end disp('With the final answer for Ax=b as x equals:') x = x' # And if you want to check it all out # borig - Aorig*x
 A = 2 -1 3 0 4 3 4 1 -1 1 -2 -3 5 0 0 4 b = 6 9 -12 37 After getting the first pivot: A = 1.00000 -0.50000 1.50000 0.00000 4.00000 3.00000 4.00000 1.00000 -1.00000 1.00000 -2.00000 -3.00000 5.00000 0.00000 0.00000 4.00000 b = 3 9 -12 37 After zeroing out the rest of column one: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 5.00000 -2.00000 1.00000 0.00000 0.50000 -0.50000 -3.00000 0.00000 2.50000 -7.50000 4.00000 b = 3 -3 -9 22 After getting the second pivot: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 0.00000 0.50000 -0.50000 -3.00000 0.00000 2.50000 -7.50000 4.00000 b = 3.00000 -0.60000 -9.00000 22.00000 After zeroing out the rest of column two: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 0.00000 0.00000 -0.30000 -3.10000 0.00000 0.00000 -6.50000 3.50000 b = 3.00000 -0.60000 -8.70000 23.50000 After getting the third pivot: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 -0.00000 -0.00000 1.00000 10.33333 0.00000 0.00000 -6.50000 3.50000 b = 3.00000 -0.60000 29.00000 23.50000 k = 4 After zeroing out the rest of column three: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 -0.00000 -0.00000 1.00000 10.33333 0.00000 0.00000 0.00000 70.66667 b = 3.00000 -0.60000 29.00000 212.00000 With the final answer for Ax=b as x equals: x = 5 -2 -2 3 A = 2 -1 3 0 4 3 4 1 -1 1 -2 -3 5 0 0 4 b = 6 9 -12 37 After getting the first pivot: A = 1.00000 -0.50000 1.50000 0.00000 4.00000 3.00000 4.00000 1.00000 -1.00000 1.00000 -2.00000 -3.00000 5.00000 0.00000 0.00000 4.00000 b = 3 9 -12 37 After zeroing out the rest of column one: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 5.00000 -2.00000 1.00000 0.00000 0.50000 -0.50000 -3.00000 0.00000 2.50000 -7.50000 4.00000 b = 3 -3 -9 22 After getting the second pivot: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 0.00000 0.50000 -0.50000 -3.00000 0.00000 2.50000 -7.50000 4.00000 b = 3.00000 -0.60000 -9.00000 22.00000 After zeroing out the rest of column two: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 0.00000 0.00000 -0.30000 -3.10000 0.00000 0.00000 -6.50000 3.50000 b = 3.00000 -0.60000 -8.70000 23.50000 After getting the third pivot: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 -0.00000 -0.00000 1.00000 10.33333 0.00000 0.00000 -6.50000 3.50000 b = 3.00000 -0.60000 29.00000 23.50000 k = 4 After zeroing out the rest of column three: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 -0.00000 -0.00000 1.00000 10.33333 0.00000 0.00000 0.00000 70.66667 b = 3.00000 -0.60000 29.00000 212.00000 With the final answer for Ax=b as x equals: x = 5 -2 -2 3
# Now to put the elimination part all together into one big loop that can be generalized to any size... # Gaussian Elimination (without Pivoting) clear A = [2 -1 3 0; 4 3 4 1; -1 1 -2 -3; 5 0 0 4] b = [6; 9; -12; 37] n = size(b)(1); for step = 1:n-1 pivot = A(step,step); A(step,:) = A(step,:)/pivot; b(step) = b(step)/pivot; for k=step+1:n multiplier = -A(k,step); A(k,:) = A(k,:) + multiplier*A(step,:); # kth row, all columns b(k) = b(k) + multiplier*b(step); end end disp('The result after using Gaussian Elimination:') A b # And again, backward substitution x(n) = b(n)/A(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end disp('With the final answer for Ax=b as x equals:') x = x'
 A = 2 -1 3 0 4 3 4 1 -1 1 -2 -3 5 0 0 4 b = 6 9 -12 37 The result after using Gaussian Elimination: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 -0.00000 -0.00000 1.00000 10.33333 0.00000 0.00000 0.00000 70.66667 b = 3.00000 -0.60000 29.00000 212.00000 With the final answer for Ax=b as x equals: x = 5 -2 -2 3 A = 2 -1 3 0 4 3 4 1 -1 1 -2 -3 5 0 0 4 b = 6 9 -12 37 The result after using Gaussian Elimination: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -0.40000 0.20000 -0.00000 -0.00000 1.00000 10.33333 0.00000 0.00000 0.00000 70.66667 b = 3.00000 -0.60000 29.00000 212.00000 With the final answer for Ax=b as x equals: x = 5 -2 -2 3

For the fun of it let's create a bigger A and use a given x to determine the vector b.  Then, pretend you don't know x and solve $Ax=b$.

clear A = floor(10*rand(6,6)+1) xexact = floor(8*rand(6,1)-3) b = A*xexact # Gaussian Elimination (without Pivoting) n = size(b)(1); for step = 1:n-1 pivot = A(step,step); A(step,:) = A(step,:)/pivot; b(step) = b(step)/pivot; for k=step+1:n multiplier = -A(k,step); A(k,:) = A(k,:) + multiplier*A(step,:); # kth row, all columns b(k) = b(k) + multiplier*b(step); end end disp('The result after using Gaussian Elimination:') A b # And again, backward substitution x(n) = b(n)/A(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end disp('With the final answer for Ax=b as x equals:') x' # And to see how close we got... error = x'-xexact
 A = 6 5 7 9 7 3 10 9 3 8 3 1 1 10 2 9 7 9 3 7 3 4 10 7 1 4 1 8 1 3 4 1 8 9 5 6 xexact = -1 3 0 3 1 -1 b = 40 43 54 33 33 25 The result after using Gaussian Elimination: A = 1.00000 0.83333 1.16667 1.50000 1.16667 0.50000 0.00000 1.00000 -13.00000 -10.50000 -13.00000 -6.00000 0.00000 0.00000 1.00000 0.86458 1.04167 0.52917 -0.00000 -0.00000 -0.00000 1.00000 -1.34969 -0.53252 0.00000 0.00000 0.00000 0.00000 1.00000 0.51123 0.00000 0.00000 0.00000 0.00000 0.00000 4.95569 b = 6.66667 -35.50000 3.10625 2.18282 0.48877 -4.95569 With the final answer for Ax=b as x equals: ans = -1.0000e+00 3.0000e+00 -3.5527e-15 3.0000e+00 1.0000e+00 -1.0000e+00 error = -8.8818e-16 0.0000e+00 -3.5527e-15 8.8818e-16 3.7748e-15 -1.9984e-15 A = 6 5 7 9 7 3 10 9 3 8 3 1 1 10 2 9 7 9 3 7 3 4 10 7 1 4 1 8 1 3 4 1 8 9 5 6 xexact = -1 3 0 3 1 -1 b = 40 43 54 33 33 25 The result after using Gaussian Elimination: A = 1.00000 0.83333 1.16667 1.50000 1.16667 0.50000 0.00000 1.00000 -13.00000 -10.50000 -13.00000 -6.00000 0.00000 0.00000 1.00000 0.86458 1.04167 0.52917 -0.00000 -0.00000 -0.00000 1.00000 -1.34969 -0.53252 0.00000 0.00000 0.00000 0.00000 1.00000 0.51123 0.00000 0.00000 0.00000 0.00000 0.00000 4.95569 b = 6.66667 -35.50000 3.10625 2.18282 0.48877 -4.95569 With the final answer for Ax=b as x equals: ans = -1.0000e+00 3.0000e+00 -3.5527e-15 3.0000e+00 1.0000e+00 -1.0000e+00 error = -8.8818e-16 0.0000e+00 -3.5527e-15 8.8818e-16 3.7748e-15 -1.9984e-15

Notice that we haven't used one of the allowable row operations...that is, we haven't ever swapped rows.  Sometimes this might be necessary or at least "adivsable".  See what happens when using the matrix provided below.  What is the problem and why does this fail?

# Now to put the elimination part all together into one big loop... # Gaussian Elimination (without Pivoting) clear A = [2 -1 3 0; 0 0 4 1; -1 1 -2 -3; 5 0 0 4] b = [6; 9; -12; 37] n = size(b)(1); for step = 1:n-1 pivot = A(step,step); A(step,:) = A(step,:)/pivot; b(step) = b(step)/pivot; for k=step+1:n multiplier = -A(k,step); A(k,:) = A(k,:) + multiplier*A(step,:); # kth row, all columns b(k) = b(k) + multiplier*b(step); end end disp('The result after using Gaussian Elimination:') A b # And again, backward substitution x(n) = b(n)/A(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end disp('With the final answer for Ax=b as x equals:') x
 A = 2 -1 3 0 0 0 4 1 -1 1 -2 -3 5 0 0 4 b = 6 9 -12 37 warning: division by zero warning: called from /home/sageserver/.sage/temp/travis-PowerEdge-T430/278042/interface/tmp27\ 8048 at line 15 column 13 warning: division by zero warning: called from /home/sageserver/.sage/temp/travis-PowerEdge-T430/278042/interface/tmp27\ 8048 at line 16 column 11 The result after using Gaussian Elimination: A = 1.00000 -0.50000 1.50000 0.00000 NaN NaN Inf Inf NaN NaN NaN NaN NaN NaN NaN NaN b = 3 Inf NaN NaN With the final answer for Ax=b as x equals: x = NaN NaN NaN NaN A = 2 -1 3 0 0 0 4 1 -1 1 -2 -3 5 0 0 4 b = 6 9 -12 37 warning: division by zero warning: called from /home/sageserver/.sage/temp/travis-PowerEdge-T430/278042/interface/tmp278048 at line 15 column 13 warning: division by zero warning: called from /home/sageserver/.sage/temp/travis-PowerEdge-T430/278042/interface/tmp278048 at line 16 column 11 The result after using Gaussian Elimination: A = 1.00000 -0.50000 1.50000 0.00000 NaN NaN Inf Inf NaN NaN NaN NaN NaN NaN NaN NaN b = 3 Inf NaN NaN With the final answer for Ax=b as x equals: x = NaN NaN NaN NaN

This matrix did not perform so well due to the fact that one of the pivot elements turned out to be zero. Since we have to divide by that element to proceed, the algorithm bombs.

If we proceed blindly through Gaussian Elimination without considering swapping rows then it is possible that we can end up with a pivot element that equals 0 for which dividing in the algorithm is not possible.  However, if we can swap the row with a zero in the desired pivot spot with another row below that has a non-zero then we can be good to go.

Using the problem presented above, let's swap rows 2 and 3 (from the beginning) and see if that works.

To see what makes this different, write out the equations from the matrix problem given above and then write out the equations for the matrix problem given below.

# Now to put the elimination part all together into one big loop... # Gaussian Elimination (without Pivoting) clear A = [2 -1 3 0; -1 1 -2 -3; # this was previously row 3 0 0 4 1; # this was previously row 2 5 0 0 4] b = [6;-12; 9; 37] # we also swapped the 9 and -12 since they correspond to rows 2 and 3 as well n = size(b)(1); for step = 1:n-1 pivot = A(step,step); A(step,:) = A(step,:)/pivot; b(step) = b(step)/pivot; for k=step+1:n multiplier = -A(k,step); A(k,:) = A(k,:) + multiplier*A(step,:); # kth row, all columns b(k) = b(k) + multiplier*b(step); end end disp('The result after using Gaussian Elimination:') A b # And again, backward substitution x(n) = b(n)/A(n,n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(k,j)*x(j); end x(k) = (b(k)-temp)/A(k,k); end disp('With the final answer for Ax=b as x equals:') x = x'
 A = 2 -1 3 0 -1 1 -2 -3 0 0 4 1 5 0 0 4 b = 6 -12 9 37 The result after using Gaussian Elimination: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -1.00000 -6.00000 0.00000 0.00000 1.00000 0.25000 0.00000 0.00000 0.00000 20.25000 b = 3.0000 -18.0000 2.2500 78.2500 With the final answer for Ax=b as x equals: x = 4.3086 6.4691 1.2840 3.8642 A = 2 -1 3 0 -1 1 -2 -3 0 0 4 1 5 0 0 4 b = 6 -12 9 37 The result after using Gaussian Elimination: A = 1.00000 -0.50000 1.50000 0.00000 0.00000 1.00000 -1.00000 -6.00000 0.00000 0.00000 1.00000 0.25000 0.00000 0.00000 0.00000 20.25000 b = 3.0000 -18.0000 2.2500 78.2500 With the final answer for Ax=b as x equals: x = 4.3086 6.4691 1.2840 3.8642

And so the Gaussian Elimination algorithm worked if we were to know ahead of time where the zeros might pop up.

In reality and especially for larger systems, there is no way that we might know when this zero pivot event might occur.  Indeed, the algorithm itself should determine by itself when such a swap is necessary.  To do this will in theory require rows to be swapped dynamically as the algorithm progresses and doing such a swap is tedious.  An easier way to do swaps is by using pointers to the rows and simply swapping pointers instead of entire rows.

Initially, set up pointers so that p(1) = 1, p(2) = 2, p(3) = 3, ... , p(n) = n so that the kth row points to the physical kth row.  Then, to swap (say) rows 4 and 7 would require setting r(4)=7 and r(7)=4.  Watch how this works in the Gaussian Elimination with Pivoting algorithm below.

Note, the final output for the pointers shows where the kth row ended up logically.

clear A = floor(10*rand(7,7)-4); A(2,1) = 0; A(2,2) = 0; # let's make a pivot needing swapping in row 2 A(4,1) = 0; A(4,2) = 0; A(4,3) = 0;A(4,4) = 0; # and in row 4 A xexact = floor(8*rand(7,1)-3) b = A*xexact # Gaussian Elimination (with Partial Pivoting) n = size(b)(1); for k = 1:n p(k) = k; # initialize the pointers end for step=1:n-1 temp = A(p(step),step); if temp == 0 nonpivot = step-1; do nonpivot = nonpivot + 1; until ( A(p(nonpivot),step) ~=0) holdp = p(step); # swap pointers which means step row swapped with nonpivot row p(step) = p(nonpivot); p(nonpivot) = holdp; end pivot = A(p(step),step); for k=step+1:n multiplier = -A(p(k),step); A(p(k),:) = A(p(k),:) + multiplier*A(p(step),:)/pivot; # kth row, all columns b(p(k)) = b(p(k)) + multiplier*b(p(step))/pivot; end end # And again, backward substitution x(n) = b(p(n))/A(p(n),n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(p(k),j)*x(j); end x(k) = (b(p(k))-temp)/A(p(k),k); end disp('With the final answer for Ax=b as x equals:') x' # And to see how close we got... error = x'-xexact disp('The ending state of the pointers:') p
 A = 3 -2 1 -4 0 5 -3 0 0 -1 4 -4 5 2 3 -1 3 0 0 -1 -4 0 0 0 0 4 -3 5 -3 2 4 0 5 3 1 2 -2 0 -1 5 5 3 3 -4 0 5 0 -4 1 xexact = 2 -3 0 4 1 4 1 b = 13 34 1 -3 6 34 23 With the final answer for Ax=b as x equals: ans = 2 -3 -0 4 1 4 1 error = -3.55271e-15 -7.10543e-15 -0 1.77636e-15 1.55431e-15 0 -1.22125e-15 The ending state of the pointers: p = 1 3 2 5 4 6 7 A = 3 -2 1 -4 0 5 -3 0 0 -1 4 -4 5 2 3 -1 3 0 0 -1 -4 0 0 0 0 4 -3 5 -3 2 4 0 5 3 1 2 -2 0 -1 5 5 3 3 -4 0 5 0 -4 1 xexact = 2 -3 0 4 1 4 1 b = 13 34 1 -3 6 34 23 With the final answer for Ax=b as x equals: ans = 2 -3 -0 4 1 4 1 error = -3.55271e-15 -7.10543e-15 -0 1.77636e-15 1.55431e-15 0 -1.22125e-15 The ending state of the pointers: p = 1 3 2 5 4 6 7

Notice that when a zero is found in the (logically) next diagonal position, we swap with the first row below that row that has a non-zero in that same column.

Gaussian Elimination with Scaled Partial Pivoting looks in the step-th column at ALL the rows below the step-th and to minimize roundoff selects the swap row to be the one with the largest magnitude term in that column.

Revise the algorithm to find the largest in absolute value element in the step-th column below the step-th and swap regardless if one with larger magnitude is found.  No need then to look for zeros.

# Gaussian Elimination with Scaled Partial Pivoting # Hint Snippet: point = step; temp = abs(A(p(point),step)); for j=step+1:n if abs(A(p(j),step))>temp temp = abs(A(p(j),step)) point = j end holdp = p(step); # swap pointers which means step row swapped with nonpivot row p(step) = p(point); p(point) = holdp; end
# Here is the complete scaled partial pivoting Gaussian Elimination. clear A = floor(10*rand(7,7)-4); A(2,1) = 0; A(2,2) = 0; # let's make a pivot needing swapping in row 2 A(4,1) = 0; A(4,2) = 0; A(4,3) = 0;A(4,4) = 0; # and in row 4 A xexact = floor(8*rand(7,1)-3) b = A*xexact n = size(b)(1); for k = 1:n p(k) = k; # initialize the pointers end for step=1:n-1 point = step; temp = abs(A(p(point),step)); for j=step+1:n if abs(A(p(j),step))>temp temp = abs(A(p(j),step)) point = j end if temp == 0 disp("Matrix is non-invertible. Need to address issues...") end holdp = p(step); # swap pointers which means step row swapped with nonpivot row p(step) = p(point); p(point) = holdp; end pivot = A(p(step),step); for k=step+1:n multiplier = -A(p(k),step); A(p(k),:) = A(p(k),:) + multiplier*A(p(step),:)/pivot; # kth row, all columns b(p(k)) = b(p(k)) + multiplier*b(p(step))/pivot; end step,A end # And again, backward substitution x(n) = b(p(n))/A(p(n),n); # the last equation for k = n:-1:1 # each successive equation temp = 0; for j = n:-1:k+1 # accululate the previously determined terms for row k temp += A(p(k),j)*x(j); end x(k) = (b(p(k))-temp)/A(p(k),k); end disp('With the final answer for Ax=b as x equals:') x' # And to see how close we got... error = x'-xexact disp('The ending state of the pointers:') p
 WARNING: Output truncated! full_output.txt A = -2 4 -2 -4 4 -1 4 0 0 -1 5 2 5 4 2 3 -1 0 2 -3 1 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 -4 5 2 0 -2 5 -4 -4 0 -2 -2 -2 1 -1 xexact = -2 -1 1 0 0 3 2 b = 3 22 -15 14 -21 12 7 temp = 5 point = 5 step = 1 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 1 -0.6 0.8 0 -1.8 0.2 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 9 1.2 -1.6 2 2.6 -2.4 0 4 -2.8 -3.6 2 -1.4 0.6 temp = 1 point = 3 temp = 6 point = 5 temp = 9 point = 6 step = 2 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 -0.2 1.6 -1 -1.43333 -0.6 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 4.8 5.6 -7 5.9 -9.6 0 0 -1.2 -0.4 -2 0.0666667 -2.6 ... A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 0 0.6 -1.4 -2.43333 -1.4 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 0 -3.55271e-15 71.6667 149.944 78.6667 0 0 0 0 -19.3333 -31.8889 -22.3333 temp = 71.6667 point = 6 step = 5 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 0 0.6 -1.4 -2.43333 -1.4 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 0 -3.55271e-15 0 185.778 150.333 0 0 0 0 0 -41.5556 -41.6667 step = 6 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 0 0.6 -1.4 -2.43333 -1.4 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 0 -3.55271e-15 0 185.778 150.333 0 0 0 -7.94686e-16 0 0 -8.03947 With the final answer for Ax=b as x equals: ans = -2 -1 1 2.96059e-15 -0 3 2 error = -4.44089e-15 8.10463e-15 1.42109e-14 2.96059e-15 -0 0 0 The ending state of the pointers: p = 5 1 2 3 4 6 7 WARNING: Output truncated! full_output.txt A = -2 4 -2 -4 4 -1 4 0 0 -1 5 2 5 4 2 3 -1 0 2 -3 1 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 -4 5 2 0 -2 5 -4 -4 0 -2 -2 -2 1 -1 xexact = -2 -1 1 0 0 3 2 b = 3 22 -15 14 -21 12 7 temp = 5 point = 5 step = 1 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 1 -0.6 0.8 0 -1.8 0.2 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 9 1.2 -1.6 2 2.6 -2.4 0 4 -2.8 -3.6 2 -1.4 0.6 temp = 1 point = 3 temp = 6 point = 5 temp = 9 point = 6 step = 2 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 -0.2 1.6 -1 -1.43333 -0.6 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 4.8 5.6 -7 5.9 -9.6 0 0 -1.2 -0.4 -2 0.0666667 -2.6 ... A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 0 0.6 -1.4 -2.43333 -1.4 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 0 -3.55271e-15 71.6667 149.944 78.6667 0 0 0 0 -19.3333 -31.8889 -22.3333 temp = 71.6667 point = 6 step = 5 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 0 0.6 -1.4 -2.43333 -1.4 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 0 -3.55271e-15 0 185.778 150.333 0 0 0 0 0 -41.5556 -41.6667 step = 6 A = 0 6 -2.4 -4.8 6 -2.2 4.8 0 0 -1 5 2 5 4 0 0 0 0.6 -1.4 -2.43333 -1.4 0 0 0 0 -4 2 4 5 5 -1 -2 5 -3 2 0 0 0 -3.55271e-15 0 185.778 150.333 0 0 0 -7.94686e-16 0 0 -8.03947 With the final answer for Ax=b as x equals: ans = -2 -1 1 2.96059e-15 -0 3 2 error = -4.44089e-15 8.10463e-15 1.42109e-14 2.96059e-15 -0 0 0 The ending state of the pointers: p = 5 1 2 3 4 6 7