7. (LU factorization) (a) write a program that takes the output A and p from pro

Question

7. (LU factorization) (a) write a program that takes the output A and p from problem # 4, along with a righthand side b, and computes the solution of Az -b by performing the forward and backward substitution steps. If you are using MATLAB, name your m-file lusolve.m. The first line of your code lusolve.m should be as follows: function x=1usolve (a , p ,b) (your code here!) Turn in a copy of your code. (b) The famous Hilbert matrices are given by Hy 1/(i+j - 1). The n x n Hilbert matrix Hn is easily produced in MATLAB using hilb(n). Assume the true solution of Hnx = b for a given n is x = [1, ,1. Hence the righthand side b is simply the row sums of Hn, and b is easily computed in MATLAB using b-sum(hilb(n))' Use your codes mylu.m and lusolve.m to solve the system Hn-b for n 5,10, 15, 20. For each n, using the oo- norm, compute the relative error and the relative residual. Discuss what is happening here. You may find it useful to look at the cond command in MATLAB.

Explanation / Answer

ANS (A) :
Here is mylusolve.txt:
type mylusolve.m
function x = mylusolve(a,p,b)
%
n = length(p);
x = zeros(n,1); z = zeros(n,1);
% forward sub to solve Lz=b
z(1) = b(p(1));
for i = 2:n
s = 0;
for j = 1:i-1
s = s+a(p(i),j)*z(j);
end
z(i) = b(p(i))-s;
end
% backsub to solve Ux=z
x(n) = z(n)/a(p(n),n);
for i = n-1:-1:1
s = 0;
for j = i+1:n
s = s+a(p(i),j)*x(j);
end
x(i) = (z(i)-s)/a(p(i),i);
end
a = rand(4); b = sum(a,2);
[a,p] = mylu(a);
x = mylusolve(a,p,b)
x =
1.0000
1.0000
1.0000
1.0000
norm(x-[1 1 1 1]’,’inf’)
ans =
8.8818e-16
diary off

ANS (B) :
Here is a code to do the computations and the output:
n = [5 10 15 20]’;
rel_err = zeros(4,1);
rel_resid = zeros(4,1);
cond_H = zeros(4,1);

for i = 1:length(n)
H = hilb(n(i)); % generate matrix
b = sum(H,2); % sum rows for RHS b
[H,p] = mylu(H); % solve system
x = mylusolve(H,p,b);
rel_err(i) = max(abs(ones(n(i),1)-x)); % note ||x||_inf = 1
H = hilb(n(i)); % we need the original H
rel_resid(i) = max(abs(b-H*x))/max(abs(b));
cond_H(i) = cond(H,’inf’);
end

disp(’ n rel_err rel_resid cond ’)
disp(’ -------------------------------------------------’)
disp([n rel_err rel_resid cond_H]);

Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 9.144066e-19.
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.291654e-19.

n rel_err rel_resid cond
-------------------------------------------------
5.0000e+00 5.5467e-13 1.9449e-16 9.4366e+05
1.0000e+01 1.3982e-04 1.5162e-16 3.5353e+13
1.5000e+01 2.4345e+01 4.0150e-16 1.0975e+18
2.0000e+01 4.7309e+01 3.7031e-16 9.2117e+18

So we see that while for each n the relative residual is on the order of machine precision, the
relative error is growing. This is directly attributable to the growth in the condition number
of the Hilbert matrices. In fact, I’ve included the warnings that MATLAB displayed for
n = 15, 20 to note that it is having trouble accurately computing Hn-1

Related Questions

Navigate

• Browse (All)
• Browse 7
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.