Difference between revisions of "User:Kevindunn"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
m
m
Line 1: Line 1:
<div class="section" id="bonus-question-1">
<html><div class="section" id="bonus-question-1">
<h1>Bonus question [1]</h1>
<h1>Bonus question [1]</h1>
<p>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 <span class="math">\(A\)</span> and <span class="math">\(b\)</span> as inputs, and return vector <span class="math">\(x\)</span>.</p>
<p>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 <span class="math">\(A\)</span> and <span class="math">\(b\)</span> as inputs, and return vector <span class="math">\(x\)</span>.</p>
Line 142: Line 142:
</div>
</div>
</div>
</div>
</html>

Revision as of 07:33, 10 February 2017

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]\).