< M A T L A B > Copyright 1984-2001 The MathWorks, Inc. Version 6.1.0.450 Release 12.1 May 18 2001 To get started, select "MATLAB Help" from the Help menu. >> rand('state',sum(100*clock)); ----- Problem 1 >> A=rand(7,7); >> [L U] = lu(A); a) >> norm(A-L*U,1) ans = 3.0531e-16 % So we see that A = LU b) >> b=A*ones(7,1); >> y=L\b; >> x=U\y x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 % We know that [1 1 1 1 1 1 1]' is a solution by construction % and so x is as we expected. % Unfortunately the problem was not worded clearly since it mixed % mathematical notation with matlab notation. % Many of you took it literally and computed b as I do below. % I still gave full credit as long as you did something sensible. >> b=A([1 1 1 1 1 1 1]') b = 0.9246 0.9246 0.9246 0.9246 0.9246 0.9246 0.9246 % The above b is not A*[1 1 1 1 1 1 1]' % Instead it takes the first element of A (the upper right entry) % and repeats it seven times. >> A(1) ans = 0.9246 >> A A = 0.9246 0.4479 0.4334 0.2512 0.8541 0.3490 0.0141 0.2992 0.0139 0.8091 0.0843 0.5571 0.7954 0.1695 0.5658 0.0317 0.1657 0.4558 0.8730 0.2161 0.0660 0.5048 0.5196 0.4786 0.0291 0.3568 0.3789 0.6109 0.7541 0.9072 0.7182 0.0793 0.5272 0.9391 0.2312 0.1007 0.4888 0.3340 0.0456 0.8512 0.8123 0.9514 0.8028 0.7609 0.8006 0.1802 0.6039 0.5180 0.8915 c) >> [Lt Ut] = lu(A'); % Let's see if this LU decomposition of A' is U' L' >> norm(Lt-U') ans = 2.2419 >> norm(Ut-L') ans = 1.9327 % No it is not. let us investigate further by examining the matrices. >> Lt, U', Ut, L' Lt = 1.0000 0 0 0 0 0 0 0.4845 -0.1921 -0.6928 0.3089 -0.7958 -0.6394 1.0000 0.4688 0.9800 -0.2923 0.0046 1.0000 0 0 0.2717 0.0044 0.8647 0.0899 -0.5891 1.0000 0 0.9238 0.4113 1.0000 0 0 0 0 0.3775 1.0000 0 0 0 0 0 0.0152 0.2417 0.1625 1.0000 0 0 0 ans = 0.9246 0 0 0 0 0 0 0.4479 0.5418 0 0 0 0 0 0.4334 0.3647 0.7570 0 0 0 0 0.2512 -0.1256 -0.0274 0.2482 0 0 0 0.8541 -0.1695 0.2397 0.2544 0.7759 0 0 0.3490 0.6544 0.8408 0.2248 0.1446 -0.4522 0 0.0141 0.2197 0.2181 0.1374 0.7077 0.7689 0.2096 Ut = 0.9246 0.2992 0.5658 0.5048 0.7541 0.1007 0.8028 0 0.6824 0.0025 0.1883 0.6544 0.7743 0.2150 0 0 0.3493 -0.1870 -0.4386 0.4397 -0.2261 0 0 0 0.5881 0.1329 0.6913 0.8641 0 0 0 0 -0.4055 -0.3467 0.1436 0 0 0 0 0 -0.6318 0.1635 0 0 0 0 0 0 0.2085 ans = 1.0000 0.3236 0.6120 0.5460 0.8157 0.1089 0.8683 0 -0.2420 -0.4475 0.5076 1.0000 0.8122 0.6865 0 1.0000 0.0840 0.0751 0 -0.0124 0.2298 0 0 1.0000 -0.1702 0 0.4832 0.2198 0 0 0 0.0023 0 1.0000 -0.1706 0 0 0 0.3738 0 0 1.0000 0 0 0 1.0000 0 0 0 % So we see that, as was stated in the matlab #3 text % when matlab finds the LU decomposition it may do some row switches, % and then L will not be lower triangular, but will be a row switched lower triangular matrix. % On the other hand U is a genuine upper triangular matrix. % This is one source of the difference. % If matlab had done no row switches, then Lt and U' would both be lower triangular. % But they still would probably not be equal since Lt would have 1 on the diagonal, % but U' probably would not. ------- Problem 2 >> A=rand(8,8); B=rand(8,8); a) >> a=det(A*B); b=det(A)*det(B); [a b] ans = -0.0035 -0.0035 >> a-b ans = -2.1684e-18 % Since the difference is small compared to a and b, the identity seems to hold. b) >> a=det(inv(A)); b=1/det(A); [a b], a-b ans = 32.2968 32.2968 ans = 7.1054e-15 % Since the difference is small compared to a and b, the identity seems to hold. c) >> a = det(A+B); b=det(A)+det(B); [a b], a-b ans = -2.0251 -0.0826 ans = -1.9425 % a and b differ by quite a bit so the identity does not hold. d) >> Z = zeros(8,8); >> C=[A Z; Z B]; D=[inv(A) Z; Z inv(B)]; >> norm(inv(C)-D) ans = 5.6948e-16 >> norm(D) ans = 10.6925 % Since the difference is small compared to D, the identity seems to hold. e) >> a = det([A eye(8);Z B]); b = det(A)*det(B); [a b], a-b ans = -0.0035 -0.0035 ans = -1.3010e-18 % Since the difference is small compared to a and b, the identity seems to hold. f) >> a = det(A); b=prod(diag(A)); [a b], a-b ans = 0.0310 0.0007 ans = 0.0303 % The difference is small, but so are a and b. % In fact they are quite different so the identity does not hold. g) >> a = det(diag(diag(A))); [a b], a-b ans = 1.0e-03 * 0.6858 0.6858 ans = 0 % Since the difference is small compared to a and b, the identity seems to hold. --------- Problem 3 >> V=rand(10,10); U=eye(10)+1000*triu(V,1); >> det(U), det(U'), det(U*U'), det(U)*det(U') ans = 1 ans = 1.0000 ans = -6.3975e+33 ans = 1.0000 % Most versions of matlab gave an inaccurate answer for det(U') but this didn't. % However, all gave a wildly inaccurate answer for det(U*U'). % the problem is, since the diagonal entries which give the determinent are so much smaller % than the off-diagonal entries the roundoff error is exagerated.