Difference between revisions of "User:Kevindunn"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
m
m (Replaced content with "<html> <div class="math"> \[\begin{split}C^1= \left[ \begin{array}{ccccc} 1 & 1 & 1 & 1 & 0\\ 0 & -3 & -3 & -3 &...")
 
(2 intermediate revisions by the same user not shown)
Line 1: Line 1:
<div class="section" id="bonus-question-1">
<html>
<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>


<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>
<div class="math">
 
<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">&#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">
\[\begin{split}C^1=
\[\begin{split}C^1=
      \left[ \begin{array}{ccccc}
1  &amp;  1  &amp;  1  &amp;  1  &amp;  0\\
        0  &amp;  -3  &amp;  -3  &amp;  -3  &amp;  4\\
        0  &amp;  3  &amp;  -5  &amp;  0  &amp;  2\\
0  &amp;  -6 &amp;  -6  &amp;  -4  &amp;  -4
      \end{array} \right]\end{split}\]</div>
<div class="math">
\[\begin{split}C^2=
\left[ \begin{array}{ccccc}
\left[ \begin{array}{ccccc}
    1  &amp;  1  &amp;  1  &amp;  1 &amp;  0\\
  1  &amp;  1  &amp;  1  &amp;  1 &amp;  0\\
    0  &amp;  -3  &amp;  -3  &amp;  -3 &amp; 4\\
  0  &amp;  -3  &amp;  -3  &amp;  -3 &amp;   4\\
    0  &amp;  0 &amp;  -8  &amp;  -3 &amp;  6\\
  0  &amp;  3 &amp;  -5 &amp;  0  &amp;  2\\
    0  &amp;  0 &amp;  0  &amp;  2 &amp;  -12
   0  &amp;  -6 &amp;   -6 &amp;  -4  &amp;  -4        \end{array} \right]\end{split}\]</div>
\end{array} \right]\end{split}\]</div>
<script type='text/javascript' src='https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML-full,local/mwMathJaxConfig'> </script>
<div class="math">
</html>
\[\begin{split}C^3=
      \left[ \begin{array}{ccccc}
1  &amp;  1  &amp;  1  &amp;   1  &amp;  0\\
0  &amp;  -3  &amp; -3 &amp;  -3  &amp;  4\\
        0 &amp;  0  &amp;  -8  &amp;  -3  &amp;  6\\
        0  &amp;  0  &amp;  0  &amp;  2  &amp; -12
      \end{array} \right]\end{split}\]</div>
</div></blockquote>
<p>The final result for <span class="math">\(x = [1.333, 3.167, 1.5, -6]\)</span>.</p>
</div>
</div>

Latest revision as of 07:46, 10 February 2017

\[\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}\]

Retrieved from "https://learnche.org/wiki_3E4/index.php?title=User:Kevindunn&oldid=1072"