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:
''' Contact details '''
<div class="section" id="bonus-question-1">
<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 your function to solve this set of equations from question 1, and print the update <span class="math">\(A\)</span> matrix and <span class="math">\(b\)</span> vector after every pivot step.</p>
<div class="section" id="solution">
<h2>Solution</h2>
<p>A function that implements the Gauss elimination without pivoting is provided below.</p>
<table border="1" class="docutils">
<thead valign="bottom">
<tr><th>
MATLAB code</th><th>
Python code</th></tr>
</thead>
<td>
<div class="highlight-matlab"><div class="highlight"><pre><span class="k">function</span><span class="w"> </span>x <span class="p">=</span><span class="w"> </span><span class="nf">gauss</span><span class="p">(</span>A,b<span class="p">)</span><span class="w"></span>
<span class="c">% This function performs the Gauss elimination without pivoting</span>
<span class="c">% </span>
<span class="c">% x = GAUSS(A, b)</span>


I don't have an office on campus.  The best way to contact me is by email:
<span class="p">[</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">]</span> <span class="p">=</span> <span class="nb">size</span><span class="p">(</span><span class="n">A</span><span class="p">);</span>


* [mailto:dunnkg@mcmaster.ca dunnkg@mcmaster.ca] (preferred)
<span class="c">% Check for zero diagonal elements</span>
* [mailto:kevin.dunn@connectmv.com kevin.dunn@connectmv.com] (alternate)
<span class="k">if</span> <span class="n">any</span><span class="p">(</span><span class="nb">diag</span><span class="p">(</span><span class="n">A</span><span class="p">)</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">error</span><span class="p">(</span><span class="s">&#39;Division by zero will occur; pivoting not supported&#39;</span><span class="p">)</span>
<span class="k">end</span>


<span class="c">% Forward elimination</span>
<span class="k">for</span> <span class="n">row</span><span class="p">=</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span>
    <span class="k">for</span> <span class="nb">i</span><span class="p">=</span><span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span>
        <span class="nb">factor</span> <span class="p">=</span> <span class="n">A</span><span class="p">(</span><span class="nb">i</span><span class="p">,</span><span class="n">row</span><span class="p">)</span> <span class="o">/</span> <span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">);</span>
        <span class="k">for</span> <span class="nb">j</span> <span class="p">=</span> <span class="n">row</span><span class="p">:</span><span class="n">n</span>
            <span class="n">A</span><span class="p">(</span><span class="nb">i</span><span class="p">,</span><span class="nb">j</span><span class="p">)</span> <span class="p">=</span> <span class="n">A</span><span class="p">(</span><span class="nb">i</span><span class="p">,</span><span class="nb">j</span><span class="p">)</span> <span class="o">-</span> <span class="nb">factor</span><span class="o">*</span><span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="nb">j</span><span class="p">);</span>
        <span class="k">end</span>
        <span class="n">b</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="p">=</span> <span class="n">b</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="o">-</span> <span class="nb">factor</span><span class="o">*</span><span class="n">b</span><span class="p">(</span><span class="n">row</span><span class="p">);</span>
    <span class="k">end</span>
  <span class="n">A_and_b</span> <span class="p">=</span>  <span class="p">[</span><span class="n">A</span> <span class="n">b</span><span class="p">]</span>
<span class="k">end</span>
<span class="c">% Backward substitution</span>
<span class="n">x</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="p">=</span> <span class="n">b</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="o">/</span> <span class="n">A</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">);</span>
<span class="k">for</span> <span class="n">row</span> <span class="p">=</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span><span class="mi">1</span>
    <span class="n">sums</span> <span class="p">=</span> <span class="n">b</span><span class="p">(</span><span class="n">row</span><span class="p">);</span>
    <span class="k">for</span> <span class="nb">j</span> <span class="p">=</span> <span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span> <span class="n">n</span>
        <span class="n">sums</span> <span class="p">=</span> <span class="n">sums</span> <span class="o">-</span> <span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="nb">j</span><span class="p">)</span> <span class="o">*</span> <span class="n">x</span><span class="p">(</span><span class="nb">j</span><span class="p">);</span>
    <span class="k">end</span>
    <span class="n">x</span><span class="p">(</span><span class="n">row</span><span class="p">)</span> <span class="p">=</span> <span class="n">sums</span> <span class="o">/</span> <span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">);</span>
<span class="k">end</span>
<span class="c">%Now run this in the command window or an M-file:</span>
<span class="c">% &gt;&gt; A = [1 1 1 1;1 -2 -2 -2;1 4 -4 1;1 -5 -5 -3];</span>
<span class="c">% &gt;&gt; b = [0 4 2 -4]&#39;;</span>
<span class="c">% &gt;&gt; x = gauss(A,b)</span>
</pre></div>
</div>
</td><td>
<div class="highlight-python"><div class="highlight"><pre><span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>
<span class="k">def</span> <span class="nf">forward_elimination</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculates the forward part of Gaussian elimination.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">for</span> <span class="n">row</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
            <span class="n">factor</span> <span class="o">=</span> <span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">row</span><span class="p">]</span> <span class="o">/</span> <span class="n">A</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">]</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">row</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
                <span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">factor</span> <span class="o">*</span> <span class="n">A</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>
            <span class="n">b</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span> <span class="n">factor</span> <span class="o">*</span> <span class="n">b</span><span class="p">[</span><span class="n">row</span><span class="p">]</span>
        <span class="k">print</span><span class="p">(</span><span class="s">&#39;A = </span><span class="se">\n</span><span class="si">%s</span><span class="s"> and b = </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">b</span><span class="p">))</span>
    <span class="k">return</span> <span class="n">A</span><span class="p">,</span> <span class="n">b</span>
<span class="k">def</span> <span class="nf">back_substitution</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;&quot;</span>
<span class="sd">    Does back substitution, returns the Gauss result.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
    <span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">/</span> <span class="n">a</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
    <span class="k">for</span> <span class="n">row</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">):</span>
        <span class="n">sums</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">row</span><span class="p">]</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
            <span class="n">sums</span> <span class="o">=</span> <span class="n">sums</span> <span class="o">-</span> <span class="n">a</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">*</span> <span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
        <span class="n">x</span><span class="p">[</span><span class="n">row</span><span class="p">]</span> <span class="o">=</span> <span class="n">sums</span> <span class="o">/</span> <span class="n">a</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">]</span>
    <span class="k">return</span> <span class="n">x</span>
<span class="k">def</span> <span class="nf">gauss</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    This function performs Gauss elimination without pivoting.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">n</span> <span class="o">=</span> <span class="n">A</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="c"># Check for zero diagonal elements</span>
    <span class="k">if</span> <span class="nb">any</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">A</span><span class="p">)</span><span class="o">==</span><span class="mi">0</span><span class="p">):</span>
        <span class="k">raise</span> <span class="ne">ZeroDivisionError</span><span class="p">((</span><span class="s">&#39;Division by zero will occur; &#39;</span>
                                  <span class="s">&#39;pivoting currently not supported&#39;</span><span class="p">))</span>
    <span class="n">A</span><span class="p">,</span> <span class="n">b</span> <span class="o">=</span> <span class="n">forward_elimination</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">back_substitution</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span>
<span class="c"># Main program starts here</span>
<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&#39;__main__&#39;</span><span class="p">:</span>
    <span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span>  <span class="mi">1</span><span class="p">,</span>  <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                  <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span>  <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">],</span>
                  <span class="p">[</span><span class="mi">1</span><span class="p">,</span>  <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                  <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">5</span><span class="p">,</span>  <span class="o">-</span><span class="mi">5</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">]])</span>
    <span class="n">b</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">4</span><span class="p">])</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">gauss</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">)</span>
    <span class="k">print</span><span class="p">(</span><span class="s">&#39;Gauss result is x = </span><span class="se">\n</span><span class="s"> </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> <span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</td></tr>
</table>
<p>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: <span class="math">\(x =[1.3333, 3.1667,1.5000,-6.0000]\)</span>.
Also, the <span class="math">\(A\)</span> and <span class="math">\(b\)</span> (denoted as <span class="math">\(C=[A\quad b]\)</span> below) are updated as:</p>
<blockquote>
<div><div class="math">
<div><div class="math">
\[\begin{split}C^1=
\[\begin{split}C^1=

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.

<thead valign="bottom"> </thead>
MATLAB code Python code
<span class="k">function</span><span class="w"> </span>x <span class="p">=</span><span class="w"> </span><span class="nf">gauss</span><span class="p">(</span>A,b<span class="p">)</span><span class="w"></span>
<span class="c">% This function performs the Gauss elimination without pivoting</span>
<span class="c">% </span>
<span class="c">% x = GAUSS(A, b)</span>

<span class="p">[</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">]</span> <span class="p">=</span> <span class="nb">size</span><span class="p">(</span><span class="n">A</span><span class="p">);</span> 

<span class="c">% Check for zero diagonal elements</span>
<span class="k">if</span> <span class="n">any</span><span class="p">(</span><span class="nb">diag</span><span class="p">(</span><span class="n">A</span><span class="p">)</span><span class="o">==</span><span class="mi">0</span><span class="p">)</span>
    <span class="n">error</span><span class="p">(</span><span class="s">'Division by zero will occur; pivoting not supported'</span><span class="p">)</span>
<span class="k">end</span>

<span class="c">% Forward elimination</span>
<span class="k">for</span> <span class="n">row</span><span class="p">=</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span>
    <span class="k">for</span> <span class="nb">i</span><span class="p">=</span><span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span><span class="n">n</span>
        <span class="nb">factor</span> <span class="p">=</span> <span class="n">A</span><span class="p">(</span><span class="nb">i</span><span class="p">,</span><span class="n">row</span><span class="p">)</span> <span class="o">/</span> <span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">);</span>
        <span class="k">for</span> <span class="nb">j</span> <span class="p">=</span> <span class="n">row</span><span class="p">:</span><span class="n">n</span>
            <span class="n">A</span><span class="p">(</span><span class="nb">i</span><span class="p">,</span><span class="nb">j</span><span class="p">)</span> <span class="p">=</span> <span class="n">A</span><span class="p">(</span><span class="nb">i</span><span class="p">,</span><span class="nb">j</span><span class="p">)</span> <span class="o">-</span> <span class="nb">factor</span><span class="o">*</span><span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="nb">j</span><span class="p">);</span>
        <span class="k">end</span>
        <span class="n">b</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="p">=</span> <span class="n">b</span><span class="p">(</span><span class="nb">i</span><span class="p">)</span> <span class="o">-</span> <span class="nb">factor</span><span class="o">*</span><span class="n">b</span><span class="p">(</span><span class="n">row</span><span class="p">);</span>
    <span class="k">end</span>
  <span class="n">A_and_b</span> <span class="p">=</span>  <span class="p">[</span><span class="n">A</span> <span class="n">b</span><span class="p">]</span>
<span class="k">end</span>

<span class="c">% Backward substitution</span>
<span class="n">x</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="p">=</span> <span class="n">b</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="o">/</span> <span class="n">A</span><span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">);</span>
<span class="k">for</span> <span class="n">row</span> <span class="p">=</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span><span class="mi">1</span>
    <span class="n">sums</span> <span class="p">=</span> <span class="n">b</span><span class="p">(</span><span class="n">row</span><span class="p">);</span>
    <span class="k">for</span> <span class="nb">j</span> <span class="p">=</span> <span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">:</span> <span class="n">n</span>
        <span class="n">sums</span> <span class="p">=</span> <span class="n">sums</span> <span class="o">-</span> <span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="nb">j</span><span class="p">)</span> <span class="o">*</span> <span class="n">x</span><span class="p">(</span><span class="nb">j</span><span class="p">);</span>
    <span class="k">end</span>
    <span class="n">x</span><span class="p">(</span><span class="n">row</span><span class="p">)</span> <span class="p">=</span> <span class="n">sums</span> <span class="o">/</span> <span class="n">A</span><span class="p">(</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">);</span>
<span class="k">end</span>

<span class="c">%Now run this in the command window or an M-file:</span>
<span class="c">% >> A = [1 1 1 1;1 -2 -2 -2;1 4 -4 1;1 -5 -5 -3];</span>
<span class="c">% >> b = [0 4 2 -4]';</span>
<span class="c">% >> x = gauss(A,b)</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="kn">as</span> <span class="nn">np</span>

<span class="k">def</span> <span class="nf">forward_elimination</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
    <span class="sd">"""</span>
<span class="sd">    Calculates the forward part of Gaussian elimination.</span>
<span class="sd">    """</span>
    <span class="k">for</span> <span class="n">row</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
            <span class="n">factor</span> <span class="o">=</span> <span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">row</span><span class="p">]</span> <span class="o">/</span> <span class="n">A</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">]</span>
            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">row</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
                <span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">factor</span> <span class="o">*</span> <span class="n">A</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">j</span><span class="p">]</span>

            <span class="n">b</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span> <span class="n">factor</span> <span class="o">*</span> <span class="n">b</span><span class="p">[</span><span class="n">row</span><span class="p">]</span>

        <span class="k">print</span><span class="p">(</span><span class="s">'A = </span><span class="se">\n</span><span class="si">%s</span><span class="s"> and b = </span><span class="si">%s</span><span class="s">'</span> <span class="o">%</span> <span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="n">b</span><span class="p">))</span>
    <span class="k">return</span> <span class="n">A</span><span class="p">,</span> <span class="n">b</span>

<span class="k">def</span> <span class="nf">back_substitution</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
    <span class="sd">""""</span>
<span class="sd">    Does back substitution, returns the Gauss result.</span>
<span class="sd">    """</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
    <span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">/</span> <span class="n">a</span><span class="p">[</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
    <span class="k">for</span> <span class="n">row</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">):</span>
        <span class="n">sums</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">row</span><span class="p">]</span>
        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">row</span><span class="o">+</span><span class="mi">1</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
            <span class="n">sums</span> <span class="o">=</span> <span class="n">sums</span> <span class="o">-</span> <span class="n">a</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">j</span><span class="p">]</span> <span class="o">*</span> <span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
        <span class="n">x</span><span class="p">[</span><span class="n">row</span><span class="p">]</span> <span class="o">=</span> <span class="n">sums</span> <span class="o">/</span> <span class="n">a</span><span class="p">[</span><span class="n">row</span><span class="p">,</span><span class="n">row</span><span class="p">]</span>
    <span class="k">return</span> <span class="n">x</span>

<span class="k">def</span> <span class="nf">gauss</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
    <span class="sd">"""</span>
<span class="sd">    This function performs Gauss elimination without pivoting.</span>
<span class="sd">    """</span>
    <span class="n">n</span> <span class="o">=</span> <span class="n">A</span><span class="o">.</span><span class="n">shape</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>

    <span class="c"># Check for zero diagonal elements</span>
    <span class="k">if</span> <span class="nb">any</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">A</span><span class="p">)</span><span class="o">==</span><span class="mi">0</span><span class="p">):</span>
        <span class="k">raise</span> <span class="ne">ZeroDivisionError</span><span class="p">((</span><span class="s">'Division by zero will occur; '</span>
                                  <span class="s">'pivoting currently not supported'</span><span class="p">))</span>

    <span class="n">A</span><span class="p">,</span> <span class="n">b</span> <span class="o">=</span> <span class="n">forward_elimination</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">back_substitution</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span>

<span class="c"># Main program starts here</span>
<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">'__main__'</span><span class="p">:</span>
    <span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([[</span><span class="mi">1</span><span class="p">,</span>   <span class="mi">1</span><span class="p">,</span>   <span class="mi">1</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                  <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span>  <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">],</span>
                  <span class="p">[</span><span class="mi">1</span><span class="p">,</span>   <span class="mi">4</span><span class="p">,</span> <span class="o">-</span><span class="mi">4</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span>
                  <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">5</span><span class="p">,</span>  <span class="o">-</span><span class="mi">5</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">]])</span>
    <span class="n">b</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span> <span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">4</span><span class="p">])</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">gauss</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">)</span>
    <span class="k">print</span><span class="p">(</span><span class="s">'Gauss result is x = </span><span class="se">\n</span><span class="s"> </span><span class="si">%s</span><span class="s">'</span> <span class="o">%</span> <span class="n">x</span><span class="p">)</span>

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