MATLAB QUESTION================================= for n=2:12, use GEPiv, LTriSol,

Question

MATLAB QUESTION=================================

for n=2:12, use GEPiv, LTriSol, and UTriSol to compute d(1:n), the diagonal of the inverse of the n-by-n Hilbert matrix. The build-in function hilb(n) can be used to set up those matrices. The exact inverse is obtainable via invhilb(n). Print a table that reports

for each n-value and the condition number of hilb(n). Submit output and the script used to produce it.

I believe all three GEPiv LTriSol UTriSol are matlab functions.

GEPiv is General LU Factorization with pivoting.

LTriSol is lower triangle solution, solving for Lx = b.

UTriSol is UpperTriangle Solution, solving for Ux = b

Explanation / Answer

close all
z = linspace(.25,1);
fz = sqrt(z);
for m = [2 100 ]
x = linspace(.25,1,m)';
A = [ones(m,1) x];
b = sqrt(x);
xLS = A;
alpha = xLS(1);
beta = xLS(2);
figure
plot(z,fz,z,alpha+beta*z,'--')
title(sprintf('m = %2.0f, alpha = %10.6f, beta = %10.6f',m,alpha,beta))
end
m = input('Enter m: ');
n = input('Enter n: ');
logcond = input('Enter the log of the desired condition number: ');
condA = 10^logcond;
[Q1,R1] = qr(randn(m,m));
[Q2,R2] = qr(randn(n,n));
A = Q1(:,1:n)*diag(linspace(1,1/condA,n))*Q2;
clc
disp(sprintf('m = %2.0f, n = %2.0f, cond(A) = %5.3e',m,n,condA))
b = A*ones(n,1);
[x,res] = LSq(A,b);
disp(' ')
disp('Nonzero residual problem.')
disp(' ')
disp(' Exact Solution Computed Solution')
disp('-----------------------------------------')
for i=1:n
disp(sprintf(' %20.16f %20.16f',1,x(i)))
end
b = b+Q1(:,m);
[x,res] = LSq(A,b);
disp(' ')
disp(' ')
disp(' ')
disp('Zero residual problem.')
disp(' ')
disp(' Exact Solution Computed Solution')
disp('-----------------------------------------')
for i=1:n
disp(sprintf(' %20.16f %20.16f',1,x(i)))
end

clc
n = 50;
disp(' m n full A flops sparse A flops')
disp('----------------------------------------')
for m = 100:100:1000
A = tril(triu(rand(m,m),-5),5);
p = m/n;
A = A(:,1:p:m);
A_sparse = sparse(A);
b = rand(m,1);
% Solve an m-by-n LS problem where the A matrix has about
% 10 nonzeros per column. In column j these nonzero entries
% are more or less A(j*m/n+k,j), k=-5:5.
flops(0)
x = A;   
f1 = flops;
flops(0)
x = A_sparse;
f2 = flops;
disp(sprintf('%4d %4d %10d %10d ',m,n,f1,f2))
end

m = 5;
n = 3;
disp(' ')
disp('Show upper triangularization of ')
A = rand(m,n);
for j=1:n
for i=m:-1:j+1
input('Strike Any Key to Continue');
clc
%Zero A(i,j)
Before = A
[c,s] = Rotate(A(i-1,j),A(i,j));
A(i-1:i,j:n) = [c s; -s c] *A(i-1:i,j:n);
disp(sprintf('Zero A(%g,%g)',i,j))
After = A
end
end
disp('Upper Triangularization Complete.')

close all
n = 100;
x = linspace(0,pi,n)';
hx = pi/(n-1);
d = 2*ones(n-2,1) + hx^2;
e = -ones(n-2,1);
[g,h] = CholTrid(d,e);
b = hx^2*( 2*x(2:n-1).*sin(x(2:n-1)) - 2*cos(x(2:n-1)));
umid = CholTridSol(g,h,b);
u = [0;umid;0];
plot(x,u,x,x.*sin(x))
err = norm(u - x.*sin(x),'inf');
title('Solution to -D^2 u + u = 2xsin(x) - 2cos(x), u(0)=u(pi)=0')
xlabel(sprintf(' n = %3.0f norm(u - xsin(x),''inf'') = %10.6f',n,err))

clc
n = input('Enter n: ');
X = rand(2*n,n);
A = X'*X;
disp(' ');
disp(sprintf('n = %2.0f',n))
disp(' ');
disp('Algorithm Relative Time Flops')
disp('-----------------------------------')

flops(0); tic; G = CholScalar(A); t1 = toc; f = flops;
disp(sprintf('CholScalar %6.3f %5.0f',t1/t1,f));

flops(0); tic; G = CholDot(A); t2 = toc; f = flops;
disp(sprintf('CholDot %6.3f %5.0f',t2/t1,f));

flops(0); tic; G = CholSax(A); t3 = toc; f = flops;
disp(sprintf('CholSax %6.3f %5.0f',t3/t1,f));

flops(0); tic; G = CholGax(A); t4 = toc; f = flops;
disp(sprintf('CholGax %6.3f %5.0f',t4/t1,f));

flops(0); tic; G = CholRecur(A); t5 = toc; f = flops;
disp(sprintf('CholRecur %6.3f %5.0f',t5/t1,f));

disp(' ')
disp('Relative time = Time/CholScalar Time')
clc
disp(' n cond(A) relerr ')
disp('----------------------------------')
for n = 2:12
A = hilb(n);
b = randn(n,1);
G = CholScalar(A);
y = LTriSol(G,b);
x = UTriSol(G',y);
condA = cond(A);
xTrue = invhilb(n)*b;
relerr = norm(x-xTrue)/norm(xTrue);
disp(sprintf(' %2.0f %10.3e %10.3e ',n,condA,relerr))
end

[n,n] = size(A);
G = zeros(n,n);
for i=1:n
% Compute G(i,1:i)
for j=1:i
if j==1
s = A(j,i);
else
s = A(j,i) - G(j,1:j-1)*G(i,1:j-1)';
end
if j < i
G(i,j) = s/G(j,j);
else
G(i,i) = sqrt(s);
end
end
end

[n,n] = size(A);
m = n/p;
A = MakeBlock(A,p);
G = cell(m,m);
for i=1:m
for j=1:i
S = A{i,j};
for k=1:j-1
S = S - G{i,k}*G{j,k}';
end
if j < i
G{i,j} = (G{j,j}S')';
else
G{i,i} = CholScalar(S);
end
end
end

G = MakeScalar(G);

[n,n] = size(A);
piv = 1:n;
for k=1:n-1
[maxc,r] = max(abs(A(k:n,k)));
q = r+k-1;
piv([k q]) = piv([q k]);
A([k q],:) = A([q k],:);
if A(k,k) ~=0
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end
end
L = eye(n,n) + tril(A,-1);
U = triu(A);

>> A = [17 24 1 8 15; 23 5 7 14 16; 4 6 13 20 22; 10 12 19 21 3; ...
2 11 18 25 2 9]
3
4 A =
5
6 17 24 1 8 15
7 23 5 7 14 16
8 4 6 13 20 22
9 10 12 19 21 3
10 11 18 25 2 9
11
12 >> b = [1; 3; 5; 7; 9]
13
14 b =
15
16 1
17 3
18 5
19 7
20 9
21
22 >> [L,U,piv] = GePiv(A);
23 >> L
24
25 L =
26
27 1.0000 0 0 0 0
28 0.7391 1.0000 0 0 0
29 0.4783 0.7687 1.0000 0 0
30 0.1739 0.2527 0.5164 1.0000 0
31 0.4348 0.4839 0.7231 0.9231 1.0000
32
33 >> U
34
35 U =
36
37 23.0000 5.0000 7.0000 14.0000 16.0000
38 0 20.3043 -4.1739 -2.3478 3.1739
39 0 0 24.8608 -2.8908 -1.0921
40 0 0 0 19.6512 18.9793
41 0 0 0 0 -22.2222
42
43 >> piv
44
45 piv =
46
47 2 1 5 3 4
48
49 >> y = LTriSol(L,b(piv));
50 >> y
51
52 y =
53
54 3.0000
55 -1.2174
56 8.5011
57 0.3962
58 -0.2279
59
60 >> x = UTriSol(U,y);
61 >> x
62
63 x =
64
65 0.0103
66 0.0103
67 0.3436
68 0.0103
69 0.0103
70
71 >> c = A*x
72
73 c =
74
75 1.0000
76 3.0000
77 5.0000
78 7.0000
79 9.0000

Related Questions

Navigate

• Browse (All)
• Browse M
• Subjects

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