Ma 4501: Numerical Applied Mathematics Model Solutions to Homework Assignment 3 Prof. Wickerhauser Exercise 1(iv) of Section 3.1, p.108. < M A T L A B > Version of 25 May 1982 Ported to NeXT 30 Nov 1991 HELP is available <>x=[1,-2,4,2]; y=[3,-5,-4,0]; <>x+y ANS = 4. -7. 0. 2. <>x-y ANS = -2. 3. 8. 2. <>3*x ANS = 3. -6. 12. 6. <>norm(x) ANS = 5.0000 <>7*y-4*x ANS = 17. -27. -44. -8. <>x*y' ANS = -3. <>norm(7*y-4*x) ANS = 54.9363 Exercise 4(c) of Section 3.1, p.108. <>A=[-1 9 4;2 -3 -6; 0 5 7] A = -1. 9. 4. 2. -3. -6. 0. 5. 7. <>B=[-4 9 2; 3 -5 7; 8 1 -6] B = -4. 9. 2. 3. -5. 7. 8. 1. -6. <>3*A-2*B ANS = 5. 9. 8. 0. 1. -32. -16. 13. 33. Exercise 5(b) of Section 3.1, p.108. <>B=[4 9 2; 3 5 7; 8 1 6] B = 4. 9. 2. 3. 5. 7. 8. 1. 6. <>b' ANS = 4. 3. 8. 9. 5. 1. 2. 7. 6. Exercise 2 of Section 3.2, p.118. <>A=[1 -2 3; 2 0 5]; B=[3 0;-1 5; 3 -2]; <>A*B = 14. -16. 21. -10. <>B*A = 3. -6. 9. 9. 2. 22. -1. -6. -1. Exercise 5(a*) of Section 3.2, p.118. Use Matlab: <> a=[-1 -7; 5 2]; det(a) ANS = 33. or else compute det(a)= (-1)(2)-(5)(-7)= -2+35 = 33. Exercise 5(b) of Section 3.2, p.118. <>A=[2 0 6; -1 5 -4; 3 -5 2]; det(A) = -80. Algorithm 2 of Section 3.2, p.119. Use help plot3 to get some documentation on plotting. The following command draws a polygon through the points (x(1),y(1),z(1)) ,..., (x(n),y(n),z(n)), where x,y,z are column vectors of length n: <> plot3(x,y,z) To get a cube, imagine tracing the edges in sequence and then put the vertex coordinates into x,y, and z in the order they are encountered. It is OK to retrace edges; the following solution with n=16 retraces 3 edges: % Vertices in the order they are encountered: u0 = [ 0 0 0; 1 0 0; 1 0 1; 1 1 1; 1 1 0; 0 1 0; 0 1 1; 0 0 1; 0 0 0; 0 1 0; 0 1 1; 1 1 1; 1 1 0; 1 0 0; 1 0 1; 0 0 1; 0 0 0] % Clear the current figure clf % Copy columns into x,y,z and plot the original cube with black lines: x=u0(:,1); y=u0(:,2); z=u0(:,3); plot3(x,y,z,'Color', 'black') % Set the axes big enough to see the cube in context: lims = [-2 2]; axis([lims lims lims]); axis manual; hold on % Make the rotation matrices alpha = pi/12; c=cos(alpha); s=sin(alpha); Rx=[1 0 0; 0 c -s; 0 s c]; beta = 0; c=cos(beta); s=sin(beta); Ry=[c 0 s; 0 1 0; -s 0 c]; gamma = pi/6; c=cos(gamma); s=sin(gamma); Rz=[c -s 0; s c 0; 0 0 1]; % Apply the rotation matrices to the original cube's vertices and plot % the results using three different colors: % NOTE: here we use the lemma: (R*u')' = u*R' u = u0*Rx'; x=u(:,1); y=u(:,2); z=u(:,3); plot3(x,y,z,'Color', 'red') u = u*Ry'; x=u(:,1); y=u(:,2); z=u(:,3); % plot3(x,y,z,'Color', 'green') % nothing new to plot u = u*Rz'; x=u(:,1); y=u(:,2); z=u(:,3); plot3(x,y,z,'Color', 'blue') Exercise 3 of Section 3.3, p.124. Use back substitution in Matlab: <> x5=6/3 X5 = 2. <> x4=(10+x5)/(-2) X4 = -6. <> x3=3+x4+2*x5 X3 = 1. <> x2=(0-6*x3-2*x4-7*x5)/(-2) X2 = 4. <> x1=(4+x2-2*x3-2*x4+x5)/4 X1 = 5. Exercise 6 of Section 3.3, p.124. <>A=[5 0 0 0;1 3 0 0;3 4 2 0;-1 3 -6 -1] = 5. 0. 0. 0. 1. 3. 0. 0. 3. 4. 2. 0. -1. 3. -6. -1. <>b=[-10 4 2 5]' = -10. 4. 2. 5. <>det(A) = -30.0000 <>inv(A)*b = -2.0000 2.0000 .0000 3.0000 Algorithm 2 of Section 3.3, p.125. MATLAB PROGRAM: Write the following into the file "forsub.m": % Expect lower triangular matrix A and vector b: function x = forsub(A, b) N = length(b); % expect NxN matrix A, Nx1 vector b for i=1:N x(i) = b(i); % start with i'th right-hand side for j=1:(i-1) x(i) = x(i) -A(i,j)*x(j); % subtract prior x values end x(i) = x(i)/A(i,i); % divide by x(i)'s row-i coefficient end INPUT: Test this program with the matrix of Exercise 3 in this section, reverse-indexed so that it is lower rather than upper triangular:: >> a=[3 0 0 0 0; -1 -2 0 0 0; -2 -1 1 0 0; 7 2 6 -2 0; -1 2 2 -1 4] a = 3 0 0 0 0 -1 -2 0 0 0 -2 -1 1 0 0 7 2 6 -2 0 -1 2 2 -1 4 >> b=[6 10 3 0 4]' b = 6 10 3 0 4 OUTPUT: Check agreement with (x1,x2,x3,x4,x5) from Exercise 3: >> forsub(a,b) ans = 2 -6 1 4 5 Algorithm 3 of Section 3.3, p.125. Apply forsub() to a matrix constructed on the command line: >> for i=1:20 b(i,1)=i; for j=1:i a(i,j)=i+j; end end >> b b = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 >> a a = 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 5 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 6 7 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 7 8 9 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 7 8 9 10 11 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 9 10 11 12 13 14 0 0 0 0 0 0 0 0 0 0 0 0 0 9 10 11 12 13 14 15 16 0 0 0 0 0 0 0 0 0 0 0 0 10 11 12 13 14 15 16 17 18 0 0 0 0 0 0 0 0 0 0 0 11 12 13 14 15 16 17 18 19 20 0 0 0 0 0 0 0 0 0 0 12 13 14 15 16 17 18 19 20 21 22 0 0 0 0 0 0 0 0 0 13 14 15 16 17 18 19 20 21 22 23 24 0 0 0 0 0 0 0 0 14 15 16 17 18 19 20 21 22 23 24 25 26 0 0 0 0 0 0 0 15 16 17 18 19 20 21 22 23 24 25 26 27 28 0 0 0 0 0 0 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 0 0 0 0 0 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 0 0 0 0 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 0 0 0 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 0 0 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 0 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 >> forsub(a,b) ans = Columns 1 through 15 0.5000 0.1250 0.0625 0.0391 0.0273 0.0205 0.0161 0.0131 0.0109 0.0093 0.0080 0.0070 0.0062 0.0055 0.0050 Columns 16 through 20 0.0045 0.0041 0.0038 0.0035 0.0032