{[ promptMessage ]}

Bookmark it

{[ promptMessage ]}

HW3Solns

# HW3Solns - 12.2 Suppose n-vector of data is D and the n n...

This preview shows page 1. Sign up to view the full content.

This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: 12.2 Suppose n-vector 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 m-vector 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*n-1; 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(U-U2) norm(V-V2) norm(A-U2*S2*V2')] end % norms of U-U2 and V-V2 are about 2, since the columns are % nearly the same up to a factor of -1. Thus an SVD [u,s,v] of U-U2 % will be essentially u = U2, s is diagonal with entries nearly 2 or nearly 0, % and v = identity. Hence the 2-norm of U-U2, 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(U-U2) norm(V-V2) norm(A-U2*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(U-U2) norm(V-V2) norm(A-U2*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(U-U2) norm(V-V2) norm(A-U2*S2*V2') cond(A)] end % Now when A is ill-conditioned, 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 lower-triangular, and U is upper-triangular 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 lower-triangular with all 1s on the diagonal, D is diagonal, and U is upper-triangular 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

{[ snackBarMessage ]}

Ask a homework question - tutors are online