CSC336-Assignment 2 Solved

1.        Consider the matrix A and its inverse A−1,

                                                                           6        13    −17               6      − 4     − 1 

                                                                 A =  13       29    −38   A−1 =  − 4      11         7 .

                                                                                                                                 

                                                                        −17    −38       50             − 1        7         5 

(a)What is the condition number of A in the infinity norm?

(b)Suppose we solve Ax = b for some b, and obtain xˆ, so that ||b − Axˆ||∞ ≤ 0. 01. How small an upper bound can be found for the absolute error ||x − xˆ||∞? Giv e the bound as a numerical value.

|| ˆx − x||∞

(c)With the same situation as in (b), how small an upper bound can be found for the relative error           ?

||x||∞ Give the bound in terms of ||b||∞.

 2. A group of n parachutists each with given mass mi and drag coefficient ci, i = 1,..., n, are connected by a weightless cord, and are falling at a given velocity v. We would like to calculate the tension ti, i = 1,..., n − 1, in each section of the cord and the acceleration a of the whole group. Let’s index the parachutists from top (i = 1) to bottom (i = n), and let g = 9. 81 be the acceleration of gravity. For the top parachutist (i = 1), Newton’s second law giv es the equation

t1 + m1g − c1v = m1a

which can be written as

                                                                                                                                                                                    t1 − m1a = c1v − m1g.                                                                                        (1)

For an arbitrary ‘‘interior’’ parachutist indexed i, i = 2, . . . , n − 1, Newton’s second law giv es the equation

−ti−1 + ti + mig − civ = mia

which can be written as

                                                                                                                                                                               −ti−1 + ti − mia = civ − mig                                                                                   (2)

For the bottom parachutist (i = n), Newton’s second law giv es the equation

−tn−1 + mng − cnv = mna

which can be written as

                                                                                                                                                                                 −tn−1 − mna = cnv − mng                                                                                     (3)

For convenience, let’s denote the unknown a by tn. Writing equations (1), (2) for i = 2,..., n − 1 and (3) civ (in that order), we get a linear system of equations At = b, with respect to the unknowns ti, i = 1,..., n. Note that the matrix of the system is very sparse; it has at most 3 non-zero entries per row, independently of the size of n. (However, it is not tridiagonal.)

(a) Write a MATLAB (or equivalent) script which, for n = 8, 16, 32, 64, generates the matrix and right-hand side vector of the linear system, then solves the linear system (using backslash), and outputs n, the maximum and minimum tension, the acceleration and the condition number of A, using format fprintf(’%3d %10.3e %10.3e %10.3e %10.3e0, .... After the loop of n, the script plots the tension vectors components (not including the acceleration), versus their normalized (by the respective n) index, in one plot (four lines plotted). Do this twice, (i)

                                                        i − 1                                 i − 1

           with v = 6, mi = 50 + 50  , and ci = 50 − 20    , and (ii) with v = 6, mi drawn randomly in the interval (50,100),

                                                        n − 1                                n − 1

CSC336                 Assignment 2       C. Christara - 2 -

then sorted from smallest to largest, and ci drawn randomly in the interval (30,50), then sorted from largest to smallest. Use appropriate labels and a legend. In latex, the two plots must be placed side-to-side with captions and subcaptions. At the end of the loop for n, and for the case (i) only, plot in log-log scale (loglog) the condition numbers versus n, using a solid line and thick dots for the data (’k.-’). Use appropriate labels.

Based on the numerical results (including plots), comment on how the acceleration and the maximum and minimum tensions behave with n. How do the components of the tension vectors vary with their index? Where (for which i) does the max tension occur? Also comment on how the condition numbers behave with n.

Also outside the loop for n, and for the case (i) and n = 16 only, plot the sparsity patterns of A, P, L,U, in a 2 × 2 format, either using latex or using subplot in matlab. Use appropriate titles and caption. Comment about whether they agree with what you expected. (These comments will be elaborated further in (b).)

Notes: For (i), use m = linspace(50, 100, n)’ and c = 50 - 20*linspace(0, 1, n)’; For (ii), use m = sort(50 + 50*rand(n, 1), ’ascend’); and c = sort(30 + 20*rand(n, 1),

’descend’);

Because the matrix A is sparse, we use sparse matrix techniques to generate it and store it. E.g e = ones(n, 1); A = spdiags([-e, e], [-1, 0], n, n); A(:, n) = -m; Note that you can visualize the sparsity pattern of a sparse matrix A by spy(A).

To get (an estimate of) the condition number of a sparse matrix A, use condest.

If you have four vectors of ni, i = 1,...,4, components respectively, to plot their components versus their normalized index use plot([1:ni(1)-1]/ni(1), t(1:ni(1)-1, 1), ’r-’, ...

[1:ni(2)-1]/ni(2), t(1:ni(2)-1, 2), ’g--’, ...

[1:ni(3)-1]/ni(3), t(1:ni(3)-1, 3), ’b-.’, ...

[1:ni(4)-1]/ni(4), t(1:ni(4)-1, 4), ’k.’);

(b)What would the forms of the L and U factors and of the permutation matrix P be, if LU factorization with row piv oting was applied to the matrix A? Your answer should be given in terms of n and mi (summation notation is acceptable). Note that this is a mathematical question, but MATLAB could help you get ideas.

(c)Find (mathematically) a closed form formula for the acceleration, and for the tensions ti, i = 1,..., n − 1, in terms of n, mi, ci, v and g. Justify mathematically where (for which i) the maximum tension occurs.

                                                                                                                                                          1                                ..., n, j = 1,..., n. The

3.                                                                                                                                              The Hilbert matrix of order n is defined to have entries (Hn)ij = , for i = 1, i + j − 1

Hilbert matrices are known to be very ill-conditioned once n increases beyond a certain value. Use MATLAB (or equivalent) to carry out the following experiments.

(a)        For n = 2,...,13 do the following:

− construct the Hilbert matrix Hn; you may use the built-in MATLAB function hilb(n) or your own function;

                                                                                                   n               j

− construct the vector bn, such that (bn)i = Σ  , for i = 1,..., n; j=1 i + j − 1

− solve Hn xn = bn by Gauss elimination; use the backslash \ operator to obtain the solution (help slash); let xˆn be the computed solution; note that the exact solution is xn = (1, 2, 3,..., n)T ; let rn = bn − Hxˆn be the residual; − compute the condition number of Hn, the computational relative error ||xn − xˆn|| / ||xn||, the residual ||rn||, and ||bn||, where || ⋅ || is the same norm as that used for the condition number; use the built-in MATLAB functions cond() for computing the condition number of Hn and norm() for computing the norms of vectors (or matrices); (note that both cond() and norm() use Euclidean norms); store the condition number of Hn in the nth component of a vector c, and the respective relative error, norm of residual and of right side vector in the nth component of vectors e, r, v, respectively.

− output n and the four quantities computed above with format fprintf(’%2d %9.2e %9.2e %9.2e

%9.2e\n’, ...

− Comment on the output and possible warnings by matlab.

(b)     After exiting from the above for-loop,

− In one graph, plot, in semilogy scale, c (the condition number of Hn) versus n, e (||xn − xˆn|| / ||xn||) versus n, and c(n)r(n)/v(n) versus n. (Thus n will be in the x axis, and will be in regular scale, and the y axis will be in log scale.) To distinguish between the three lines in the plot a solid one for the condition numbers, a dashed line for the relative errors, and a dotted line for the third line. Use appropriate labels and a legend.

− Interpret what you observe in the graph, and discuss whether it agrees with the theory.

Note: You may assume that MATLAB uses double precision. In other language, you must use double precision.