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.

<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=[Ab]$ below) are updated as:

$$\begin{array}{r}{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{array}$$

$$\begin{array}{r}{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{array}$$

$$\begin{array}{r}{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{array}$$

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