Please write clearly, typing is also ok. use matlab 2.11 (a) Write programs impl
ID: 2267652 • Letter: P
Question
Please write clearly, typing is also ok.
use matlab
2.11 (a) Write programs implementing Gaus- sian elimination with no pivoting, partial piv- oting, and complete pivoting. (b) Generate several linear systems with ran- dom matrices (i.e., use a random number gen- craior to olialll the airixculines) and riplii. hand sides chosen so that the solutions are known, and compare the accuracy, residuals, and performance of the three implementations. c) Can you devise a (nonrandom) matrix for which complete pivoting is significantly more te than partial pivoting?Explanation / Answer
ANS) a)
Gauss elimination with complete pivoting.
G)aussian (E)limination (C)omplete (P)ivoting
Input
A nxn matrix
Output
L = Lower triangular matrix with ones as diagonals
U = Upper triangular matrix
P and Q permutations matrices so that P*A*Q = L*U
examples :
[L U] = gecp(A);
[L U P] = gecp(A);
[L U P Q] = gecp(A);
MATLAB PROGRAM
function [L, U, P, Q] = gecp(A)
%GECP calculate Gauss elimination with complete pivoting
%
% (G)aussian (E)limination (C)omplete (P)ivoting
% Input : A nxn matrix
% Output
% L = Lower triangular matrix with ones as diagonals
% U = Upper triangular matrix
% P and Q permutations matrices so that P*A*Q = L*U
%
% See also LU
%
% written by : Cheilakos Nick
[n, n] = size(A);
p = 1:n;
q = 1:n;
for k = 1:n-1
[maxc, rowindices] = max( abs(A(k:n, k:n)) );
[maxm, colindex] = max(maxc);
row = rowindices(colindex)+k-1; col = colindex+k-1;
A( [k, row], : ) = A( [row, k], : );
A( :, [k, col] ) = A( :, [col, k] );
p( [k, row] ) = p( [row, k] ); q( [k, col] ) = q( [col, k] );
if A(k,k) == 0
break
end
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
i = k+1:n;
A(i,i) = A(i,i) - A(i,k) * A(k,i);
end
L = tril(A,-1) + eye(n);
U = triu(A);
P = eye(n);
P = P(p,:);
Q = eye(n);
Q = Q(:,q);
Gaussian elimination with partial pivoting, which produces a factorization of P*A into the product L*U where P is a permutation matrix, and L and U are lower and upper triangular, respectively.
%
% Example:
% x = gee_its_simple (A,b) ; % x = A using Gaussian elimination
% x = gee_its_short (A,b) ; % same as gee_its_simple, just shorter
%
% For production use:
% Use x=A instead of x = gee_its_simple (A,b)
% Use x=A instead of x = gee_its_short (A,b)
% Use x=L instead of x = gee_its_simple_forwardsolve (L,b)
% Use x=U instead of x = gee_its_simple_backsolve (U,b)
% Use [L,U,p]=lu(A) and rcond(A) for [LU,p,rcnd] = gee_its_simple_factorize(A)
%
% Primary Files:
% gee_its_short - x=A, no comments or error checking (just for line count)
% gee_its_simple - solves A*x=b using a Gaussian Elimination Example (it's simple!)
%
IN ANOTHER METHOD FOR PARTIAL PIVOTING
%create the augmented matrix A|B
Aug=[A B];
n=rank(A);
%initialize the nrow vector
for i=1:n
nrow(i)=i;
end
nrow=nrow';
for k=1:n-1;
max=0;
index=0;
%find the maximum value in the column under the current checked element and
%return its row position
for j=k:n
if abs(Aug(nrow(j),k))>max
max=abs(Aug(nrow(j),k));
index=j;
end
end
%perform row exchange in the nrow vector
if nrow(k)~=nrow(index)
ncopy=nrow(k);
nrow(k)=nrow(index);
nrow(index)=ncopy;
disp(sprintf ('row changed '))
else
disp(sprintf ('no change '))
end
%Gaussian elimination
for i=(k+1):n
m(nrow(i),k)=Aug(nrow(i),k)/Aug(nrow(k),k);
for j=k:n+1
Aug(nrow(i),j)=Aug(nrow(i),j)-m(nrow(i),k)*Aug(nrow(k),j);
end
end
end
%backward subsitution
x(n)=0;
x=x';
x(n)=Aug(nrow(n),n+1)/Aug(nrow(n),n);
i=n-1;
while i>0
x(i)=(Aug(nrow(i),n+1)-Aug(nrow(i),i+1:n)*x(i+1:n))/(Aug(nrow(i),i));
i=i-1;
end
x_soln=x;
A_aug=Aug;
Gaussian elimination without Pivoting
function x = Gauss(A, b)
% Solve linear system Ax = b
% using Gaussian elimination without pivoting
% A is an n by n matrix
% b is an n by k matrix (k copies of n-vectors)
% x is an n by k matrix (k copies of solution vectors)
[n, n] = size(A); % Find size of matrix A
[n, k] = size(b); % Find size of matrix b
x = zeros(n,k); % Initialize x
for i = 1:n-1
m = -A(i+1:n,i)/A(i,i); % multipliers
A(i+1:n,:) = A(i+1:n,:) + m*A(i,:);
b(i+1:n,:) = b(i+1:n,:) + m*b(i,:);
end;
% Use back substitution to find unknowns
x(n,:) = b(n,:)/A(n,n);
for i = n-1:-1:1
x(i,:) = (b(i,:) - A(i,i+1:n)*x(i+1:n,:))/A(i,i);
end
% Secondary Files:
% gee_its_simple_factorize - Gaussian Elimination Example, with partial pivoting
% gee_its_simple_forwardsolve - computes x=L where L is unit lower triangular
% gee_its_simple_backsolve - x=U where U is upper triangular
% gee_its_simple_resid - compute the relative residual for x=A
% gee_its_simple_test - tests the "Gee! It's Simple!" package
%
% Just for fun:
% gee_its_sweet - solves Ax=b with just x=A; it doesn't get sweeter than this
% gee_its_too_short - x=A, no pivoting (thus unstable!), just bare bones