User:Kevindunn

Bonus question [1]

Use the pseudo code developed in the course notes to write a MATLAB or Python function that implements Gauss elimination, without pivoting. The function should take $A$ and $b$ as inputs, and return vector $x$.

Use your function to solve this set of equations from question 1, and print the update $A$ matrix and $b$ vector after every pivot step.

Solution

A function that implements the Gauss elimination without pivoting is provided below.

MATLAB code

Python code

function x = gauss(A,b) % This function performs the Gauss elimination without pivoting % % x = GAUSS(A, b) [n,n] = size(A); % Check for zero diagonal elements if any(diag(A)==0) error('Division by zero will occur; pivoting not supported') end % Forward elimination for row=1:n-1 for i=row+1:n factor = A(i,row) / A(row,row); for j = row:n A(i,j) = A(i,j) - factor*A(row,j); end b(i) = b(i) - factor*b(row); end A_and_b = [A b] end % Backward substitution x(n) = b(n) / A(n,n); for row = n-1:-1:1 sums = b(row); for j = row+1: n sums = sums - A(row,j) * x(j); end x(row) = sums / A(row,row); end %Now run this in the command window or an M-file: % >> A = [1 1 1 1;1 -2 -2 -2;1 4 -4 1;1 -5 -5 -3]; % >> b = [0 4 2 -4]'; % >> x = gauss(A,b)

import numpy as np def forward_elimination(A, b, n): """ Calculates the forward part of Gaussian elimination. """ for row in range(0, n-1): for i in range(row+1, n): factor = A[i,row] / A[row,row] for j in range(row, n): A[i,j] = A[i,j] - factor * A[row,j] b[i] = b[i] - factor * b[row] print('A = \n%s and b = %s' % (A,b)) return A, b def back_substitution(a, b, n): """" Does back substitution, returns the Gauss result. """ x = np.zeros((n,1)) x[n-1] = b[n-1] / a[n-1, n-1] for row in range(n-2, -1, -1): sums = b[row] for j in range(row+1, n): sums = sums - a[row,j] * x[j] x[row] = sums / a[row,row] return x def gauss(A, b): """ This function performs Gauss elimination without pivoting. """ n = A.shape[0] # Check for zero diagonal elements if any(np.diag(A)==0): raise ZeroDivisionError(('Division by zero will occur; ' 'pivoting currently not supported')) A, b = forward_elimination(A, b, n) return back_substitution(A, b, n) # Main program starts here if __name__ == '__main__': A = np.array([[1, 1, 1, 1], [1, -2, -2, -2], [1, 4, -4, 1], [1, -5, -5, -3]]) b = np.array([0, 4, 2, -4]) x = gauss(A, b) print('Gauss result is x = \n %s' % x)

Now we use this function to solve the system of equations in Question 1. The commented lines in the MATLAB code show how the function is used. The results are: $x =[1.3333, 3.1667,1.5000,-6.0000]$. Also, the $A$ and $b$ (denoted as $C=[A\quad b]$ below) are updated as:

$$\begin{split}C^1= \left[ \begin{array}{ccccc} 1 & 1 & 1 & 1 & 0\\ 0 & -3 & -3 & -3 & 4\\ 0 & 3 & -5 & 0 & 2\\ 0 & -6 & -6 & -4 & -4 \end{array} \right]\end{split}$$

$$\begin{split}C^2= \left[ \begin{array}{ccccc} 1 & 1 & 1 & 1 & 0\\ 0 & -3 & -3 & -3 & 4\\ 0 & 0 & -8 & -3 & 6\\ 0 & 0 & 0 & 2 & -12 \end{array} \right]\end{split}$$

$$\begin{split}C^3= \left[ \begin{array}{ccccc} 1 & 1 & 1 & 1 & 0\\ 0 & -3 & -3 & -3 & 4\\ 0 & 0 & -8 & -3 & 6\\ 0 & 0 & 0 & 2 & -12 \end{array} \right]\end{split}$$

The final result for $x = [1.333, 3.167, 1.5, -6]$.