Problem(s) of the Day

| CS 421 Home | 

These problems are designed to get you to think about lecture content.

Nov 29. How would you choose c = cos(theta) and s = sin(theta) so that the diagonal entries in [c s;-s;c]*A are equal where A = [w x; y z] is given.

Nov 27. Consider the initial value problem

        r'(t)  = 2r - .01*r*w      r(0) = 100;

        w'(t) = -w + .01*r*w   w(0) = 50

Describe the computations that must be performed for a single step of (a) Euler's method, (b) 2nd order Runge Kutta method, (c) 2nd order Adams-Moulton method, and (d) 2nd order Adams-Bashforth method.

Nov 20. Let A and B be n-by-n matrices. How would you use ODE23 to solve

              y''(t) = Ay'(t) + By     y(0) = y0 (given vector), y'(0) = w0 (given vector)

Nov 17. Why might Euler's method be better than the second order Adams-Bashforth method for the problem

                                                            y'(t) = ( 1+sin(t*y) )^{5/2}

Nov 15.  Function f has a 2nd derivative that satisfies  |f"| <= M2 and a 4th derivative that satisfies |f""|<=M4. We wish to integrate f from a to b. Under what circumstances would the composite trapezoidal rule with equal spacing be more efficient than the composite simpson rule with equal spacing? FYI, when applied to an integral from c to d, the trapezoidal rule error is bounded by (1/12)*M2*(d-c)^3   while the Simpson error is (1/2880)*M4*(d-c)^5

Nov 13. Suppose u(t) is an n-vector for all scalars t and that it is available through an m-file. Assume A is a nonsingular n-by-n matrix. How would you use QUAD8 to integrate f(t) = u(t)'*inv(A)*u(t) from t=0 to t=1?

Nov 10 We want to compute an integral I. A quadrature rule Q1 satisfies I = Q1 + c*h^k and a quadrature rule Q2 satisfies I = Q2 + c*(h/4)^k  for some k. How would you estimate |I - Q2|?

Nov 8. Take a look at the adaptive piecewise linear interpolation code pwLAdapt. Given an upper bound M on the second derivative of the function f on [xL,xR],  find an upper bound on the number required f-evaluations.

Nov 6.  Find a polynomial p(x) = a + bx + cx^2 + dx^3 so that p interpolates the function f(x) = sin(x) and its derivative at the points x = 0 and x = pi.

Nov 3. In the Levenberg-Marquardt approach to nonlinear least squares, the xNext = xCurrent + s where the n-vector s solves a linear system of the form (J'*J+ lambda*I)s = -J'r where J is m-by-n and r is m-by-1. Assume that the SVD of J is available. Give an efficient method for computing  f(lambda) = (norm(s,2)^2.

Nov 1. B is nonsingular and symmetric. We solve the linear system Bs = -g where g is the gradient and B is the Hessian. Suppose s is NOT a descent direction. Show that there is a nonzero vector q so that (a) Bq = mu*q with mu<0 and (b)  q'*g is not zero.

Oct 30.   Line searches often utilize parabolic interpolation. Suppose the quadratic q(t) interpolates a function f(t) at t= 0, t=.5, and t = 1. Under what conditions will we be able to find a tstar, 0<tstar<=1 , so q(tstar)<f(0).

Oct 27.  Suppose t and y are column m-vectors. We wish to fit data points (t(i),y(i)), i=1:m with a function

                                                                          f(t) = a*exp(bt). 

Taking the least squares approach, the idea is to choose a and b so that norm(a*exp(bt) - y,2) is minimized. Show how this problem could be solved using the 1-D minimizer  fminbnd. Hint: for a fixed b, what is the best a? Suggest a heuristic for picking the initial bracketing interval.

Oct 25. Suppose p(x)= a(1)+a(2)x + a(3)x^2 + a(4)x^3 + x^4 has 4 real roots a<b<c<d. We apply golden section search with initial interval [L,R] . Under what conditions  will the process find a local minimum?

Oct 23.  Suppose F(x) is a function from Rn to Rn with the property that its Jacobian J(x) is tridiagonal for all x. Explain how one can generate a divided difference approximation of J with just four F-evaluations. Hint: for a cleverly chosen vector v, (F(x+hv)-F(x))/h will provide an approximation to J(:,1), J(:,4), J(:,7),...

Oct 20. Consider applying the secant method to the function f(x) = (x^2 - A) where A>0. Show that

                        e(k+1)  =  e(k) * e(k -1) / (x(k) + x(k-1) )            where e(j) = x(j) - sqrt(A)

From this we conclude that  |e(k+1) / ( e(k) * e(k-1)) |  converges to a constant C>0.

Oct 18.  The act of finding sqrt(A) is the act of "building a square" with area A. If we have a rectangle with base x and height A/x, then we can make it "more square" by replacing x with the average of the "current" base and current height, i.e., xNew = (x + A/x)/2. This is precisely the update recipe obtained when we apply Newton's method to get a zero of the  function f(x) = x^2 - A.

To find cube roots we can apply Newton's method to   f(x) = x^3 -A   giving

                                                              xNew = x - (x^3-A)/(3x^2) =   ( 2x^3  +  A )/( 3x^2 )

Can you derive this recipe using an intuitive geometric argument similar to the one above for square roots? Start with the fact that we want to build a cube with volume A. Develop an update  that makes the "current cube" more cubical!

Oct 16.  Let A = [lambda_1  m ; 0 lambda_2] and assume that the eigenvalues lambda_1 and lambda_2 are distinct. Determine X so inv(A)*A*X = diag(lambda_1,lambda_2). Recall that if Ax = lambda* x and y'*A = lambda*y' then the condition of lambda is  1/abs(x'*y) assuming that x and y have unit 2-norm. What is the condition of lambda_1 and lambda_2 in the above 2-by-2 example? More generally, how would you compute the condition number of all the eigenvalues using [X,D] = eig(A) assuming that the eigenvalues of A are distinct?

Oct 13.  Assume that A has distinct eigenvalues and that [Q,T] = schur(A). How would you compute the condition of the eigenvalue T(k,k)?

Oct 11.  A is a real matrix and mu is a complex eigenvalue. Suppose Az = mu*z. Show that the real and imaginary parts of the (complex) eigenvector z define a real, 2-dimensional invariant subspace. Suppose the dominant eigenvalue lambda_1 of a real matrix is complex. Note that the conjugate of lambda_1 is also an eigenvalue so the dominant eigenvalue is not unique. What happens if we apply the power method with a real starting vector? Assume that all the other eigenvalues are smaller in maginitude.

Oct 6. Suppose A = [w x; y z] has complex eigenvalues alpha + i *beta and alpha - i*beta. Show how to compute a rotation Q = [c s; -s c] so that the diagonal entries of B = Q'*A*Q are equal. Relate the entries of B to alpha and beta.

Oct 4.  Suppose Ax = lambda x + r with norm(x,2) <= delta. Show that there is a matrix E with norm(E,2)<=tol so that (A+E)x = lambda*x. That is, lambda and x are an exact eigenpair of a matrix near to A. Hint: There is a rank-1 matrix E that works. You'll have to prove that norm(uv',2) = norm(u,2)*norm(v,2).

Oct 2.   Suppose A is symmetric  and that mu is very close to an eigenvalue lambda of A. Use the Schur decomposition to explain why  the solution to (A - mu*I)x = v is "rich" in the direction of the eigenvector associated with lambda.

Sept 29.  Determine c and s so that [c s ; -s c] is as close to [w x; y z] as possible in the Frobenius norm.

Sept 27.   The orthogonal complement of a subspace S is the the set of all vectors that are perpendicular to every vector in S. b is an m-vector and A is m-by-n with rank equal to n.  What is the closest vector to b that is in the orthogonal complement of the range of A? Give a solution using the SVD and give a solution using QR.

Sept 25.  A is m-by-n but rank(A) = m < n. A show how to obtain a solution to Ax = b given that you have the QR factorization of A' .

Sept 22. Reduce the number of required flops in this function from 3n to 2n by thinking carefully about the highlighted line:

   function v = House(x)
   % x a nonzero column n-vector
   % v is a unit column n-vector so that if P = I - wvv' then
   % Px  is a multiple of I(:,1)
    
   v = x;
   if x(1) > 0
      v(1) = x(1) + norm(x,2);
   else
      v(1) = x(1) - norm(x,2);
   end
   beta = norm(v);    % do this with O(1) flops!
   v = v/beta;
      

Sept 20   Here is how to solve the full rank min ||Ax - b|| problem using the QR factorization:

   [m,n] = size(A);
   [Q,R] = qr(A);
   x = R(1:n,1:n)\(Q(:,1:n)'*b);     % The minimizing x
   del = norm(Q(:,n+1:m)'*b,2)       % The minimum sum-of-squares

   function del = MinSumSq(A,b)
   % A is m-by-n with m>n = rank(A)
   % del is a column n-vector with del(k) the minimum value of
   % norm(A(:,1:k)z - b,2) where z is a column k-vector.

Sept 18. (a) Assume 0 <= theta < 2pi. Draw a picture that illustrates the rotation y = [c s;-s c]' x where                    c = cos(theta) and s = sin(theta). Draw a picture that illustrates the reflection y = [c s ; s -c]x.  (b) How would you solve the square, nonsingular linear system  A' x = b  given that you have the QR factorization A = QR?

Sept 15. Given column n-vectors x and y with distinct x(i),  write a Matlab script that determines scalars a and b so that y = ax+b is the least squares fitting line of the data (x(i),y(i)), i=1:n. Use the method of normal equations.

Sep 13. Suppose B is symmetric m-by-m and v is a column m-vector. Write a Matlab script that overwrites the lower triangular portion of B with the lower triangular portion of B - v*v'.

Sep 11. (a) Compute the Cholesky factor of A = [ 1 -2 3 ; -2 8 -6 ; 3 -6 31] by hand. (b) Show that if norm(F)<1 then I + F is positive definite.

Sep 8. Suppose M is n-by-n. What is the 1-norm condition of A = [eye(n,n) M ; zeros(n,n) eye(n,n) ].

Sep 6. (a) Suppose x and y are column n-vectors. The Cauchy-Schwartz inequality says

                                                              abs(x'*y) <= norm(x,2)*norm(y,2)

Show that norm(x,1) <= sqrt(n)*norm(x,2). The vector of all ones is useful!. (b) A is an n-by-n real matrix with         norm(A,1) = 1 and it is stored on a computer that supports IEEE double precision (p = 52). Suppose fl(A) = A + E. What is the largest possible value of norm(E,1)? Hmmm. that might be hard. Somebody send me the answer!

Sep 4. (a) Let p be the number of bits in the mantissa. If 127 is the nearest floating point number to 128, what is p? (b) The distance from 7 to the next largest floating point number is 2^{-12). What is the distance from 70 to the next largest floating point number? (c) What's the output when you run this:

   h = 1./2.;
   x = 2./3. - h;
   y = 3./5. - h;
   e = (x+x+x) - h;
   f = (y+y+y+y+y)-h;
   z = f/e;
   disp(sprintf('\nz = %10.5f',z))

% Plot (x-1)^6 = x^6 -6x^5 +15x^4 - 20x^3 + 15x^2 - 6x + 1 near x = 1...
x = linspace(.995,1.005);
y1 = (x-1).^6;
y2 = x.^6 - 6*x.^5 + 15*x.^4 - 20*x.^3 + 15*x.^2 - 6*x + 1;
subplot(2,1,1)
plot(x,y1)
subplot(2,1,2)
plot(x,y2)

Sep 1. (a) Let B be a 4-by-4 matrix to which we apply the following operations in succession:

  - double column 1 
  - halve row 3
  - interchange columns 1 and 4
  - subtract row 2 from each of the other rows
  - replace column 4 by column 3
  - delete column 1.

(b) Here is the LU factorization w/o pivoting:

     function [L,U] = OurLU(A)
   % A is an n-by-n matrix which has an LU factorization
   % A = LU where L is unit lower triangular and U is upper triangular.

   [n,n] = size(A);
   L = eye(n,n); m = zeros(n,1);
   for k=1:n-1
      % Generate the k-th Gauss Transform
      m(k+1:n) = A(k+1:n,k)/A(k,k);
      % Apply the k-th Gauss transform, i.e., update A(k+1:n,k+1:n), i.e.,
      % eliminate the k-th unknown from equations k+1 through n
      for i=k+1:n
         for j = k+1:n
            A(i,j) = A(i,j) - m(i)*A(k,j);
         end
      end
      L(k+1:n,k) = m(k+1:n);
   end
   U = A;

Aug 30 (a) How would you solve the linear systems Ax = b and A'z = c assuming that A=LU? (b) Suppose v is an n-vector with no zero components. How would you determine a vector m so that if M = I - m*ek' then Mv is zero everywhere except in component k? Here, ek = I(:,k).

Aug 28  Here is a lower triangular system solver:

   function y = LSolve(L,b)
% b is n-by-1 and L is an n-by-n nonsingular lower triangular matrix.
% y solves Ly = b

   n = length(b); y = zeros(n,1);
   for k=1:n
      y(k) = b(k)/L(k,k);
      b(k+1:n) = b(k+1:n) - y(k)*L(k+1:n,k);
   end

Aug 25. (a) A and B are n-by-n lower triangular matrices. Prove that C = AB is also lower triangular. (b) How many floating point operations are required to solve a 3-by-3 lower triangular system Ly = b?