# Difference between revisions of "Assignment 2 - 2010 - Solution/Bonus question"

 ← Question 5 (previous step) Back to all questions

# Bonus question 

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

# 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]$$.

 ← Question 5 (previous step) Back to all questions