jacobi method, gauss siedel for solving linear equations

EC 313: Numerical Methods Mayank Awasthi(2006033) MATLAB Assignment – 2 B.Tech 3rd Year(ECE) Pdpm IIITDM, JABALPUR Ques1: Gauss Siedel Method for Solving the Linear Equations: % Implementing Gauss Siedel method % a is matrix whose diagonal elements are non-zero else rearrange it function x = gauss_siedel(a,b,x) n=length(a); l=zeros(n,n); u=zeros(n,n); d=zeros(n,n); id=[1,0,0;0,1,0;0,0,1]; for i=1:n for j=1:n a(i,j)= a(i,j)/a(i,i); %making the diagonal elements 1. % Breaking the matrix as a sum of ( L+U+I ) if i>j l(i,j)=a(i,j); end if i<j u(i,j)=a(i,j); end if i==j d(i,j)=a(i,j); end end end %Implementing Norm of c where c = inverse(I+L)*U norm2(inv2(id+l)*u); if norm2(inv2(id+l)*u)>1 fprintf('Norm of c is greater than 1, Solution will diverge'); return; end for i=1:n b(i,1)=b(i,1)/a(i,i); end % Calculating x using FORWARD DIFFERENCE CONCEPT x=[0;0;0]; for i=1:100 x= forward_subs((l+d),(b-u*x)); m=a*x-b; % Calculating(dis(A*x,b)<1e-6)to stop iteration at the given tolerance temp=0; for j=1:n temp=temp+power(m(j,1),2); end temp=sqrt(temp); if temp<0.0000001 break; end end fprintf('\nThe maximum no. of iteration required is %d',i); end MATLAB routine for Inverse of 3x3 matrix : % Calculating Inverse of 3x3 matrix function x = inv2(A) I=[1,0,0;0,1,0;0,0,1]; n=length(A); a21=A(2,1)/A(1,1); for k = 1:n-1 % Elimination Phase for i= k+1:n if A(i,k) ~= 0 lambda = A(i,k)/A(k,k); A(i,k:n) = A(i,k:n) - lambda*A(k,k:n); I(i,k:n)= I(i,k:n) - lambda*I(k,k:n); I(3,1)=I(2,1)*I(3,2)+I(3,1); end end end x1=back_subs(A,I(:,1)); % Backward Substitution x2=back_subs(A,I(:,2)); x3=back_subs(A,I(:,3)); x=[x1,x2,x3]; end MATLAB routine for Norm: % Calculating Norm of matrix function n0= norm2(c) l=length(c); n0=0; for m=1:l for n=1:l n0 = n0 + c(m,n)*c(m,n); end end n0=sqrt(n0); end MATLAB routine for Forward Substitution: % Implementing Forward Substitution function x=forward_subs(A,b) n=length(A); for i = 1:3 t=0; for j = 1:(i-1) t=t+A(i,j)*x(j); end x(i,1)=(b(i,1)-t)/A(i,i); end Part1: Part2: >> a=[1 ,.2 .2;.2, 1, .2;.2, .2, 1] >> a=[1 ,.5 .5;.5, 1, .5;.5, .5, 1] a = a = 1.0000 0.2000 0.2000 1.0000 0.5000 0.5000 0.2000 1.0000 0.2000 0.5000 1.0000 0.5000 0.2000 0.2000 1.0000 0.5000 0.5000 1.0000 >> b=[2;2;2] >> b=[2;2;2] b = b = 2 2 2 2 2 2 >> x=[0;0;0] >> x=[0;0;0] x = x = 0 0 0 0 0 0 >> gauss_siedel(a,b,x) >> gauss_siedel(a,b,x) The maximum no. of iteration required is 8 The maximum no of iteration required is 17 ans = ans = 1.4286 1.0000 1.4286 1.0000 1.4286 1.0000 Part3: >> a=[1 ,.9 .9;.9, 1, .9;.9, .9, 1] a = 1.0000 0.9000 0.9000 0.9000 1.0000 0.9000 0.9000 0.9000 1.0000 >> b=[2;2;2] b = 2 2 2 >> x=[0;0;0] x = 0 0 0 >> gauss_siedel(a,b,x) Norm of c is greater than 1, Solution will diverge ans = 0 0 0 ************************************************************************ Spectral radius of C: Spectral radius of C is maximum eigenvalue of C*CT. In this case C=(I+L)-1 * U 1). >> a=[1 ,.2 .2;.2, 1, .2;.2, .2, 1] a = 1.0000 0.2000 0.2000 0.2000 1.0000 0.2000 0.2000 0.2000 1.0000 >> L=[0,0,0;.2,0,0;.2,.2,0] L = 0 0 0 0.2000 0 0 0.2000 0.2000 0 >> I=[1,0,0;0,1,0;0,0,1] I = 1 0 0 0 1 0 0 0 1 >> U=[0,.2,.2;0,0,.2;0,0,0] U = 0 0.2000 0.2000 0 0 0.2000 0 0 0 >> C=inv2(I+L)*U C = 0 0.2000 0.2000 0 -0.0400 0.1600 0 -0.0320 -0.0720 >> lamda=eig(C*C') lamda = 0.0000 0.0181 0.0953 >> SR=max(lamda) SR = Spectral Radius 0.0953 2). >> a=[1 ,.5 .5;.5, 1, .5;.5, .5, 1] a = 1.0000 0.5000 0.5000 0.5000 1.0000 0.5000 0.5000 0.5000 1.0000 >> L=[0,0,0;.5,0,0;.5,.5,0] L = 0 0 0 0.5000 0 0 0.5000 0.5000 0 >> I=[1,0,0;0,1,0;0,0,1] I = 1 0 0 0 1 0 0 0 1 >> U=[0,.5,.5;0,0,.5;0,0,0] U = 0 0.5000 0.5000 0 0 0.5000 0 0 0 >> C=inv2(I+L)*U C = 0 0.5000 0.5000 0 -0.2500 0.2500 0 -0.1250 -0.3750 >> lamda=eig(C*C') lamda = 0.0000 0.1481 0.6332 >> SR=max(lamda) SR = Spectral Radius 0.6332 3). >> a=[1 ,.9 .9;.9, 1, .9;.9, .9, 1] a = 1.0000 0.9000 0.9000 0.9000 1.0000 0.9000 0.9000 0.9000 1.0000 >> L=[0,0,0;.9,0,0;.9,.9,0] L = 0 0 0 0.9000 0 0 0.9000 0.9000 0 >> I=[1,0,0;0,1,0;0,0,1] I = 1 0 0 0 1 0 0 0 1 >> U=[0,.9,.9;0,0,.9;0,0,0] U = 0 0.9000 0.9000 0 0 0.9000 0 0 0 >> C=inv2(I+L)*U C = 0 0.9000 0.9000 0 -0.8100 0.0900 0 -0.0810 -0.8910 >> lamda=eig(C*C') lamda = 0.0000 0.7301 2.3546 >> SR=max(lamda) SR = Spectral Radius 2.3546 Steps vs t : Spectral radius vs t : Ques2: Jacobi Method for Solving the Linear Equations: %Implementing Jacobi Method function x = jacobi(A,b,x) I=[1 0 0; 0 1 0; 0 0 1]; c=I-A; %Checking Norm if norm2(c) > 1 fprintf('\nNorm is greater than one so it will diverge'); return; end for j=1:50 x=b+(I-A)*x; n=A*x-b; %checking the condition of tolerance to stop the iterations temp=0; for i = 1:3 temp= temp+ power(n(i,1),2); end temp=sqrt(temp); if temp < 0.0001 break; end end fprintf('\nThe iteration required is %d',j); end MATLAB routine for Norm: % Calculating Norm of matrix function n0= norm2(c) l=length(c); n0=0; for m=1:l for n=1:l n0 = n0 + c(m,n)*c(m,n); end end n0=sqrt(n0); end Part1: Part2 >> a=[1 ,.2 .2;.2, 1, .2;.2, .2, 1] >> a=[1 ,.5 .5;.5, 1, .5;.5, .5, 1] a = a = 1.0000 0.2000 0.2000 1.0000 0.5000 0.5000 0.2000 1.0000 0.2000 0.5000 1.0000 0.5000 0.2000 0.2000 1.0000 0.5000 0.5000 1.0000 >> b=[2;2;2] >> b=[2;2;2] b = b = 2 2 2 2 2 2 >> x=[0;0;0] >> x=[0;0;0] x = x = 0 0 0 0 0 0 >> jacobi(a,b,x) >> jacobi(a,b,x) The iteration required is 12 Norm is greater than one so it will diverge ans = ans = 1.4285 0 1.4285 0 1.4285 0