Numerical analisis

Numerical Methods in Chemical Engineering

by Roberto Lopez Ruiz Chemical Engineering Department Phoenix University – USA

Solution
a) Gauss elimination % GAUSS ELIMINATION CODE % Forward Elimination A=[0.4096 0.1234 0.3678 0.2943 0.4043;0.2246 0.3872 0.4015 0.1129 0.1550;0.3645 0.1920 0.3781 0.0643 0.4240;0.1784 0.4002 0.2786 0.3927 0.2557] n=4; for k=1:n-1; for i=k+1:n; em=A(i,k)/A(k,k);A(i,n+1)=A(i,n+1)-em*A(k,n+1); for j=1:n; A(i,j)=A(i,j)-em*A(k,j); end end end

% Back substitution x(n)=A(n,n+1)/A(n,n); for i=n-1:-1:1; x(i)=A(i,n+1); for j=n:-1:i+1; x(i)=x(i)-A(i,j)*x(j); end x(i)=x(i)/A(i,i); end x

Result
A= 0.4096 0.1234 0.3678 0.2943 0.4043 0.2246 0.3872 0.4015 0.1129 0.1550 0.3645 0.1920 0.3781 0.0643 0.4240 0.1784 0.4002 0.2786 0.3927 0.2557 x= 3.4606 1.5610-2.9342 -0.4301 b) Gauss elimination with pivoting %GAUSS ELIMINATION WITH PIVOTING CODE % Forward Elimination A=[0.4096 0.1234 0.3678 0.2943 0.4043;0.2246 0.3872 0.4015 0.1129 0.1550;0.3645 0.1920 0.3781 0.0643 0.4240;0.1784 0.4002 0.2786 0.3927 0.2557] n=4; for k=1:n-1; % Pivoting piv=abs(A(k,k)); p=k;

for i=k+1:n; if abs(A(i,k))>piv; piv=abs(A(i,k)); p=i; end end if p>k; matrix=A(k,:);A(k,:)=A(p,:); A(p,:)=matrix; end for i=k+1:n; em=A(i,k)/A(k,k); A(i,n+1)=A(i,n+1)-em*A(k,n+1); for j=1:n; A(i,j)=A(i,j)-em*A(k,j); end end end % Back substitution x(n)=A(n,n+1)/A(n,n); for i=n-1:-1:1; x(i)=A(i,n+1); for j=n:-1:i+1; x(i)=x(i)-A(i,j)*x(j); end x(i)=x(i)/A(i,i); end

Result
A= 0.4096 0.1234 0.3678 0.2943 0.4043 0.2246 0.3872 0.4015 0.1129 0.1550 0.3645 0.1920 0.3781 0.0643 0.4240 0.17840.4002 0.2786 0.3927 0.2557

x= 3.4606 1.5610 -2.9342 -0.4301 c) Gauss-Jordan elimination with pivoting %GAUSS_JORDAN ELIMINATION WITH PIVOTING CODE % Forward Elimination A=[0.4096 0.1234 0.3678 0.2943 0.4043;0.2246 0.3872 0.4015 0.1129 0.1550;0.3645 0.1920 0.3781 0.0643 0.4240;0.1784 0.4002 0.2786 0.3927 0.2557] n=4; for k=1:n-1; % Pivoting piv=abs(A(k,k)); p=k; for i=k+1:n; if abs(A(i,k))>piv;piv=abs(A(i,k)); p=i; end end if p>k; matrix=A(k,:); A(k,:)=A(p,:);

A(p,:)=matrix; end A(k,:)=A(k,:)/A(k,k); for i=k+1:n; A(i,:)=A(i,:)-A(i,k)*A(k,:); end end % Back Substitution for k=n:-1:2; A(k,:)=A(k,:)/A(k,k); for i=k-1:-1:1; A(i,:)=A(i,:)-A(i,k)*A(k,:); end for j=1:n x(j)=A(j,n+1); end end x

Result
A= 0.4096 0.1234 0.3678 0.2943 0.4043 0.2246 0.3872 0.4015 0.1129 0.1550 0.3645 0.19200.3781 0.0643 0.4240 0.1784 0.4002 0.2786 0.3927 0.2557

x= 3.4606 1.5610 -2.9342 -0.4301

d) Gauss-Seidel iteration %GAUSS-SEIDEL ITERATION CODE % Reordering matrix A=[0.4096 0.1234 0.3678 0.2943 0.4043;0.3645 0.1920 0.3781 0.0643 0.4240;0.2246 0.3872 0.4015 0.1129 0.1550;0.1784 0.4002 0.2786 0.3927 0.2557] %initial value x=[0; 0; 0; 0]; %Rows number n=4; %iterations number k=0; %NormN=0.000005; while N>0.000001 %convergence criterion k=k+1; for i=1:n R(i)=A(i,n+1); for j=1:n R(i)=R(i)-A(i,j)*x(j); %residual end x(i)=x(i)+R(i)/A(i,i); end N=((R(1)^2)+(R(2)^2)+(R(3)^2)+(R(4)^2))^0.5; end x k

Result
A (reordering)= 0.4096 0.3645 0.2246 0.1784 x= 3.4606 1.5610 -2.9342 -0.4301 k = 32 e) Matrix inversion techniques using Matlab %MATRIX INVERSION TECHNIQUES A=[0.4096 0.1234 0.36780.2943;0.2246 0.3872 0.4015 0.1129;0.3645 0.1920 0.3781 0.0643;0.1784 0.4002 0.2786 0.3927] b=[0.4043 0.1550 0.4240 0.2557] 0.1234 0.1920 0.3872 0.4002 0.3678 0.3781 0.4015 0.2786 0.2943 0.0643 0.1129 0.3927 0.4043 0.4240 0.1550 0.2557

B=inv(A) x=B*b'

Result
A= 0.4096 0.1234 0.3678 0.2943 0.2246 0.3872 0.4015 0.1129 0.3645 0.1920 0.3781 0.0643 0.1784 0.4002 0.2786 0.3927 B= -8.1360 -12.865016.3274 7.1226 -8.1860 -6.2722 9.9782 6.3042

11.3186 15.3595 -17.2531 -10.0733 4.0085 x= 3.4606 1.5610 -2.9342 -0.4301 1.3397 -5.3460 0.0325

Obtain the solutions to 4 decimal places. Compare the solutions by tabulating the results. x1
Gauss Elimination(G-E) G-E with pivoting Gauss-Jordan with pivot. Gauss-Seidel iteration Matrix inversion (Matlab)

x2 1.5610 1.5610 1.5610 1.5610 1.5610...