# Worksheet 1


```
%matplotlib inline
```


```
import ngcm_utils as ngcm
```

The first worksheet covers basic topics in linear algebra. There is also a basic question on nonlinear root-finding.

## Question 1

Write down the 1, 2 and $\infty$ vector norms of

\begin{equation}
  {\bf v}_1 = \begin{pmatrix} 1 \\ 3 \\ -1 \end{pmatrix}, \quad {\bf v}_2 = \begin{pmatrix} 1 \\ -2 \end{pmatrix}, \quad {\bf v}_3 = \begin{pmatrix} 1 \\ 6 \\ -3 \\ 1 \end{pmatrix}.
\end{equation}

### Answer Question 1

We know that the 1-norm is the sum of the absolute values, so

\begin{align}
  | {\bf v}_{1} |_1 &= 1 + 3 + 1 = 5, \\
  | {\bf v}_{2} |_1 &= 1 + 2 = 3, \\
  | {\bf v}_{2} |_1 &= 1 + 6 + 3 + 1 = 11.
\end{align}

The 2-norm is the square root of the sum of the squares, so

\begin{align}
  | {\bf v}_{1} |_2 &= \sqrt{1^2 + 3^2 + 1^2} = \sqrt{11} \approx 3.3166, \\
  | {\bf v}_{2} |_2 &= \sqrt{1^2 + 3^2} = \sqrt{10} \approx 2.2361, \\
  | {\bf v}_{2} |_2 &= \sqrt{1^2 + 6^2 + 3^2 + 1^2} = \sqrt{47} \approx 6.8557.
\end{align}

The $\infty$-norm is the maximum absolute value, so 

\begin{align}
  | {\bf v}_{1} |_{\infty} &= 3, \\
  | {\bf v}_{2} |_{\infty} &= 2, \\
  | {\bf v}_{2} |_{\infty} &= 6.
\end{align}


Next we calculate the condition numbers numerically using python, as a cross-check. Use numpy for this, via its linear algebra subpackage `numpy.linalg`.


```
import numpy as np
import numpy.linalg as la

v1 = np.array([1.0,  3.0, -1.0])
v2 = np.array([1.0, -2.0])
v3 = np.array([1.0,  6.0, -3.0, 1.0])

for norm in [1, 2, np.inf]:
    for v in [v1, v2, v3]:
        print("The {0} norm of {1} is {2:.3}".\
              format(norm, v, la.norm(v, norm)))
```

    The 1 norm of [ 1.  3. -1.] is 5.0
    The 1 norm of [ 1. -2.] is 3.0
    The 1 norm of [ 1.  6. -3.  1.] is 11.0
    The 2 norm of [ 1.  3. -1.] is 3.32
    The 2 norm of [ 1. -2.] is 2.24
    The 2 norm of [ 1.  6. -3.  1.] is 6.86
    The inf norm of [ 1.  3. -1.] is 3.0
    The inf norm of [ 1. -2.] is 2.0
    The inf norm of [ 1.  6. -3.  1.] is 6.0


To show how this would be done in Matlab, we use the Octave mode here (note that Octave is a free software Matlab clone: not as full featured, but using very similar syntax).


```
%load_ext octavemagic
```


```octave
%%octave

v1 = [1,  3, -1];
v2 = [1, -2];
v3 = [1,  6, -3, 1];

printf("The 1 norm of v1 is %f.\n", norm(v1, 1))
printf("The 1 norm of v2 is %f.\n", norm(v2, 1))
printf("The 1 norm of v3 is %f.\n", norm(v3, 1))
printf("The 2 norm of v1 is %f.\n", norm(v1, 2))
printf("The 2 norm of v2 is %f.\n", norm(v2, 2))
printf("The 2 norm of v3 is %f.\n", norm(v3, 2))
printf("The inf norm of v1 is %f.\n", norm(v1, "inf"))
printf("The inf norm of v2 is %f.\n", norm(v2, "inf"))
printf("The inf norm of v3 is %f.\n", norm(v3, "inf"))
```


    The 1 norm of v1 is 5.000000.
    The 1 norm of v2 is 3.000000.
    The 1 norm of v3 is 11.000000.
    The 2 norm of v1 is 3.316625.
    The 2 norm of v2 is 2.236068.
    The 2 norm of v3 is 6.855655.
    The inf norm of v1 is 3.000000.
    The inf norm of v2 is 2.000000.
    The inf norm of v3 is 6.000000.


## Question 2

Find the 1 and $\infty$ norms of

\begin{equation}
  A_1 = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}, \quad A_2 = \begin{pmatrix} -3 & 2 \\ 3 & 6 \end{pmatrix}.
\end{equation}

### Answer Question 2

The 1-norm of a matrix is the maximum of the 1-norms of the **column** vectors. For $A_1$ the 1-norms of the column vectors are 4 and 6 respectively. For $A_2$ they are 6 and 8 respectively. So we have

\begin{align}
  \| A_1 \|_1 &= 6, \\
  \| A_2 \|_1 & = 8.
\end{align}

The $\infty$ norm of a matrix is the maximum of the 1-norms of the **row** vectors. For $A_1$ the 1-norms of the row vectors are 3 and 7 respectively. For $A_2$ they are 5 and 9 respectively. So we have

\begin{align}
  \| A_1 \|_{\infty} &= 7, \\
  \| A_2 \|_{\infty} & = 9.
\end{align}


It is equally straightforward to repeat this calculation with python, using the same approach.


```
A1 = np.array([[ 1.0, 2.0],[3.0, 4.0]])
A2 = np.array([[-3.0, 2.0],[3.0, 6.0]])

for norm in [1, np.inf]:
    for A in [A1, A2]:
        print("The {0} norm of \n{1}\nis {2:.3}\n".\
              format(norm, A, la.norm(A, norm)))
```

    The 1 norm of 
    [[ 1.  2.]
     [ 3.  4.]]
    is 6.0
    
    The 1 norm of 
    [[-3.  2.]
     [ 3.  6.]]
    is 8.0
    
    The inf norm of 
    [[ 1.  2.]
     [ 3.  4.]]
    is 7.0
    
    The inf norm of 
    [[-3.  2.]
     [ 3.  6.]]
    is 9.0
    


Again, using Matlab syntax we would do the following.


```octave
%%octave

A1 = [ 1 2; 3 4];
A2 = [-3 2; 3 6];

printf("The 1 norm of A1 is %f.\n", norm(A1, 1))
printf("The 1 norm of A2 is %f.\n", norm(A2, 1))
printf("The inf norm of A1 is %f.\n", norm(A1, "inf"))
printf("The inf norm of A2 is %f.\n", norm(A2, "inf"))
```


    The 1 norm of A1 is 6.000000.
    The 1 norm of A2 is 8.000000.
    The inf norm of A1 is 7.000000.
    The inf norm of A2 is 9.000000.


## Question 3

Find the condition number of the above matrices.

### Answer Question 3

The analytic calculation of the condition number requires the norm of the inverse matrices, which is **not** the inverse of the norm of the matrix. This is, of course, not useful for the numerical work.

The inverse matrices are

\begin{align}
  A_1^{-1} & = - \frac{1}{2} \begin{pmatrix} 4 & -2 \\ -3 & 1 \end{pmatrix}, \\
  A_2^{-1} & = - \frac{1}{24} \begin{pmatrix} 6 & -2 \\ -3 & -3 \end{pmatrix}.
\end{align}

It follows that the matrix norms of the inverse matrices are

\begin{align}
  \| A_1^{-1} \| & = \frac{7}{2}, & \| A_2^{-1} \|_1 & = \frac{3}{8}, \\
  \| A_1^{-1} \|_{\infty} & = 3, & \| A_2^{-1} \|_{\infty} & = \frac{1}{3}.
\end{align}

Therefore the condition numbers with respect to the 1-norm are

\begin{align}
  K(A_1) & = \| A_1 \|_1 \| A_1^{-1} \|_1 \\ & = 21, \\
  K(A_2) & = \| A_2 \|_1 \| A_2^{-1} \|_1 \\ & = 3,
\end{align}

and the condition numbers with respect to the $\infty$-norm are

\begin{align}
  K(A_1) & = \| A_1 \|_{\infty} \| A_1^{-1} \|_{\infty} \\ & = 21, \\
  K(A_2) & = \| A_2 \|_{\infty} \| A_2^{-1} \|_{\infty} \\ & = 3.
\end{align}

Note that it is chance that, in this case, they are identical: this will not be true in general.

This suggests that, if the solution of a linear system is needed as part of a numerical method, any intrinsic errors will be increased more by the matrix $A_1$ than by the matrix $A_2$, by a factor $\sim 10$.

When calculating the condition number numerically, the Singular Value Decomposition (SVD) is typically used, which does not require inverting the matrix. Thus calculating the condition number of fast, and should be done before any matrix operation to check the condition of the matrix. Within python, this is straightforward.


```
A1 = np.array([[ 1.0, 2.0],[3.0, 4.0]])
A2 = np.array([[-3.0, 2.0],[3.0, 6.0]])

for norm in [1, 2, np.inf]:
    for A in [A1, A2]:
        print("The condition number with respect to the"\
              " {0} norm of\n{1}\nis {2:.3}\n"\
              .format(norm, A, la.cond(A, norm)))
```

    The condition number with respect to the 1 norm of
    [[ 1.  2.]
     [ 3.  4.]]
    is 21.0
    
    The condition number with respect to the 1 norm of
    [[-3.  2.]
     [ 3.  6.]]
    is 3.0
    
    The condition number with respect to the 2 norm of
    [[ 1.  2.]
     [ 3.  4.]]
    is 14.9
    
    The condition number with respect to the 2 norm of
    [[-3.  2.]
     [ 3.  6.]]
    is 1.89
    
    The condition number with respect to the inf norm of
    [[ 1.  2.]
     [ 3.  4.]]
    is 21.0
    
    The condition number with respect to the inf norm of
    [[-3.  2.]
     [ 3.  6.]]
    is 3.0
    


Again, similar functions exist in Matlab.


```octave
%%octave

A1 = [ 1 2; 3 4];
A2 = [-3 2; 3 6];

printf("The condition number with respect to the 1 norm of A1 is %f.\n", \
       cond(A1, 1))
printf("The condition number with respect to the 1 norm of A2 is %f.\n", \
       cond(A2, 1))
printf("The condition number with respect to the 2 norm of A1 is %f.\n", \
       cond(A1, 2))
printf("The condition number with respect to the 2 norm of A2 is %f.\n", \
       cond(A2, 2))
printf("The condition number with respect to the inf norm of A1 is %f.\n", \
       cond(A1, "inf"))
printf("The condition number with respect to the inf norm of A2 is %f.\n", \
       cond(A2, "inf"))

```


    The condition number with respect to the 1 norm of A1 is 21.000000.
    The condition number with respect to the 1 norm of A2 is 3.000000.
    The condition number with respect to the 2 norm of A1 is 14.933034.
    The condition number with respect to the 2 norm of A2 is 1.886618.
    The condition number with respect to the inf norm of A1 is 21.000000.
    The condition number with respect to the inf norm of A2 is 3.000000.


## Question 4

Explain the difference between direct and indirect methods for solving linear systems. Give an example of when the latter may be more useful.

*This is a standard part of an exam question: see, e.g., 07/08*.

### Answer Question 4

*Direct methods* consist of a finite list of transformations of the original matrix of the coefficients that reduce the linear systems to one that is easily solved. 

*Indirect* or *iterative methods*, consist of algorithms that specify a series of steps, possibly infinite, that lead closer and closer to the solution; there may not be a guarantee that they ever exactly reach it. This may not seem a very desirable feature until we remember that we cannot in any case perfectly represent an exact solution: most iterative methods provide us with a highly accurate solution in relatively few iterations.

Large, sparse matrices are ideally solved using iterative methods.

## Coding Question 1

For each of the matrices above, work out their transpose and inverse using standard python commands.

### Answer Coding Question 1


```
A1 = np.array([[ 1.0, 2.0],[3.0, 4.0]])
A2 = np.array([[-3.0, 2.0],[3.0, 6.0]])

for A in [A1, A2]:
    print("The transpose of \n{0}\nis \n{1}\n".format(A, np.transpose(A)))
    print("The inverse of \n{0}\nis \n{1}\n".format(A, la.inv(A)))
```

    The transpose of 
    [[ 1.  2.]
     [ 3.  4.]]
    is 
    [[ 1.  3.]
     [ 2.  4.]]
    
    The inverse of 
    [[ 1.  2.]
     [ 3.  4.]]
    is 
    [[-2.   1. ]
     [ 1.5 -0.5]]
    
    The transpose of 
    [[-3.  2.]
     [ 3.  6.]]
    is 
    [[-3.  3.]
     [ 2.  6.]]
    
    The inverse of 
    [[-3.  2.]
     [ 3.  6.]]
    is 
    [[-0.25        0.08333333]
     [ 0.125       0.125     ]]
    


Once again, we can use Matlab functions to do this.


```octave
%%octave

A1 = [ 1 2; 3 4];
A2 = [-3 2; 3 6];

printf("The transpose of A1 is \n")
disp(A1.')
printf("The inverse of A1 is \n")
disp(inv(A1))
printf("The transpose of A2 is \n")
disp(A2.')
printf("The inverse of A2 is \n")
disp(inv(A2))
```


    The transpose of A1 is
            1        3
            2        4
    The inverse of A1 is
      -2.00000  1.00000
      1.50000  -0.50000
    The transpose of A2 is
           -3        3
            2        6
    The inverse of A2 is
      -0.25000  0.08333
      0.12500  0.12500


## Coding Question 3

Note that coding question 2 has been effectively answered above, in the "theory" section.

Write a function that takes a matrix, computes its condition number, and reports whether the matrix is suitably well-conditions. The choice of criteria is up to you.

### Answer Coding Question 3

As we are now starting to use functions, we could define them in the notebook, or in a separate file. Here I have used a separate file : `Worksheet1_Functions.py`. The function is defined as follows.


```
ngcm.highlight_source("Worksheet1_Functions", "MatrixConditionCheck")
```




<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=None">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #408080; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #008000; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
body .cp { color: #BC7A00 } /* Comment.Preproc */
body .c1 { color: #408080; font-style: italic } /* Comment.Single */
body .cs { color: #408080; font-style: italic } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #008000 } /* Keyword.Pseudo */
body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #B00040 } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BA2121 } /* Literal.String */
body .na { color: #7D9029 } /* Name.Attribute */
body .nb { color: #008000 } /* Name.Builtin */
body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #0000FF } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #19177C } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #BA2121 } /* Literal.String.Backtick */
body .sc { color: #BA2121 } /* Literal.String.Char */
body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BA2121 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BA2121 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BA2121 } /* Literal.String.Single */
body .ss { color: #19177C } /* Literal.String.Symbol */
body .bp { color: #008000 } /* Name.Builtin.Pseudo */
body .vc { color: #19177C } /* Name.Variable.Class */
body .vg { color: #19177C } /* Name.Variable.Global */
body .vi { color: #19177C } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span class="k">def</span> <span class="nf">MatrixConditionCheck</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">MaxConditionNumber</span> <span class="o">=</span> <span class="mf">10.0</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Check the condition number of a matrix.</span>
<span class="sd">    Only write output to screen if the condition number is too high.</span>
<span class="sd">    Should return something, really.&quot;&quot;&quot;</span>

    <span class="n">ConditionNumber</span> <span class="o">=</span> <span class="n">la</span><span class="o">.</span><span class="n">cond</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">ConditionNumber</span> <span class="o">&gt;</span> <span class="n">MaxConditionNumber</span><span class="p">:</span>
        <span class="k">print</span><span class="p">(</span><span class="s">&quot;The condition number of the matrix</span><span class="se">\n</span><span class="s">{0}</span><span class="se">\n</span><span class="s">&quot;</span>\
              <span class="s">&quot;is too large (i.e., it is {1:.4} which is larger&quot;</span>\
              <span class="s">&quot; than {2:.4}).</span><span class="se">\n</span><span class="s">&quot;</span><span class="o">.</span>\
              <span class="n">format</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">ConditionNumber</span><span class="p">,</span> <span class="n">MaxConditionNumber</span><span class="p">))</span>
</pre></div>
</body>
</html>




We can then import the module and check the function using the earlier input matrices.


```
import Worksheet1_Functions as w1


A1 = np.array([[ 1.0, 2.0],[3.0, 4.0]])
A2 = np.array([[-3.0, 2.0],[3.0, 6.0]])

w1.MatrixConditionCheck(A1)
w1.MatrixConditionCheck(A2)
```

    The condition number of the matrix
    [[ 1.  2.]
     [ 3.  4.]]
    is too large (i.e., it is 14.93 which is larger than 10.0).
    


Note the Matlab implementation in `Worksheet1_MyConditionNumber.m`.


```
ngcm.highlight_source_matlab("Worksheet1_MyConditionNumber.m")
```




<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=None">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #408080; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #008000; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
body .cp { color: #BC7A00 } /* Comment.Preproc */
body .c1 { color: #408080; font-style: italic } /* Comment.Single */
body .cs { color: #408080; font-style: italic } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #008000 } /* Keyword.Pseudo */
body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #B00040 } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BA2121 } /* Literal.String */
body .na { color: #7D9029 } /* Name.Attribute */
body .nb { color: #008000 } /* Name.Builtin */
body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #0000FF } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #19177C } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #BA2121 } /* Literal.String.Backtick */
body .sc { color: #BA2121 } /* Literal.String.Char */
body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BA2121 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BA2121 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BA2121 } /* Literal.String.Single */
body .ss { color: #19177C } /* Literal.String.Symbol */
body .bp { color: #008000 } /* Name.Builtin.Pseudo */
body .vc { color: #19177C } /* Name.Variable.Class */
body .vg { color: #19177C } /* Name.Variable.Global */
body .vi { color: #19177C } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span class="c">%</span>
<span class="c">% function ConditionNumber = Worksheet1_MyConditionNumber(A)</span>
<span class="c">%</span>
<span class="c">% Computes the condition number of the input matrix A.</span>
<span class="c">%</span>
<span class="c">% Expects a square 2-D matrix A. It will complain if the condition number</span>
<span class="c">% is larger than 10 (this is a ludicrously low bound, but tests that the</span>
<span class="c">% function actually works with the default input).</span>
<span class="c">%</span>
<span class="k">function</span><span class="w"> </span>ConditionNumber <span class="p">=</span><span class="w"> </span><span class="nf">Worksheet1_MyConditionNumber</span><span class="p">(</span>A<span class="p">)</span><span class="w"></span>

<span class="c">% Check the input is reasonable</span>
<span class="k">if</span> <span class="p">(</span><span class="n">not</span><span class="p">(</span><span class="n">isnumeric</span><span class="p">(</span><span class="n">A</span><span class="p">)))</span>
    <span class="n">error</span><span class="p">(</span><span class="s">&#39;The input must be a numerical array!&#39;</span><span class="p">);</span>
<span class="k">elseif</span> <span class="p">(</span><span class="nb">ndims</span><span class="p">(</span><span class="n">A</span><span class="p">)</span><span class="o">~=</span><span class="mi">2</span><span class="p">)</span>
    <span class="n">error</span><span class="p">(</span><span class="s">&#39;The input must be a 2-d array!&#39;</span><span class="p">);</span>
<span class="k">elseif</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="mi">1</span><span class="p">)</span><span class="o">~=</span><span class="nb">size</span><span class="p">(</span><span class="n">A</span><span class="p">,</span><span class="mi">2</span><span class="p">))</span>
    <span class="n">error</span><span class="p">(</span><span class="s">&#39;The input must be a square 2-d array!&#39;</span><span class="p">);</span>
<span class="k">end</span>

<span class="c">% Set up the maximum condition number</span>
<span class="n">Max_ConditionNumber</span> <span class="p">=</span> <span class="mi">10</span><span class="p">;</span>

<span class="c">% Compute the condition number</span>
<span class="n">ConditionNumber</span> <span class="p">=</span> <span class="n">cond</span><span class="p">(</span><span class="n">A</span><span class="p">);</span>

<span class="k">if</span> <span class="p">(</span><span class="n">ConditionNumber</span> <span class="o">&gt;</span> <span class="n">Max_ConditionNumber</span><span class="p">)</span>
    <span class="nb">disp</span><span class="p">(</span><span class="n">sprintf</span><span class="p">(</span><span class="s">&#39;The condition number, %g, of the matrix makes it not suitably well-conditioned (bigger than %g)&#39;</span><span class="p">,</span> <span class="c">...</span>
        <span class="n">ConditionNumber</span><span class="p">,</span> <span class="n">Max_ConditionNumber</span><span class="p">));</span>
<span class="k">end</span>

<span class="k">end</span>
</pre></div>
</body>
</html>





```octave
%%octave

A1 = [ 1 2; 3 4];
A2 = [-3 2; 3 6];

Worksheet1_MyConditionNumber(A1);
Worksheet1_MyConditionNumber(A2);
```


    The condition number, 14.933, of the matrix makes it not suitably well-conditioned (bigger than 10)


## Coding Question 4

Write a bisection method to find the root of

\begin{equation}
  f(x) = \tan (x) - e^{-x}, \quad x \in [0, 1].
\end{equation}

### Answer Coding Question 4

Again we have defined the function in the separate file.


```
ngcm.highlight_source("Worksheet1_Functions", "bisection")
```




<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=None">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #408080; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #008000; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
body .cp { color: #BC7A00 } /* Comment.Preproc */
body .c1 { color: #408080; font-style: italic } /* Comment.Single */
body .cs { color: #408080; font-style: italic } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #008000 } /* Keyword.Pseudo */
body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #B00040 } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BA2121 } /* Literal.String */
body .na { color: #7D9029 } /* Name.Attribute */
body .nb { color: #008000 } /* Name.Builtin */
body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #0000FF } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #19177C } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #BA2121 } /* Literal.String.Backtick */
body .sc { color: #BA2121 } /* Literal.String.Char */
body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BA2121 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BA2121 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BA2121 } /* Literal.String.Single */
body .ss { color: #19177C } /* Literal.String.Symbol */
body .bp { color: #008000 } /* Name.Builtin.Pseudo */
body .vc { color: #19177C } /* Name.Variable.Class */
body .vg { color: #19177C } /* Name.Variable.Global */
body .vi { color: #19177C } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span class="k">def</span> <span class="nf">bisection</span><span class="p">(</span><span class="n">f</span><span class="p">,</span> <span class="n">interval</span><span class="p">,</span> <span class="n">tolerance</span> <span class="o">=</span> <span class="mf">1.e-10</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;General bisection method for a function f of one variable.     </span>
<span class="sd">    There must be at least one root within the interval.     </span>
<span class="sd">    Default tolerance (width of the interval) is 1e-10.&quot;&quot;&quot;</span>

    <span class="k">assert</span> <span class="nb">len</span><span class="p">(</span><span class="n">interval</span><span class="p">)</span> <span class="o">==</span> <span class="mi">2</span>

    <span class="c"># Get the endpoints of the interval</span>
    <span class="p">[</span><span class="n">x_min</span><span class="p">,</span> <span class="n">x_max</span><span class="p">]</span> <span class="o">=</span> <span class="n">interval</span>

    <span class="c"># Values at the ends of the domain</span>
    <span class="n">f_min</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x_min</span><span class="p">)</span>
    <span class="n">f_max</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x_max</span><span class="p">)</span>

    <span class="c"># Check that at least one root lies within the interval</span>
    <span class="k">assert</span><span class="p">(</span><span class="n">f_min</span> <span class="o">*</span> <span class="n">f_max</span> <span class="o">&lt;</span> <span class="mf">0.0</span><span class="p">)</span>

    <span class="c"># The loop</span>
    <span class="n">x_c</span> <span class="o">=</span> <span class="p">(</span><span class="n">x_min</span> <span class="o">+</span> <span class="n">x_max</span><span class="p">)</span> <span class="o">/</span> <span class="mf">2.0</span>
    <span class="n">f_c</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x_c</span><span class="p">)</span>
    <span class="n">iteration</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">while</span> <span class="p">((</span><span class="n">x_max</span> <span class="o">-</span> <span class="n">x_min</span> <span class="o">&gt;</span> <span class="n">tolerance</span><span class="p">)</span> <span class="ow">and</span> \
               <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">f_c</span><span class="p">)</span> <span class="o">&gt;</span> <span class="n">tolerance</span><span class="p">)</span> <span class="ow">and</span> \
               <span class="p">(</span><span class="n">iteration</span> <span class="o">&lt;</span> <span class="mi">100</span><span class="p">)):</span>
        <span class="n">iteration</span> <span class="o">=</span> <span class="n">iteration</span><span class="o">+</span><span class="mi">1</span>    
        <span class="k">if</span> <span class="n">f_min</span> <span class="o">*</span> <span class="n">f_c</span> <span class="o">&lt;</span> <span class="mf">0.0</span><span class="p">:</span>
            <span class="n">x_max</span> <span class="o">=</span> <span class="n">x_c</span>
            <span class="n">f_max</span> <span class="o">=</span> <span class="n">f_c</span>
        <span class="k">else</span><span class="p">:</span>
            <span class="n">x_min</span> <span class="o">=</span> <span class="n">x_c</span>
            <span class="n">f_min</span> <span class="o">=</span> <span class="n">f_c</span>
        <span class="n">x_c</span> <span class="o">=</span> <span class="p">(</span><span class="n">x_min</span> <span class="o">+</span> <span class="n">x_max</span><span class="p">)</span> <span class="o">/</span> <span class="mf">2.0</span>
        <span class="n">f_c</span> <span class="o">=</span> <span class="n">f</span><span class="p">(</span><span class="n">x_c</span><span class="p">)</span>

    <span class="k">print</span><span class="p">(</span><span class="s">&quot;The root is approximately {0} where &quot;</span>\
          <span class="s">&quot;f is {1:.4} (tolerance {2:.4})&quot;</span><span class="o">.</span><span class="n">format</span><span class="p">(</span><span class="n">x_c</span><span class="p">,</span> <span class="n">f_c</span><span class="p">,</span> <span class="n">tolerance</span><span class="p">))</span>
    <span class="k">return</span> <span class="n">x_c</span>
</pre></div>
</body>
</html>




We have already imported the module. Also within that module is the function whose root we are trying to find.


```
ngcm.highlight_source("Worksheet1_Functions", "fn_worksheet1_q4")
```




<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01//EN"
   "http://www.w3.org/TR/html4/strict.dtd">

<html>
<head>
  <title></title>
  <meta http-equiv="content-type" content="text/html; charset=None">
  <style type="text/css">
td.linenos { background-color: #f0f0f0; padding-right: 10px; }
span.lineno { background-color: #f0f0f0; padding: 0 5px 0 5px; }
pre { line-height: 125%; }
body .hll { background-color: #ffffcc }
body  { background: #f8f8f8; }
body .c { color: #408080; font-style: italic } /* Comment */
body .err { border: 1px solid #FF0000 } /* Error */
body .k { color: #008000; font-weight: bold } /* Keyword */
body .o { color: #666666 } /* Operator */
body .cm { color: #408080; font-style: italic } /* Comment.Multiline */
body .cp { color: #BC7A00 } /* Comment.Preproc */
body .c1 { color: #408080; font-style: italic } /* Comment.Single */
body .cs { color: #408080; font-style: italic } /* Comment.Special */
body .gd { color: #A00000 } /* Generic.Deleted */
body .ge { font-style: italic } /* Generic.Emph */
body .gr { color: #FF0000 } /* Generic.Error */
body .gh { color: #000080; font-weight: bold } /* Generic.Heading */
body .gi { color: #00A000 } /* Generic.Inserted */
body .go { color: #888888 } /* Generic.Output */
body .gp { color: #000080; font-weight: bold } /* Generic.Prompt */
body .gs { font-weight: bold } /* Generic.Strong */
body .gu { color: #800080; font-weight: bold } /* Generic.Subheading */
body .gt { color: #0044DD } /* Generic.Traceback */
body .kc { color: #008000; font-weight: bold } /* Keyword.Constant */
body .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */
body .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */
body .kp { color: #008000 } /* Keyword.Pseudo */
body .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */
body .kt { color: #B00040 } /* Keyword.Type */
body .m { color: #666666 } /* Literal.Number */
body .s { color: #BA2121 } /* Literal.String */
body .na { color: #7D9029 } /* Name.Attribute */
body .nb { color: #008000 } /* Name.Builtin */
body .nc { color: #0000FF; font-weight: bold } /* Name.Class */
body .no { color: #880000 } /* Name.Constant */
body .nd { color: #AA22FF } /* Name.Decorator */
body .ni { color: #999999; font-weight: bold } /* Name.Entity */
body .ne { color: #D2413A; font-weight: bold } /* Name.Exception */
body .nf { color: #0000FF } /* Name.Function */
body .nl { color: #A0A000 } /* Name.Label */
body .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */
body .nt { color: #008000; font-weight: bold } /* Name.Tag */
body .nv { color: #19177C } /* Name.Variable */
body .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */
body .w { color: #bbbbbb } /* Text.Whitespace */
body .mf { color: #666666 } /* Literal.Number.Float */
body .mh { color: #666666 } /* Literal.Number.Hex */
body .mi { color: #666666 } /* Literal.Number.Integer */
body .mo { color: #666666 } /* Literal.Number.Oct */
body .sb { color: #BA2121 } /* Literal.String.Backtick */
body .sc { color: #BA2121 } /* Literal.String.Char */
body .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */
body .s2 { color: #BA2121 } /* Literal.String.Double */
body .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */
body .sh { color: #BA2121 } /* Literal.String.Heredoc */
body .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */
body .sx { color: #008000 } /* Literal.String.Other */
body .sr { color: #BB6688 } /* Literal.String.Regex */
body .s1 { color: #BA2121 } /* Literal.String.Single */
body .ss { color: #19177C } /* Literal.String.Symbol */
body .bp { color: #008000 } /* Name.Builtin.Pseudo */
body .vc { color: #19177C } /* Name.Variable.Class */
body .vg { color: #19177C } /* Name.Variable.Global */
body .vi { color: #19177C } /* Name.Variable.Instance */
body .il { color: #666666 } /* Literal.Number.Integer.Long */

  </style>
</head>
<body>
<h2></h2>

<div class="highlight"><pre><span class="k">def</span> <span class="nf">fn_worksheet1_q4</span><span class="p">(</span><span class="n">x</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Simple function defined in question, f(x) = tan(x) - exp(-x).&quot;&quot;&quot;</span>

    <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">tan</span><span class="p">(</span><span class="n">x</span><span class="p">)</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</body>
</html>




Test the bisection method using various tolerances.


```
# Now find the root.
x_root = w1.bisection(w1.fn_worksheet1_q4, [0, 1])
# Try changing the precision to see what difference it makes
x_root_e6 = w1.bisection(w1.fn_worksheet1_q4, [0, 1], 1e-6)
x_root_e6 = w1.bisection(w1.fn_worksheet1_q4, [0, 1], 1e-15)
```

    The root is approximately 0.531390856602 where f is -9.652e-11 (tolerance 1e-10)
    The root is approximately 0.531391143799 where f is 5.551e-07 (tolerance 1e-06)
    The root is approximately 0.531390856652 where f is 1.11e-16 (tolerance 1e-15)


A quick and dirty Matlab implementation of the bisection algorithm follows. Note that it is not implemented as a function, but directly as a script: it would be straightforward, in a similar fashion to the python above, to implement it as a function.


```octave
%%octave

% Set up the tolerance in the width of the interval
tolerance = 1e-15;

% Set up the initial bracketing interval
x_lo = 0;
x_hi = 1;

% Set up the function and the initial function values.
f = @(x)(tan(x) - exp(-x));
f_lo = f(x_lo);
f_hi = f(x_hi);

% Set up the bisection loop.
x_guess = (x_lo + x_hi) / 2;
f_guess = f(x_guess);
% Loop until the root is bracketed within the tolerance.
while (x_hi - x_lo > tolerance)
   
    if (f_lo * f_guess < 0)
       % The root is between the guessed value and the lower bound
       x_hi = x_guess;
       f_hi = f_guess;
    else
       % The root is between the guessed value and the upper bound
       x_lo = x_guess;
       f_lo = f_guess;
    end
    
    % Set the new bisection estimate
    x_guess = (x_lo + x_hi) / 2;
    f_guess = f(x_guess);
    
end

printf("The root is approximately %f where f is %g.\n", x_guess, f_guess)
```


    The root is approximately 0.531391 where f is -7.77156e-16.



```
from IPython.core.display import HTML
def css_styling():
    styles = open("../../IPythonNotebookStyles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```




<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
    }
    div.cell{
        width:80%;
        margin-left:10% !important;
        margin-right:auto;
    }
    h1 {
        font-family: Verdana, Arial, Helvetica, sans-serif;
	text-align: center;
    }
    h2 {
        font-family: Verdana, Arial, Helvetica, sans-serif;
    }
    h3 {
        font-family: Verdana, Arial, Helvetica, sans-serif;
    }
    div.text_cell_render{
        font-family: Gill, Verdana, Arial, Helvetica, sans-serif;
        line-height: 110%;
        font-size: 120%;
        width:75%;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro", source-code-pro,Consolas, monospace;
    }
/*    .prompt{
        display: None;
    }*/
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 12pt;
        color: #4057A1;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }  
</style>




> (The cell above executes the style for this notebook. It closely follows the style used in the [12 Steps to Navier Stokes](http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/) course.)

1de63cc705c1c7f9cd09cea2d65b8383ec43a74f
ipynb
Jupyter Notebook
Worksheets/Worksheet1_Notebook.ipynb
indranilsinharoy/NumericalMethods
989e0205565131057c9807ed9d55b6c1a5a38d42
MIT
2021-12-01T09:15:04.000Z
2021-12-01T09:15:04.000Z
Worksheets/Worksheet1_Notebook.ipynb
indranilsinharoy/NumericalMethods
989e0205565131057c9807ed9d55b6c1a5a38d42
MIT
Worksheets/Worksheet1_Notebook.ipynb
indranilsinharoy/NumericalMethods
989e0205565131057c9807ed9d55b6c1a5a38d42
MIT
Qwen/Qwen-72B
1. YES
2. YES