Matlab code please
Problem 7 (10 points] Consider solving the linear system Ax = b, where A e Rnxn and beR. Suppose that we have computed a numerical solution x1 x. The residual is r = b Axi = AX - Axi A(x x1). Denote dy = x X and solve Ad = r1. Then x2 = x1 + d and Ax2 = Axi + Ad = Axi + Ax Ax = Ax = b. That is, X2 solves Ax2 = b. However, in computer arithmetic, r2 = b AC2 is unlikely the zero vector, and we can repeat the above with X2 to compute x3 and so on. This is called iterative refinement. When computing the residuals in the same precision, there will be cancellations, and we have to evaluate them in higher precision. Let Aij = 1/(i+j 1), i, j = 1, ..., n. This is a Hilbert matrix. Let u be the n-vector of ones, u= (1,1, ... 1)?, and let b = Au. Then Ax = b has the exact solution u. For a given n, store this A in single precision, and denote this single precision matrix by As. Compue b = Agu in double precision and round it to single precision. Then compute the residual in double precision and iterate: %file iterrefin.m clear all; close all; %turn this on to see warning messages W = warning ('off','all'); n = 8; iterations = 100; As = hilb(n,'single'); = ones (n,1,'single'); b = double(As)*double(u); bs = single(b); u tic; x(:,1) = As\bs; for i=2:iterations r = b - double(As)*double(x(:,1-1)); r = single(r); d = As ; x(:,i) = x(:,1-1)+d; end toc; error = vecnorm (x-u).orm(u); semilogy (error,'--o') legend('$\/x-ull Alu\$','interpreter', 'latex') set (gca, 'FontSize', 14) xlabel('iteration') ylabel('error') To see the importance of evaluating the residual in higher precision, try this code with r = b double(As)*double(x(:,1-1)) replaced by bs-As*x(:,1-1). 1 = (a) (2 points] Try this code with n = 5,6,7,8. Why does the error in the solution not go below 10-8? Why does the error increase for n > 9? (b) (2 points] Modify this code so it solves the systems using the LU factorization of Ag, where this factorization is computed by Matlab's lu. Include your code in the PDF file and highlight the places in the code where the LU factorization is used. Run your modified code with n = 8. If the plots are different, can you explain why. Include the resulting plots in the PDF. (c) (4 points] Implement in Matlab the following functions function [L, U, p] = lupp(A) % [L,U,p] = lupp(A) computes the LU factorization of A % using scaled partial pivoting. % p is a permutation vector storing the row permulations % during the pivoting. function x = lusolve(L,U,b,p) % Solves Lux=b (p) by solving Ly=b(p) and then Ux=y. % L is unit lower triangular, U is upper triangular % p is a permutation vector obtained from lupp. Include the produced plots for n = 7, 8, 9 in the PDF file. Problem 7 (10 points] Consider solving the linear system Ax = b, where A e Rnxn and beR. Suppose that we have computed a numerical solution x1 x. The residual is r = b Axi = AX - Axi A(x x1). Denote dy = x X and solve Ad = r1. Then x2 = x1 + d and Ax2 = Axi + Ad = Axi + Ax Ax = Ax = b. That is, X2 solves Ax2 = b. However, in computer arithmetic, r2 = b AC2 is unlikely the zero vector, and we can repeat the above with X2 to compute x3 and so on. This is called iterative refinement. When computing the residuals in the same precision, there will be cancellations, and we have to evaluate them in higher precision. Let Aij = 1/(i+j 1), i, j = 1, ..., n. This is a Hilbert matrix. Let u be the n-vector of ones, u= (1,1, ... 1)?, and let b = Au. Then Ax = b has the exact solution u. For a given n, store this A in single precision, and denote this single precision matrix by As. Compue b = Agu in double precision and round it to single precision. Then compute the residual in double precision and iterate: %file iterrefin.m clear all; close all; %turn this on to see warning messages W = warning ('off','all'); n = 8; iterations = 100; As = hilb(n,'single'); = ones (n,1,'single'); b = double(As)*double(u); bs = single(b); u tic; x(:,1) = As\bs; for i=2:iterations r = b - double(As)*double(x(:,1-1)); r = single(r); d = As ; x(:,i) = x(:,1-1)+d; end toc; error = vecnorm (x-u).orm(u); semilogy (error,'--o') legend('$\/x-ull Alu\$','interpreter', 'latex') set (gca, 'FontSize', 14) xlabel('iteration') ylabel('error') To see the importance of evaluating the residual in higher precision, try this code with r = b double(As)*double(x(:,1-1)) replaced by bs-As*x(:,1-1). 1 = (a) (2 points] Try this code with n = 5,6,7,8. Why does the error in the solution not go below 10-8? Why does the error increase for n > 9? (b) (2 points] Modify this code so it solves the systems using the LU factorization of Ag, where this factorization is computed by Matlab's lu. Include your code in the PDF file and highlight the places in the code where the LU factorization is used. Run your modified code with n = 8. If the plots are different, can you explain why. Include the resulting plots in the PDF. (c) (4 points] Implement in Matlab the following functions function [L, U, p] = lupp(A) % [L,U,p] = lupp(A) computes the LU factorization of A % using scaled partial pivoting. % p is a permutation vector storing the row permulations % during the pivoting. function x = lusolve(L,U,b,p) % Solves Lux=b (p) by solving Ly=b(p) and then Ux=y. % L is unit lower triangular, U is upper triangular % p is a permutation vector obtained from lupp. Include the produced plots for n = 7, 8, 9 in the PDF file
