MCS 471 Lecture Eighteen

Review for Exam I (with partial answers)

The questions below are meant to stimulate the preparation of the exam. They cover some of the most important but not all topics. You must review all course materials: textbook, homework exercises, lecture notes (also consult the course web site), and quizzes.

You may use a formula sheet: postscript version , and version in pdf format .

Consider the representation of floating-point numbers with base 2 and 3 digits in the fraction part. The values for the exponents are between -8 and +8.
1. What is the machine precision of this number system?
2. Represent the numbers 8 and 3 as floating point numbers and illustrate the calculation of 8+3 and 8-3.
3. Give an example where (a+b) + c /= a + (b+c).
Answers:
1. The machine precision is B^(-p) with B = 2 and p = 3, thus 2^(-3).
2. As representations of 8 and 3 in the number system we have
```
          8 = +.100 2^4 
          3 = +.110 2^2 = +.011 2^3 = +.001 2^4
```
  We denormalize the smallest number to perform the operations:
```
          8 = +.100 2^4              8 = +.100 2^4  
        + 3 = +.001 2^4            - 3 = +.001 2^4
        ----------------           ----------------
         10 = +.101 2^4              5 = +.101 2^4
```
3. The machine precision is typically the largest number you can add to one without making a difference with one in the sum. Thus
```
      (2^(-3) + 1) + (-1) /=  2^(-3) + (1 + (-1)),
```
  since the left hand side equals 0 and the right hand side is 2^(-3).

Consider the equation f(x) = e^x - 4 x^3. Find an approximation for the root of f

using four steps of the bisection method on the interval [0,1];
using four steps of the secant method starting at x_0 = 0 and x_1 = 1;
using four steps of Newton's method starting at x = 1.

Answers:

The output of the bisection method (MATLAB):

  a = 0.5000, b = 1.0000, m = 0.5000, f(m) =  1.1487
  a = 0.7500, b = 1.0000, m = 0.7500, f(m) =  0.4295
  a = 0.7500, b = 0.8750, m = 0.8750, f(m) = -0.2808
  a = 0.8125, b = 0.8750, m = 0.8125, f(m) =  0.1080

The output of the secant method (MATLAB):

  x = 4.3827e-01, dx =  5.6173e-01, f(x) =  1.2133e+00,
  x = 7.1143e-01, dx = -2.7316e-01, f(x) =  5.9659e-01,
  x = 9.7568e-01, dx = -2.6425e-01, f(x) = -1.0623e+00,
  x = 8.0647e-01, dx =  1.6922e-01, f(x) =  1.4191e-01,

The output of Newton's method (MATLAB):

  x = 8.619094e-01, dx = 1.380906e-01, f(x) = -1.935306e-01
  x = 8.323491e-01, dx = 2.956031e-02, f(x) = -7.910105e-03
  x = 8.310340e-01, dx = 1.315075e-03, f(x) = -1.527784e-05
  x = 8.310315e-01, dx = 2.549826e-06, f(x) = -5.737322e-11

Apply four steps of the golden section search method to find the minimum of f(x) = e^x - 4 x^3 in the interval [4,8].

Answers:

The output of Golden Section Search (MATLAB):

------------------------------------------------------
    x1        x2        f(x1)      f(x2)      b - a
------------------------------------------------------
5.5279e+00 6.4721e+00 -4.2406e+02 -4.3757e+02 4.0000e+00
6.4721e+00 7.0557e+00 -4.3757e+02 -2.4555e+02 2.4721e+00
6.1115e+00 6.4721e+00 -4.6205e+02 -4.3757e+02 1.5279e+00
5.8885e+00 6.1115e+00 -4.5586e+02 -4.6205e+02 9.4427e-01

Use Newton's method to find the cube root of 27.13.

Answers:

We apply Newton's method on the equation x^3 - 27.13 = 0. A natural start point is 3. Four steps suffice (MATLAB output):

         x              |dx|     |f(x)| 
3.004814814814815e+00  4.81e-03  2.09e-04
3.004807107958115e+00  7.71e-06  5.35e-10
3.004807107938348e+00  1.98e-11  3.55e-15
3.004807107938348e+00  1.31e-16  3.55e-15

The following is the output of a Newton's method running on f(x) = (x-1)^m:

                       x              |dx|     |f(x)| 
             9.250000000000528e-01  2.50e-02  3.16e-05
             9.437499999999878e-01  1.87e-02  1.00e-05
             9.578124999999771e-01  1.41e-02  3.17e-06
             9.683593750005512e-01  1.05e-02  1.00e-06
             9.762695312457415e-01  7.91e-03  3.17e-07
             9.822021484398730e-01  5.93e-03  1.00e-07
             9.866516113385196e-01  4.45e-03  3.17e-08
             9.899887084944722e-01  3.34e-03  1.00e-08
             9.924915313578651e-01  2.50e-03  3.18e-09
             9.943686485208891e-01  1.88e-03  1.01e-09
             9.957764867505785e-01  1.41e-03  3.18e-10
             9.968323666970245e-01  1.06e-03  1.01e-10
             9.976242738902757e-01  7.92e-04  3.19e-11
             9.982182001705794e-01  5.94e-04  1.01e-11
             9.986636385013267e-01  4.45e-04  3.19e-12
             9.989976859298977e-01  3.34e-04  1.01e-12
             9.992483785941424e-01  2.51e-04  3.19e-13
             9.994362408614242e-01  1.88e-04  1.01e-13
             9.995775156995391e-01  1.41e-04  3.21e-14
             9.996838850522098e-01  1.06e-04  9.77e-15

What is m?

Answers:

We take two consecutive values for |dx| to estimate m, using the formula

     |dx(k+1)|    m - 1
     |-------| ~ -------
     | dx(k) |      m

For k somewhere in the middle, we obtain

     |dx(k+1)|    1.88e-03                    1
     |-------| = ---------- = 7.52e-01 ~ 1 - ---
     | dx(k) |    2.50e-03                    m

Solving for m leads to m = 4.03. Rounding to the nearest natural number leads to m = 4.

The following fixed-point iterations all have the same fixed point, i.e.: 5,
1. x(k+1) = 1/5 x(k)^2, k=0,1,...
2. x(k+1) = sqrt(5x(k)), k=0,1,...
3. x(k+1) = x(k)^2/(2x(k)-5), k=0,1,...
For each of the three fixed-point iterations, make a cobweb picture illustrating the convergence (or divergence), starting at x_0 = 4. Compute the convergence (or divergence) rate for each iteration. Which iteration is best?
Answers:
1. g(x) = 1/5 x^2, g'(x) = 2/5 x, at x=5: g'(5) = 2 > 1, thus divergence.
2. g(x) = sqrt(5x), g'(x) = 5/(2 sqrt(5 x)), at x=5: g'(5) = 1/2, linear convergence.
3. g(x) = x^2/(2 x-5), g'(x) = 2 x (x - 5)/(2 x - 5)^2, at x=5: g'(5) = 0, which means that the convergence is at least quadratic. The convergence is not higher than quadratic, because the method was derived from Newton's method on the equation f(x) = x(x-5). In general, one would have to evaluate also g''(5) to verify that it is different from zero.
So the third method has the fastest convergence. Follow the link for cobweb pictures.

Consider the linear system

                2.000 x_1 + 5.000 x_2 + 3.000 x_3 = 10.00 
                2.000 x_1 + 8.000 x_2 + 1.000 x_3 = 11.00 
                3.000 x_1 + 4.000 x_2 + 3.000 x_3 = 10.00

Use Gaussian elimination to compute an LU decomposition in the following four ways:
1. without partial pivoting and without scaling;
2. with partial pivoting but without scaling;
3. with explicit scaling and partial pivoting.
Solve the system three times, using the three LU decompositions obtained from above.
Compute the determinant using the three LU decompositions obtained from above.
Use the first LU decomposition you obtained to compute A^(-1) and cond(A) with ||.||_1.

Partial Answers:

Gaussian elimination:

without partial pivoting and without scaling:


       R2 = R2 - R1            [  2.000   5.000   3.000 ]
  A  ----------------------->  [  1.000   3.000  -2.000 ]
       R3 = R3 - (3/2)*R1      [  1.500  -3.500  -1.500 ]


       R3 = R3 + (3.5/3)*R2    [  2.000   5.000   3.000 ]
     ----------------------->  [  1.000   3.000  -2.000 ]
                               [  1.500  -1.167  -3.833 ]


                      [  1.000   0.000   0.000 ]
Thus A = L*U with L = [  1.000   1.000   0.000 ]
                      [  1.500  -1.167   1.000 ]

                      [  2.000   5.000   3.000 ]
              and U = [  0.000   3.000  -2.000 ]
                      [  0.000   0.000  -3.833 ]

with partial pivoting but without scaling:


              [  3.000   4.000  3.000 ] 3
  A  ---->    [  2.000   8.000  1.000 ] 2
              [  2.000   5.000  3.000 ] 1


       R2 = R2 - (2/3)*R1      [  3.000    4.000   3.000 ] 3
     ----------------------->  [  0.6667   5.333  -1.000 ] 2
       R3 = R3 - (2/3)*R1      [  0.6667   2.333   1.000 ] 1


       R3 = R3 - (2.333/5.333)*R2    [  3.000    4.000    3.000 ] 3
     ----------------------------->  [  0.6667   5.333   -1.000 ] 2 
                                     [  0.6667   0.4375   1.437 ] 1


                        [  0   0   1 ]
Thus A = P*L*U with P = [  0   1   0 ]
                        [  1   0   0 ]

                        [  1.000   0.000   0.000 ]
                    L = [  0.6667  1.000   0.000 ]
                        [  0.6667  0.4375  1.000 ]

                        [  3.000   4.000   3.000 ]
                and U = [  0.000   5.333  -1.000 ]
                        [  0.000   0.000   1.437 ]

Solving the system

with the first LU decomposition:

We solve two triangular systems:

 1) L y = b :   [  1.000   0.000   0.000 ] [ y1 ]   [ 10.00 ]
                [  1.000   1.000   0.000 ] [ y2 ] = [ 11.00 ]
                [  1.500  -1.167   1.000 ] [ y3 ]   [ 10.00 ]

     y1 = 10.00
     y2 = 11.00 - y1 = 1.000
     y3 = 10.00 + 1.167*y2 - 1.500*y1 = -3.836

 2) U x = y :   [  2.000   5.000   3.000 ] [ x1 ]   [ y1 ]
                [  0.000   3.000  -2.000 ] [ x2 ] = [ y2 ]
                [  0.000   0.000  -3.833 ] [ x3 ]   [ y3 ]

     x3 = -3.836/(-3.833) = 1.001
     x2 = (1.000 + 2.000*x3)/3.000 = 1.001
     x1 = (10.00 - 5.000*x2 - 3.000*x1)/2.000 = 0.996

with the second LU decomposition

We solve two triangular systems:

 1) P*L y = b or L y = P b :

                 [  1.000   0.000   0.000 ] [ y1 ]   [ 10.00 ]
                 [  0.6667  1.000   0.000 ] [ y2 ] = [ 11.00 ]
                 [  0.6667  0.4375  1.000 ] [ y3 ]   [ 10.00 ]

     y1 = 10.00
     y2 = 11.00 - 0.6667*y1 = 4.333 
     y3 = 10.00 - 0.6667*y1 - 0.4375*y2 = 1.437

 2) U x = y  :   [  3.000   4.000   3.000 ] [ x1 ]   [ y1 ]
                 [  0.000   5.333  -1.000 ] [ x2 ] = [ y2 ]
                 [  0.000   0.000   1.437 ] [ x3 ]   [ y3 ]

     x3 = 1.437/1.437 = 1.000
     x2 = (4.333 + 1.000*x3)/5.333 = 1.000
     x1 = (10.00 - 4.000*x2 - 3.000*x3)/3.000 = 1.000

Computing the determinant:

with the first LU decomposition:

    det(A) = det(L*U)

           = det(L)*det(U)

           = det(U)

                 [  2.000   5.000   3.000 ]
           = det [  0.000   3.000  -2.000 ]
                 [  0.000   0.000  -3.833 ]

           = 2.000*3.000*(-3.833)

           = -23.00

with the second LU decomposition:

    det(A) = det(P*L*U)

           = det(P)*det(L)*det(U)

                 [  0   0   1 ]
           = det [  0   1   0 ] * det(U)
                 [  1   0   0 ]

                   [  3.000   4.000   3.000 ]
           = - det [  0.000   5.333  -1.000 ]
                   [  0.000   0.000   1.437 ]

           = - 3.000*5.333*1.437

           = -22.99

Computing the inverse with the first LU decomposition (sketch):

   Denote the inverse of A by the unknown matrix X.  To compute X,
   we solve three linear systems, in matrix notation A*X = I, where
   I is the identity matrix, 

                       [ 1 0 0 ]
      I = [e1 e2 e3] = [ 0 1 0 ]
                       [ 0 0 1 ]
  
   With an LU decomposition, we perform three back substitutions:
   1) L y1 = e1, followed by U x1 = y1;
   2) L y2 = e2, followed by U x2 = y2;
   3) L y3 = e3, followed by U x3 = y3.
   Then the inverse consists of the columns [x1 x2 x3].

A matrix A, A^T = A, can be written as A = L D L^T, where D is a diagonal matrix and L a lower triangular matrix with ones on the diagonal. Derive the formulas to compute such a decomposition for a general 4-by-4 matrix.
Answers:
- Maple html generated output
- Maple worksheet

A tridiagonal matrix A has a special LU decomposition

       [  a_1  b_1                      ]
       [  c_2  a_2  b_2                 ]
       [       c_3  a_3  b_3            ]
       [            c_4  a_4  b_4       ]
       [                      c_5  a_5  ]

       [   1                     ] [ u_11  u_12                   ]
       [  l_1   1                ] [       u_22  u_32             ]
   =   [       l_2   1           ] [             u_33  u_34       ]
       [            l_3   1      ] [                   u_44  u_45 ]
       [                 l_4  1  ] [                         u_55 ]

assuming none of the diagonal elements of A will every become zero. Derive the formulas to compute this decomposition. Write the algorithm for general tridiagonal matrices like A.

Answers:

Please note the policy on skipping exams:

If an exam is missed, then greater weight will be placed on the final exam, especially on the material covered on the missing exam.

What this means is that if you decide not to take one midterm exam, your final exam will be weighted for one hundred points more.

What it does NOT mean is that you can drop the score of a midterm exam. If you take the midterm, then your score counts. So, please be prepared when you show up for the exam.