Math 365 Grady Wright Homework 4 Due Oct. 28, 2014 1. Functions with linear systems (20 points) Let T be a (diagonally dominant) tridiagonal matrix, A be a symmetric positive definite matrix, and B and C be full nonsingular matrices. Assume all of these matrices are of size n-by-n. Let f (x) be defined as follows f (x) = xT B −1 CT −1 A−1 x + bT B −T x, where x and b are column vectors of size n. Using the lu, linsolve, chol, spdiags, and mldivide (or \) functions show how you would efficiently evaluate the function f (x) for many different values of x in Matlab . Write some Matlab code showing exactly the sequence of steps you would use. Test your code using the following definitions for T , A, B, C, b, and x, and some value of n ≥ 10: T A A B C b x = = = = = = = spdiags(ones(n,1)*[1 -2 1],-1:1,n,n); tril(rand(n)); A*A’; rand(n); rand(n); rand(n,1); ones(n,1); Note: no where in your code should you actually compute the inverse of a matrix. 2. Conditioning of linear where  1/22  1/32  2 A=  1/42  1/5 1/62 systems (10 points) Consider the linear system of equation Ax = b,    1/32 1/42 1/52 1/62 0.10140  0.12479  1/42 1/52 1/62 1/72     2 2 2 2  0.096456  1/5 1/6 1/7 1/8  , b= .   0.071807  1/62 1/72 1/82 1/92  1/72 1/82 1/92 1/102 0.054118 Suppose the following approximate solution has been obtained in some way   −2.497  6.595     x=  −3.724  .  16.62  −15.81 It is then easy to show that the residual becomes  −1.1 · 10−5  −0.67 · 10−5  1.1 · 10−5 b − Ax =    −0.98 · 10−5 −1.0 · 10−5    .   (a) Does this imply that x is close to the exact solution x? Explain. (b) Use Matlab to obtain an accurate solution to the given system. (c) Use Matlab again to obtain a condition number for A. Use the appropriate result on perturbations of the right hand side (RHS) of a linear system to confirm that even though the residual is very small it is big enough to allow for the solution to be as far away from the correct one as occurs in this example. Hint: The A matrix can be constructed easily in Matlab using the function hankel (use Matlab help to see why). The following code constructs A: A = (1./hankel(2:6,6:10)).^2. The condition number (with respect to the two-norm) can be computed using the Matlab function cond: cond(A) 3. Conditioning of linear systems (25 points) For the heat conduction problem on the last homework assignment you solved a tridiagonal system of the form 2T1 − T2 −Tj+1 + 2Tj − Tj−1 −TN −2 + 2TN −1 = h2 g(x1 ) + a = h2 g(xj ), j = 2, . . . , N − 2 = h2 g(xN −1 ) + b where h = 1/N , xj = jh, and g(xj ) = 2f (xj ), j = 1, . . . , N − 1. This corresponds to a second order accurate finite-difference approximation of the differential equation −T 00 (x) = g(x), with boundary conditions T (0) = a and T (1) = b. If g(x) = 4k 2 π 2 sin(2kπx) then the exact solution to the differential equation is Texact (x) = sin(2kπx). In this problem, you will solve the above linear system for Tj and compare the solution to Texact (xj ) for various N and k. (a) Use spdiags to create the matrix corresponding to the linear system of equations above and call this matrix A. Use the backslash operator to solve the system and plot the solution for N = 100 and k = 4. On the same figure, plot the exact solution Texact (x). On a separate figure plot the error in the approximate and exact solutions. Add axes labels and a title to both plots. (b) Use condest to estimate the condition number of A for N = 100 · 2i , i = 0, 1, . . . , 13. Since you are storing the array as a sparse matrix you should be able to easily store (and solve) this matrix for these values of N . What happens to the conditioning of this system as N gets large? Plot a loglog plot of the condition number as a function of N . On the same graph, plot N 2 verses N . What does this second plot say about the scaling of the condition number with N ? (c) For each N in the range of values you used in part 3b and for k = 1, solve the linear system for the approximate T and compute the relative norm of the error as e(N ) = kT − Texact k∞ kTexact k∞ using the norm function with a second argument of inf. Plot this error on a loglog plot as a function of N . On the same plot, also include a plot of 1/N 2 , εκ(AN ), and 10−15 N 2 , where ε is machine epsilon (eps in Matlab ), and AN refers to the matrix from part (a) of size N -by-N . (d) What can you say about scaling of the error and the condition numbers from part 3c? Describe a potential downside to solving −T 00 (x) = g(x) using a (low-order) finite difference approximation. (e) (Extra credit: 5 points) Repeat part 3b for different values of k = 16. How does the resulting plot compare to part 3c. 2 4. Polynomial interpolation (30 points). As discussed in class, Lagrange’s interpolation formula can be rearranged into the magnificient barycentric formula (or more specifically the second (true) form of the barycentric formula): n P wj fj x − xj j=0 pn (x) = P , n wj j=0 x − xj where wj = Q n 1 (1) , j = 0, 1, . . . , n . (xj − xi ) i=0 i6=j Here xj are the given nodes and fj the corresponding function values. (a) Write a Matlab function for computing the barycentric weights {wj }nj=0 , for j = 0, . . . , n. Your function should take as input a vector containing the interpolation nodes {xj }nj=0 and return a vector containing the weights. Call your function barywghts and save it in agile called barywghts.m. Note that this can be done with 3 lines of code and no loops. Turn in a print out of your barywghts.m function and put in the DropBox folder. (b) Write a Matlab function for computing the barycentric interpolant pn (x) given by equation (1). Your function should take as input a vector containing the nodes xj , a vector containing the corresponding function values fj , and the location (or a vector of locations uj ) of where the interpolant should be evaluated. The output of the function should be the value of the interpolating polynomial at all the evaluation points. Call this function baryinterp and save it in a file called baryinterp.m. The function should begin with the line: function p = baryinterp(x,f,u) where x contains the nodes, f contains the function values, u contains the evaluation points. In this function you should first call the barywghts function from part (a) to generate the barycentric weights for the interpolating polynomial. Using these weights, you should write a loop that performs the operations in equation (1). Note that the Matlab function prod can be used to compute the product of all entries in a vector. Be careful to handle cases where you could get a divide by zero appropriately. Turn in a print out of your baryinterp.m function and put it in the Dropbox folder. (c) The three following node sets are known as equispaced, Chebyshev, and Legendre, respectively: 2j , j = 0, 1, . . . , 8 8 jπ (ii) xj = − cos 8 , j = 0, 1, . . . , 8 (iii) −0.96816023950763, −0.83603110732664, −0.61337143270059, −0.32425342340381, 0, 0.32425342340381, 0.61337143270059, 0.83603110732664, 0.96816023950763. (i) xj = −1 + (Ordering the Legendre points from smallest (-0.968...) to largest (0.968...) will make things easier). For each of these three node sets, use your baryinterp function to interpolate the 8th degree polynomial interpolant of the function f (x) = |x| + x/2 − x2 and evaluate it at 201 equally spaced points between [−1, 1]. The equally spaced points can be generated in Matlab using the command: u = linspace(-1,1,201); Plot the error, E(u) = p8 (u) − f (u), in the polynomial interpolant at these evaluation points for each of the three node sets (i.e. produce three plots). Which node set seems to produce the best result? What criteria did you use to determine what ‘best’ means? 3 5. Runge Phenomenon (15 points) (a) Do part (a) of NCM 3.9. (b) Replace line 62 of the NCM file rungeinterp.m with the line x = -cos(pi*(0:n-1)/(n-1)); This makes the distribution of interpolation nodes the Chebyshev distribution. Repeat part (a) with this distribution of points. (c) Change line 66 (v = polyinterp(x,y,u);) of the NCM file rungeinterp.m so that it instead uses your baryval.m function from problem 4. Repeat part 5(b) and verify the results have not changed (at least inside the interval [−1, 1]). Produce some kind of nice plot showing the results using the baryval function and the polyinterp function for n = 51. 4