Unformatted text preview: 12.2
Suppose nvector of data is D and the n n Vandermonde matrix of x j is Vx . 1 x1 d0 d 1 x2
1
D
Vx d n1 1 xn x1n1 n
x2 1 n
xn 1 Then the degree n 1 polynomial interpolant of the data is VxC D of which c0 c C 1 . cn1 Since A maps D to an mvector of sampled values { p( y j )} , we have AD Vy C
n 1 y1 1 y2
Vy 1 ym
1 y1n1 n y2 1 m n1 ym 1 From above, we can attain that C Vx D and A Vy CD . Therefore we have A VyVx 1 . (b)
10 Semilog of Infinity Norm of A 10 5 10 0 10 0 10 20 30 (c)
With the condition, we have D 1 1
column of Vx is 1 1 1 . Since VxC D and the first
T 1 , it can be concluded that except the first element
T of C being 1, others are zero. In particular, C 1 0 1 and D 1 1 AD Vy C 1 1 T 1 and F ( D) ( AD1 )
D1
J ( D) T 1 . If we regard AD as T F ( D) , which is AD F ( D) , we have D 0 . Then 0 AD1
D1 ( AD2 ) D2 AD2 D2 ( ADn )
Dn ADn
Dn A . 1. In addition, According to 12.6, Therefore, J ( x)
f ( x) / x . Then J ( D)
F ( D) A VyVx 1 . The condition numbers
n=1 1 n=2 1 n=3 1.25000000000000 n=4 1.62500000000000 n=5 2.17187500000000 n=6 2.99218750000000 n=7 4.26367187500000 n=8 6.29394531250000 n=9 9.61932373046876 n=10 15.1834411621094
n=11 24.6609878540046
n=12 41.0473136901850
n=13 69.7373995780222
n=14 120.509180783711
n=15 211.184145869046
n=16 374.409046122295
n=17 670.263206750751
n=18 1209.77020252113
n=19 2198.87390756793
n=20 4020.91399819004 for each n are as following J ( D) /D A . n=21 7391.69458492794
n=22 13651.7216227416
n=23 25318.1403960481
n=24 47129.2705501207
n=25 88025.1241302036
n=26 164909.552835531
n=27 309807.653886070
n=28 583512.246276936
n=29 1101542.90295973
n=30 2084427.64366012 (d) Figure 11.1 with x linearly spaced from 1 to 1
6
4
2
0
2
1 0.5 0 0.5 1 From (c) we know that when n 11, 24.66 . In Figure 11.1, the infinity
condition number is 63657.407 which is much higher than . Therefore the result
is not close to the implicit bound. Complete code is as follows
clc;
clear;
close all;
n=30;
infnA=zeros(1,30);
f or n=1:30
m=2*n1;
x=linspace(1,1,n);
y=linspace(1,1,m);
Vx=fliplr(vander(x'));
Vy=fliplr(vander(y'));
Vy=Vy(:,1:n);
A=Vy/Vx;
if n==11
x0=x;
Vx0=Vx;
A0=A;
end
infnA(n)=norm(A,Inf);
end
semilogy(infnA);
title('Semilog of Infinity Norm of A');
grid on;
z=[0,0,0,1,1,1,0,0,0,0,0];
p=polyfit(x0,z,10);
x1=linspace(1,1,1000);
y1=polyval(p,x1);
figure;
plot(x0,z,'bx',x1,y1,'r');
title('Figure 11.1 with x linearly spaced from 1 to 1');
grid on;
f printf('Bound of Fig 11.1 is %f\n',cond(Vx0,Inf)); 13.2
(a)
First it is sensible to shed light on the formula 13.2 that is x (m / ) . Since
t the sign decides whether x is positive or negative and e e decides on the range of x, we can drop them and only observe m / . From the lecture, we know that
t 1 m t or equivalently t 1 m t 1. Therefore, it can be concluded that t and t 1 belong to F. But t +1 does not belong to F because it needs t+1 %Part A
for k = 1:5;
[U,R1] = qr(randn(50));
[V,R1] = qr(randn(50));
S = diag(sort(rand(50,1), 'descend'));
A = U*S*V';
[U2,S2,V2] = svd(A);
[norm(UU2) norm(VV2) norm(AU2*S2*V2')]
end
% norms of UU2 and VV2 are about 2, since the columns are
% nearly the same up to a factor of 1. Thus an SVD [u,s,v] of UU2
% will be essentially u = U2, s is diagonal with entries nearly 2 or
nearly 0,
% and v = identity. Hence the 2norm of UU2, i.e. the largest entry in
s, is about 2.
%Part B
for k = 1:5;
[U,R1] = qr(randn(50));
[V,R1] = qr(randn(50));
S = diag(sort(rand(50,1), 'descend'));
A = U*S*V';
[U2,S2,V2] = svd(A);
for k1 = 1:50;
if norm(U2(:,k1)  U(:,k1)) > 1.95
U2(:,k1) = U2(:,k1);
V2(:,k1) = V2(:,k1);
end
end
[norm(UU2) norm(VV2) norm(AU2*S2*V2') cond(A)]
end
% now errors are much smaller. Don't see a clear correlation with cond
(A)
for k = 1:5;
lim1 = 10*rand;
[U,R1] = qr(randn(50));
[V,R1] = qr(randn(50));
S = diag(sort(10.^(lim1*rand(50,1)), 'descend'));
A = U*S*V';
[U2,S2,V2] = svd(A);
for k1 = 1:50;
if norm(U2(:,k1)  U(:,k1)) > 1.95
U2(:,k1) = U2(:,k1);
V2(:,k1) = V2(:,k1);
end
end
[norm(UU2) norm(VV2) norm(AU2*S2*V2') cond(A)]
end % Now cond(A) varies over 10 orders of magnitude,
% and it seems that the errors are proportional to cond(A)
%Part C
for k = 1:5;
[U,R1] = qr(randn(50));
[V,R1] = qr(randn(50));
S = diag(sort(rand(50,1).^6, 'descend'));
A = U*S*V';
[U2,S2,V2] = svd(A);
for k1 = 1:50;
if norm(U2(:,k1)  U(:,k1)) > 1.9
U2(:,k1) = U2(:,k1);
V2(:,k1) = V2(:,k1);
end
end
[norm(UU2) norm(VV2) norm(AU2*S2*V2') cond(A)]
end
% Now when A is illconditioned, see much larger errors in U2 and V2
than
% in the product U2*S2*V2'. Numerical Linear Algebra: Homework 4
Due on Oct. 28, 2009 Da Kuang [email protected] 18.1
(a) A∗ A =
(A∗ A)−1 = 3
3.0002
3.0002 3.00040002
150020001 −150010000
−150010000 150000000 A+ = (A∗ A)−1 A∗
= 10001 −5000 −5000
−10000 5000
5000 P = AA+ 10
= 0 0.5
0 0.5 0
0.5
0.5 (b) x = A+ b = 1
1 2
y = Ax = 2.0001
2.0001 (c) Da Kuang [email protected] Homework 4 18.1 (continued) κ(A) = 4.2429 × 104
θ = 39.2306o
η=1 (d) b
A y
1.2910
5.4775 × 104 x
5.4775 × 104
1.4699 × 109 (e1)
For κb→y , it is attained when P (δb) = P
δ b . The right singular vector of P with the largest
∗
−14
singular value is v1 = (1, 0, 0) . Set δb = 10 v1 , then computing in MATLAB with exact P yields:
P δb
Pb δb
= 1.2934
b (e2)
For κb→x , it is attained when A+ (δb) = A+ δ b . The right singular vector of A+ with the largest
singular value is v1 = (−0.8165, 0.4082, 0.4082)∗ . Set δb = 10−14 v1 , then computing in MATLAB with exact
A+ yields:
δb
A+ δb
= 5.4239 × 104
A+ b
b
(e3)
∗
For κA→y , it is most likely to be attained when δA = (δp)vn , where δp is orthogonal to range(A).
Compute the full SVD of A: −0.5773 0.8165
0
2.4496
0
−0.7071 −0.7071
A = U ΣV ∗ = −0.5773 −0.4082 −0.7071 0
5.7733 × 10−5 0.7071 −0.7071
−0.5773 −0.4082 0.7071
0
0
The right singular vector of A with the smallest singular value is v2 = (0.7071, −0.7071)∗ . The vector
orthogonal to range(A) is u3 = (0, −0.7071, 0.7071)∗ . Then 0
0
∗
u3 v2 = −0.5 0.5 0.5 −0.5 ∗
Set δA = 10−14 u3 v2 , then computing in MATLAB yields: (A + δA)(A + δA)+ b
AA+ b δA
= 3.4552 × 104
A Since for κA→y , 5.4775 × 104 is just an upper bound, this ratio of relative change 3.4552 × 104 is acceptable. 18.1 continued on next page. . . Page 2 of 5 Da Kuang [email protected] Homework 4 18.1 (continued) (e4)
For κA→x , it is most likely to be attained when δA = δA1 + δA2 , where δA2 is deﬁned to be δA in (e3).
˜
˜
˜
˜
˜
˜
For δA1 , ideally it should satisfy A−1 (δ A)x = A−1 δ A x , where A and δ A are the version of A
and δA after change of bases. Since δA2 counts for the dominant part in the expression of upper bound of
condition number κA→x , we can simply take δA as δA2 , i.e. as δA in (e3). Computing in MATLAB yields:
(A + δA)+ b
A+ b δA
= 1.4660 × 109
A 20.5
(a)
The factorization A = LU takes the form that L is lowertriangular, and U is uppertriangular with all
1s on the diagonal.
(b)
Suppose A = LU , then AD = L(U D), so the factorization takes the form that L does not change, but
the columns of U are rescaled. The solution is D−1 x0 in this case, where x0 is the solution of Ax = b.
(c)
The ﬁnal factorization takes the form that A = LDU , where L is lowertriangular with all 1s on the
diagonal, D is diagonal, and U is uppertriangular with all 1s on the diagonal. 22.1
We can prove by induction.
m = 1: A = [a11 ] = 1 · [u11 ] = LU , where a11 = u11 . In this case, ρ = 1 = 21−1 .
m > 1: In the ﬁrst step of LU factorization, before the factorization we switch the row r with the ﬁrst
row, where ar1 = max1≤i≤m . Let A(1) be the resulting matrix after the ﬁrst step factorization. Then for
any 1 ≤ r, s ≤ m,
a(1)  ≤ ars  + a1s  ≤ 2 max aij 
rs
i,j Therefore
(1) max aij  ≤ 2 max aij 
i,j i,j After this step, we are left with a submatrix with dimention m − 1 to be factorized. Since ρm−1 ≤ 2m−2 22.1 continued on next page. . . Page 3 of 5 ...
View
Full
Document
 Fall '08
 Staff
 Linear Algebra, Algebra

Click to edit the document details