HW3 sol - HW # 3 Problem 1: Page 139—Problem 15 (a) By...

Info iconThis preview shows pages 1–3. Sign up to view the full content.

View Full Document Right Arrow Icon
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
Background image of page 3
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: HW # 3 Problem 1: Page 139—Problem 15 (a) By using exact values for elements of A and applying Gaussian elimination and then back substitution one can find the exact solution to AX=B , that is 1 1/2 1/3 1/4]1| |1/2 1/3 1/4 1/510] |1/3 1/4 1/5 1/6|0| Eq(l) |1/4 1/5 1/6 1/7|0| However, not to work by hand and also not to mess with fraction, that can result in loss of sig, let' mult each row by the common denominator of row elements of A. i.e. lst row by 1*2*3*4, 2nd row by 2*3*4*5, etc doing so changes Eq(l) to: 24 12 8 6 24 a: 60 40 30 24 0 120 90 72 60 0 Eq(2) 210 168 140 120 0 Do Gaussian elimination of first column by mult Eq(2) by 1 0 0 0 —60 24 0 0 = —120 0 24 G —210 0 0 24 to get 24 12 8 6 24 0 240 240 216 —1440 b=p*a= 0 720 768 720 —2880 Eq(3) 0 1512 1680 1620 "5040 mult b by 1 0 0 0 0 1 0 0 pl: 0 —720 240 0 0 —1512 0 240 to get ‘ 24 12 8 6 24 0 240 240 216 —1440 c=p1*b 0 0 11520 17280 345600 0 0 40320 62208 967680 Since 3d row is divisible by 10 and 4th row by 8, let's reduce them to 24 12 8 6 24 0 240 240 216 —1440 d= 0 0 1152 1728 34560 0 0 5040 7776 120960 mult d by p2 to make A to become upper triangular 1 0 0 G 0 1 0 0 p2 0 0 1 0 0 0 —5040 1152 24 12 8 6 24 0 240 240 216 —1440 e=p2*d 0 0 1152 1728 34568 0 0 0 248832 —34836480 Now find exact x4, x3, x2 and x1 by back substitution (see program) i.e x4: ~34836480/248832 =—140, x3=240, x2=~120 and x1=16 (b) If one inputs elements of A with only 4 decimals, one finds the solution to AX=B to be: X=[18.7308 —149.6053 310.0628 —185.0881]' Such a large change in solution (Delta X) = [2.7308 29.6053 70.0628 45.0881]' with only small change in elements of A indicates that Hilbert matrix A is an "ill—conditioned” matrix. See Problem 3, for calculating its condition #. Probelm 2: Page 165—problem 6 The system as written, its A matrix is not strictly diagonally dominant. Thus, convergence is not guaranteed. Let't test the method using the ‘ matrix A as written and find the result for 3 iterations. Jacobi x y 2 Norm (AX—B)/10“3 k=1 5.5000 ~10.0000 0.7500 0.0808 k=2 45.8750 18.2500 4.6250 0.2221 k=3 ~65.1875 224.0000 7.6562 1.6430 Gauss~Seidel x y 2 Norm (AX—B)/10“4 k=1 5.5000 17.5000 —2.2500 0.0142 k=2 —65.6250 —340.3750 69.4375 0.2935 k=3 (1.4017 7.0680 —1.4158)*10“3 6.0753 Jacobi and Gauss—Seidel seem not to converge. Now let's write A in a. strictly diagonal form; interchange the first equation with the second equation of the system and run the same program again, to find Jacobi x y 2 Norm (AX—B) k=1 2.0000 1.3750 0.7500 3.2500 k=2 2.1250 0.9688 0.9062 0.5625 k=3 2.0125 0.9570 1.0391 0.5625 Gaussteidel x y 2 Norm (AX-B) k=1 2.0000 0.8750 1.0312 1.0312 k=2 1.9688 1.0117 0.9893 0.1787 k=3 2.0045 0 9975 1.0017 0.0267 Now that A is written in a strictly diagonally dominant form, both Jacobi and Gauss~Seidle iterations seems to be converging, with Gauss—Seidel doing a little better as expected. In fact if we do 31 iterations, we get the following Jacobi result (x,y,z)=(2.000000000000024, 1.000000000000001, 0.999999999999975) with error (AX— B): 1.203481758693670e—13, using uniform norm. However, with onl 18 iteration, we get comparable results using Gauss—Seidel method. (x, y, z)=(1.999999999999972, 1.000000000000015, 0.999999999999989) with error (AX— B): 1.669775429036235e—13, using uniform norm. Clearly both methods are converging to exact the answer x=2, y=1 and 2:1 Problem 3. By calculating |Delta A|/|A|=3.6571e—05 and mu=|A|]A|A{—1}=2.8375e+04 using the uniform norm (see matlab program below), one notices the condition number, mu is too large for the size of A, thus it should not be of any'surprise to see such a large change in the answer, i.e. Delta x3: 70.0628. The relative error |Delta X[/|X] = 0.2919, for a very small change in A. This could be understood, if we calculate expected order of relative error in X which is of the order mu*|Delta Al/|A|=2.8375e+04*3.6571e—05=1.0377. The exact formula for the bound on relative error in X if |Delta A] < 1/[A“{—1}| is |Delta X[/|X|< (mu/(1wmu|Delta A|/|A|))(|Delta A|/|A| + |Delta Bl/IBI). However for this example the condition |Delta A| < 1/1A“{—1}| is not satisfied. Thus, the above inequality can't be used. However, still the order of relative error being < (mu*relative error in A) is valid and consistent with our findings that 0.2919 < 1.0377. ——————————————— ——Matlab programs for HW3 *-———————www~—————————— % part (a) problem 1. a=[2*3*4 1*3*4 1*2*4 1*2*3 2*3*4; 3*4*5 2*4*5 2*3*5 2*3*4 0; 4*5*6 3*5*6 3*4*6 3*4*5 0 ; 5*6*7 4*6*7 4*5*7 4*5*6 0] B=[*1 0 0 0; —a(2,l) a(1,1) 0 0; —a(3,1) 0 a(1,1) 0; —a(4,1) 0 0 a(l,l)]; :p a; p1=[1 0 0 0; 0 l 0 0; 0 —b(3,2) b(2,2) 0; 0 ~b(4,2) 0 b(2,2)]; c=p1*b; q= [ 1 0 0 0; 0 1 0 0; 0 0 1/10 0; 0 0 0 1/8]; d=q*c; p2=[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 —d(4,3) d(3,3)]; e=p2*d; x4=e(4,5)/e(4,4); x3=(e(3,5)—e(3,4)*x4)/e(3,3); 'x2=(e(2,5)—e(2,4)*x4—e(2,3)*x3)/e(2,2); x1=(e(1,5)—e(1,4)*x4—e(1,3)*x3—e(1,2)*x2)/e(1,1); % part (b) problem 1. a1=1; a2=1/2; a3=1/3; a4=1/4; a5=1/5; a6=1/6; a7=1/7; A=[a1 a2 a3 a4;a2 a3 a4 a5; a3 a4 a5 a6; a4 a5 a6 a7]; B=[1 0 0 01'; IA=inV(A); X=[x1 x2 x3 x4]'; A1=[1.0000 0.5000 0.3333 0.2500 0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429]; IA1=inv(A1); X1=IA1*B; dX=abs(X1—X); %Problem 2 %We run it once for the A' and B's as they appear in the problem, which is %not a diagonally dominant matirix. We then interchange row one and row %two of the system, to make the associated A diagonally dominant and use %the Gauss—Seidel iteration method to see if it is converging. a=[ 2 8 —1; 5 —1 l; —1 1 4]; %b=[11 10 3]'; a:[5 —1 1; 2 8 “1; —1 l 4]; b=[10 11 3]'; x0=0: y0=0; 20:0; max1=31 Xj(1)= (b(1)-a(1.2)*y0 -a(1.3)*20)/a(1.1): yj(1)= (b(2)—a(2.1)*x0 —a(2.3)*20)/a 2.2); zj(1)= (b(3)—a(3.1)*xo -a(3.2)*y0)/a(3.3): x(1)= (b(1)-a(1,2)*y0 -a(1.3)*20)/a(1.1): y(1)= (b(2)—a(2.l)*X(1) —a(2.3)*20)/a(2.2); 2(1): (b(3)—a(3.1)*X(1) ma(3.2)*y(l))/a(3.3). for k=2:maxl xj(k)= (b(l)-a(1.2)*yj(kel) —a(l.3)*21(kml))/a(1,l): yJ(k)= (b(2)-a(2.1)*x1(k-l) -a(2.3)*ZJ(k-1))/a(2.2): 2j(k)= (b(3)"a(3.1)*X3(k“1) -a(3,2)*yj(k-1))/a(3.3); X(k)= (b(l)—a(l.2)*y(k-1) —a(l.3)*z(k—l))/a(1.1); y(k)= (b(2)-a(2.1)*x(k) —a(2.3)*2(k—1))/a(2.2); z(k)= (b(3)—a(3.1)*X(k) -a(3,2)*y(k))/a(3.3): end; format long for k=1:max1 Xj=[Xj(k) yj(k) zj(k)} X=[X(k) Y(k) z(k)] Cj=a*Xj' —b; Cjn(k)=norm(Cj,inf); C=a*X' —b: Cn(k)=norm(C,inf); end; % Problem 3 aX=abs(X); [ndX il]=max(dX); [nX i2]=max(aX); dA=abs(AmA1); aA=abs(A); aIA=abs(IA); for n=1:4 sumd=0; sumA=0; sumI=0; for m=1:4 sumd=sumd+dA(n,m); sumA=sumA+aA(n,m); sumI=sumI+aIA(n,m); end; nmdA(n)=sumd; nmaA(n)=sumA; nmaIA(n)=sumI; end; [ndA i3]=max(nmdA); [naA i4]=max(nmaA); [naIA 15]=max(nmaIA); mu=naA*naIA: RnA=ndA/naA; RX=ndX/nX; ...
View Full Document

Page1 / 3

HW3 sol - HW # 3 Problem 1: Page 139—Problem 15 (a) By...

This preview shows document pages 1 - 3. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online