Answers to the Review Questions of Chapters 3, 4, and 5

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

  1. Consider the polynomial p(x) = x^2 + x - 5.

    1. Construct the Newton form of p(x) by divided differences, using the points (x_i,p(x_i)), with x_i = i, for i=0,1,2,3.

    2. Explain why the last element f_0123 in the table of divided differences you constructed above is (or should have been) zero.

    3. Apply Neville's algorithm to evaluate the interpolating polynomial at 0.5.

    4. Approximate p(x) with the linear function that minimizes the squares of the errors, using the points (x_i,p(x_i)), with x_i = i, for i=0,1,2,3.

    Answers:

    1.       0  -5 
                    -3 + 5 
                    ------- = 2
                     1 - 0
                                    3 - 2
            1  -3                  ------- = 1
                     1 + 5          2 - 1
                    ------- = 3                    1 - 1
                     2 - 0                        ------- = 0
                                    4 - 2          3 - 2
            2   1                  ------- = 1
                     7 + 5          3 - 1
                    ------- = 4
                     3 - 0
      
            3   7
      
         The Newton form is
      
            f_0 + f_01*(x-x[0]) + f_012*(x-x[0])*(x-x[1])
                + f_0123*(x-x[0])*(x-x[1])*(x-x[2])
      
          = -5 + 2*(x-0) + 1*(x-0)*(x-1) + 0*(x-0)*(x-1)*(x-2)
          = -5 + 2*x + x*(x-1)
      

    2. The interpolating polynomial is unique. If f_0123 would not be zero, the interpolating polynomial would be of degree three. But this contradicts with the given polynomial, which is of course also interpolating, and of degree two. Thus we would have two different interpolating polynomials (one of degree two and one of degree three), which is impossible. Therefore, f_0123 must be zero.

    3.      0   -5
                     p
                      01
           1   -3           p
                     p       012
                      12           p
           2    1           p       0123
                     p       123
                      23
           3    7
      
      with
                  (1/2 - 1)*(-5) - (1/2 - 0)*(-3)
           p   = --------------------------------- = -4
            01                 0 - 1
      
                  (1/2 - 2)*(-3) - (1/2 - 1)*1
           p   = ------------------------------ = -5
            12                 1 - 2
      
                  (1/2 - 3)*1 - (1/2 - 2)*7
           p   = --------------------------- = -8
            23              2 - 3
      
                   (1/2 - 2)*(-4) - (1/2 - 0)*(-5)
           p    = --------------------------------- = -17/4
            012                 0 - 2
      
                   (1/2 - 3)*(-5) - (1/2 - 1)*(-8)
           p    = --------------------------------- = -17/4
            123                 1 - 3
      
                    (1/2 - 3)*(-17/4) - (1/2 - 0)*(-17/4)
           p     = --------------------------------------- = -17/4
            0123                    0 - 3
      
      
      Thus p(0.5) = -17/4. Notice that, since we are interpolating a second degree polynomial with four points, we already see the final result appearing twice in the next to last column.

    4. Let y = a*x + b denote the linear function with unknown coefficients a and b. To determine a and b we solve the overdetermined linear system:
              -5 = a*0 + b
              -3 = a*1 + b
               1 = a*2 + b
               7 = a*3 + b
      
      or in matrix notation:
      
             [ 0 1 ]         [ -5 ]
             [ 1 1 ] [ a ] = [ -3 ]
             [ 2 1 ] [ b ]   [  1 ]
             [ 3 1 ]         [  7 ]
      
                X    [ a ] = f
                     [ b ]
      
      To solve the system in the least squares sense,
      we set up the normal equations:
      
          T
         X *X = [ 0 1 2 3 ] [ 0 1 ] = [ 14  6 ]
                [ 1 1 1 1 ] [ 1 1 ]   [  6  4 ]
                            [ 2 1 ]
                            [ 3 1 ]
      
          T
         X *f = [ 0 1 2 3 ] [ -5 ] = [ 20 ]
                [ 1 1 1 1 ] [ -3 ]   [  0 ]
                            [  1 ]
                            [  7 ]
                                     T            T
      The solution to the system    X *X [ a ] = X *f
                                         [ b ]
      is a = 4 and b = -6.
      
      
      

  2. The Maclaurin expansion of e^(x^2) is
              2
             x          2        4        6         8      10
            e   =  1 + x  + 1/2 x  + 1/6 x  + 1/24 x  + O(x  )
    
    
    Use this Maclaurin expansion to construct a Padé approximation for e^(x^2) as a quotient of two quadrics.

    Answer:

    Our quotient of two quadrics takes the following form:

                            2
            a0 + a1 x + a2 x   
           -------------------
                            2
             1 + b1 x + b2 x
    
    To find the coefficients, we make the identification with the Maclaurin expansion, after multiplication with the denominator:
    
      (      2        4 ) (                2 )                   2
      ( 1 + x  + 1/2 x  )*( 1 + b1 x + b2 x  ) - a0 - a1 x - a2 x 
      (                 ) (                  )
    
    Expanding and requiring the coefficients of this fourth-degree polyomial to be zero gives five equations in the five unknowns:
        1 - a0 = 0
       b1 - a1 = 0
       b2 + 1 - a2 = 0
       b1 = 0
       b2 + 1/2 = 0
    
    So we find b2 = -1/2, b1 = 0, a2 = 1/2, and a1 = 0.

  3. The answer to the previous question is
                                                  2
                                         1 + 1/2 x
                                    r =  ----------
                                                  2
                                         1 - 1/2 x
    

    1. Compute a continued-fraction representation of r(x).

    2. Count the number of arithmetical operations you need to evaluate the continued-fraction representation of r(x). Compare this number with the number of arithmetical operations needed to evaluate r(x) if you would use the Horner form of the numerator and denominator of r(x).

    Answers:

    1. We divide numerator by the denominator:
                 2     |       2
            1/2 x  + 1 | -1/2 x  + 1
                 2     |--------------
         -( 1/2 x  - 1)|  -1 
        ---------------
                   + 2
      
      So we find that
                 2                 2
            1/2 x + 1 = -1 ( -1/2 x  + 1 ) + 2
      
      and we rewrite the quotient as
                         2                  4
           r = -1 + ------------ = -1 + ----------
                              2                2
                     1 - 1/2 x            2 - x
      
      This is the continued-fraction representation.

    2. If we write the Horner form of numerator and denominator, then we find
      
                 ( 1/2 x + 0)*x + 1
            r = --------------------
                 (-1/2 x + 0)*x + 1
      
      Evaluating the Horner form takes 4 multiplications, 2 additions (ignoring the + 0), and one division. Evaluating the continued-fraction representation takes only one multiplication, one addition, one subtraction, and one division.

  4. Consider a function whose values are tabulated below:
         +---------+----------------+
         |    x    |      f(x)      |
         ============================
         |  0.000  |  1.0000000000  |
         |  0.125  |  0.9921976672  |
         |  0.250  |  0.9689124217  |
         |  0.375  |  0.9305076219  |
         |  0.500  |  0.8775825619  |
         |  0.625  |  0.8109631195  |
         |  0.750  |  0.7316888689  |
         |  0.875  |  0.6409968582  |
         |  1.000  |  0.5403023059  |
         +---------+----------------+
    
    Suppose we are interested in the derivative of f(x) at x = 0.5.

    1. One way to compute the derivative would be to compute the interpolating polynomial through all these points and then to take its derivative for x = 0.5. Explain why numerically this is not such a good idea.

    2. Compute the most accurate approximation for f'(0.5). Estimate the accuracy in your answer.

    Answers:

    1. The interpolating polynomial is only a poor approximation for the derivative because of it oscillatory behaviour between the interpolation points. We can see this from a picture: including more interpolation points (using higher degree polynomials) will not necessarily increase the accuracy of the approximation.

    2. We take central difference formulas at x=0.5 with h = 0.5, 0.25, and 0.125. Then we extrapolate to obtain a sixth-order approximation for the derivative:
              -0.4596976941
      
              -0.4744471056    -0.4793635761
      
              -0.4781780096    -0.4794216443     -0.479425155
      
      The first column contains the approximations for h = 0.5,0.25, and 0.125, which have error of order O(h^2). The second column has an error of order O(h^4), while the error for the last number is of order O(h^6). As h = 0.125, we compute h^6 = (.125)^6 = 3.8E-6. So we expect to have about six decimal places correct.

      Note: The values in the table above were taken from sampling the cosine function. The derivative is -sin(0.5) = -0.4794255386, so we have indeed six decimal places right.

  5. Consider the operator D defined by Df(x) = df/dx, and the operator Delta defined by Delta f(x) = f(x+h) - f(x), for any h > 0. Show that D = ln(I+Delta)/h.

    Answer:

          f(x+h) = f(x) + f'(x)*h + f"(x)*h^2/2! + .. + f^(n)(x)*h^2/n! + ...
    
           Ef(x) = ( I + D*h + D^2*h^2/2! + .. + D^n*h^n/n! + .. ) f(x) 
    
                 = e^(h*D) f(x)  => E = e^(h*D)
    
          h*D = ln(E) and Delta = E - I, or E = I + Delta
    
          thus: D = ln(I+Delta)/h
    

  6. Suppose we interpolate f(x) with divided differences (i.e., with Newton interpolation) in equidistant points, at points x_i = x_0 + i*h, i=0,1,...,n. Write f[x_0,x_1,x_2] using the operator Delta.

    Answer:

                         f(x_1) - f(x_0)     Delta f(x_0)     1
           f[x_0,x_1] = ----------------- = -------------- = --- Delta f(x_0)
                            x_1 - x_0             h           h 
    
    
                             f[x_1,x_2] - f[x_0,x_1]
           f[x_0,x_1,x_2] = -------------------------
                                    x_2 - x_0 
    
    
                              Delta f(x_0+h)     Delta f(x_0)
                             ---------------- - --------------
                                    h                 h
                          = -----------------------------------
                                            2 h
    
                                  2
                             Delta f(x_0)
                          = --------------
                                    2
                                 2 h
    

  7. Explain why the formulas for Richardson extrapolation to approximate the derivative with forward differences are different from those that use central differences. In particular, why do we see only even powers of the ratio r in the extrapolation formulas that use central differences?

    Answer:

    The formula for the derivative with central differences is an even function in h. Thus the expansion of the error this formula makes with f'(x) contains only even powers of h. Therefore, if we eliminate terms in the error expansion by elimination, we only use even powers of r, when extrapolating with values for h and r*h. With Taylor series of f(x+h) we immediately get the error expansion for the forward difference formula, which contains in general all powers of h. So the extrapolation formulas contain also odd powers of r, when extrapolating with values for h and r*h.

  8. Consider
                 1
                 /  x      
                |  e  cos(2 Pi x) dx    
                / 
                 0
    

    1. Apply the composite trapezoidal rule using four subintervals of equal length to approximate this integral. Write your answer with six decimal places.

    2. Estimate the accuracy of the numerical approximation you just computed. How many decimal places can be correct?

    3. Apply Romberg integration to obtain a sixth-order approximation of the integral. Use six decimal places to denote the answers. What is the accuracy of your final answer?

    Answers:

    1. We calculate
              0.25*(1/2*e^0*cos(0) + e^(0.25)*cos(2*pi*0.25)
                                   + e^(0.5)*cos(2*pi*0.5)
                                   + 2^(0.75)*cos(2*pi*0.75)
                   +1/2*e^1*cos(2*pi))
      
              = 5.26049E-2
       
      

    2. The error of the trapezoidal rule is O(h^2). As h = 0.25, we compute (0.25)^4 = 0.0625. Observe that this estimate applies to the absolute error!

    3. The extrapolation table from Romberg integration is
                  5.26049E-2
      
                  4.47567E-2    4.21406E-2
      
                  4.30130E-2    4.24318E-2   4.24512E-2
      
      The error estimate is O(h^6), for h = 1/16, we find 5.98E-8. Also this estimate applies to the absolute error, dividing by the answer 4.2E-2, we find about six correct decimal places.

  9. Simpson's rule on an interval [a,b] approximates
                 b
                 / 
                |  f(x) dx 
                / 
                 a
    
    
    by
              b - a  (           ( a+b )        )
             ------- ( f(a) + 4 f(-----) + f(b) )
                3    (           (  2  )        )
    

    1. Write a composite formula to integrate
                   b
                   /        
                  |  f(x) dx    
                  / 
                   a
      
      with Simpon's rule, using seven function evaluations.

    2. Give the formula for the general composite Simpson's rule, over n subintervals of [a,b], of length h = (b-a)/n.

    Answers:

    1. If we have to carry out seven function evaluations, then this implies that we have consider dividing the interval [a,b] in three subintervals of equal length h = (b-a)/3. Applying Simpson's rule to each subinterval yields the following formula:
        h  [
       --- [ f(a) + 4 f(a+h/2) + 2 f(a+h) + 4 f(a+3h/2) + 2 f(a+2h)
        3  [      
                                       ]
                  + 4 f(a+5h/2) + f(b) ]
                                       ]
      

    2. With little effort we now can generalize the formula from 3 subintervals to n subintervals of length h = (b-a)/n:
                             n-1                   n-1
        h  [                 ---                   ---         ]
       --- [ f(a) + f(b) + 4 >   f(a+(i+1/2)h) + 2 >   f(a+ih) ]
        3  [                 ---                   ---         ]
                             i=0                   i=1
      

  10. Consider the approximation of
                 2a
                 /        
                |  f(x) dx    
                / 
                 0
    
    by the rule w_1 f(0) + w_2 f(a/2) + w_3 f(a).

    1. Determine the weights w_1, w_2, and w_3 so that the rule has the highest possible algebraic degree of precision.

    2. What is the highest possible algebraic degree of precision we can reach with three function evaluations? Give the nonlinear system in the weights w_1, w_2, w_3 and abscisses x_1, x_2, x_3 to determine the quadrature rule w_1 f(x_1) + w_2 f(x_2) + w_3 f(x_3).

    Answers:

    1. We want the equality
                   2a
                   /
                  |   f(x) dx = w_1 f(0) + w_2 f(a/2) + w_3 f(a)
                  /                               
                   0
      
      to hold for f = 1, x, and x^2 .  Therefore, we solve the system
      
                   2a
                   /
                  |   1 dx =   2 a  = w1  + w2 + w3
                  / 
                   0
      
                   2a
                   /              2
                  |   x dx =   2 a  = w1 0 + w2 a/2 + w3 a
                  / 
                   0
      
                   2a
                   /   2      8   3              2         2
                  |   x dx = --- a  = w1 0 + w2 a /4 + w3 a
                  /           3
                   0
      
      where a is a parameter.  The solution to the linear system is
      
               4                8             10
         w  = --- a     w  = - --- a    w  = ---- a
          1    3         2      3        3     3
      
      

    2. To define Gaussian quadrature formulas, we have to set up systems of nonlinear equations. Here we have six equations in six unknowns, of the form
                   2a          k+1    k+1
                   /   k      b    - a             k        k        k
                  |   x dx = -------------  = w1 x1  + w2 x2  + w3 x3
                  /               k+1
                   0
      
      for k from 0 to 5. If the weights and abscisses satisfy the system then every polynomial of degree five is integrated correctly, so the Gaussian quadrature rule has algebraic degree of precision equal to five.

  11. The Fourier series of a function f(t) appear in two different forms: either as
                          infinity
                          --------
                  1       \
         f(t) =  --- a   + >      a cos(2 Pi k t) + b sin(2 Pi k t)
                  2   0   /        k                 k
                          --------
                           k = 1
    
    or as
                 infinity
                 --------
                 \           i 2 Pi k t
         f(t) =   >      c e
                 /        k     
                 --------
               k = -infinity
    
    Derive the relationship between the coeffients c_k and (a_k, b_k).

    Answer:

    We first use the following relations:

             i 2 Pi k t
            e            = cos(2 Pi k t) + i sin(2 Pi k t)
    
    and
            -i 2 Pi k t
           e             = cos(2 Pi k t) - i sin(2 Pi k t)
    
    Adding the relations gives us
                          1  (  i 2 Pi k t    - i 2 Pi k t )
         cos(2 Pi k t) = --- ( e           + e             )
                          2  (                             )
    
    Subtracting the relations after multiplication by -i yields
                         -i  (   i 2 Pi k t    - i 2 Pi k t )
         sin(2 Pi k t) = --- (  e           - e             )
                          2  (                              )
    
    Now we will use these relations to transform the cosine/sine representation of the Fourier series:
               1
       f(t) = --- a  + a cos(2 Pi t) + b sin(2 Pi t) + ...
               2   0    1               1
    
               1        1     (  i 2 Pi t    - i 2 Pi t )
            = --- a  + --- a  ( e         + e           ) 
               2   0    2   1 (                         )
    
                       -i     (  i 2 Pi t    - i 2 Pi t )
                     + --- b  ( e         - e           ) + ...
                        2   1 (                         )
    
               1        1  (           )  i 2 Pi t
            = --- a  + --- ( a  - i b  ) e
               2   0    2  (  1      1 )
    
                        1  (           )  - i 2 Pi t
                     + --- ( a  + i b  ) e            + ...
                        2  (  1      1 )
    
    So we can see that
               
                    1  (           )              1  (          )
              c  = --- ( a  - i b  )  and  c   = --- ( a + i b  )
               1    2  (  1      1 )        -1    2  (  1     1 )
    
    The derivation is completely analogue for the k-th coefficient, we just have to pick the k-th term from the cosine/sine sum in the Fourier series of f(t) and perform the same transformations.

FINAL EXAM is in BH 208 on Wednesday 4 May 2005, from 8 till 10AM.