Modeling of a special physical phenomenon results into the f

Modeling of a special physical phenomenon results into the following linear system of equations: x_i-1 - 2x_i + x_i + 1 = i times P in which P = 12. Suppose that we want to solve a 5 times 5 system of equations. In this problem, x_0 is always equal to zero x_0 = 0 and x_6 = 22. Note that x_0 and x_6 are known and x_1 to x_5 are unknown. Make the system of equations and write it down in matrix format, [A][x] = [b]. Write a MATLAB program to Determine the condition number of your system using row-sum norm, column-sum norm, and Frobenius norm approaches. Using MATLAB built-in function (Cond) to find condition number is NOT allowed. Solve the system using Naive Gauss Elimination method (NaiveGaussElim.m), Gauss Elimination method with partial pivoting and without scaling (GaussEIm.m). Gauss Jordan method with partial pivoting and without scaling (Gaussjord m), and LU decomposition method without pivotmg (Ludecomp m). Validate your result by comparing them to the one obtained by using MATLAB built-in command for solving a linear system of equations (x = A\\b) Obtain the inverse of the matrix [A]^-1 using Gauss Jordan method (no pivoting). Obtain the inverse of the matrix [A]^-1 using LU decomposition method (no pivoting). Check your solutions in part (c) and (d) by confirming [A][A]^-1 = 1.Note: Write one main M-file and call other functions (e.g., NaiveGaussEIim, GoussElim GaussJord Ludecomp...)

Solution

Part has been embedded in main matlab file.

Save the files as separate M files in same folder to run

MAIN FILE

clear all
clc
%%
A=[-2,1,0,0,0;1,-2,1,0,0;0,1,-2,1,0;0,0,1,-2,1;0,0,0,1,-2];
b=[12;24;36;48;38];

B=inv(A);

%% a.

% row norm
s1=max(sum(abs(A),2)); % sum along rows of A
s2=max(sum(abs(B),2)); % sum along rows of inv(A)
condition_number_row_sum_norm=s1*s2
  
% coloumn sum norm
s1=max(sum(abs(A),1)); % sum along rows of A
s2=max(sum(abs(B),1)); % sum along rows of inv(A)
condition_number_column_sum_norm=s1*s2

% Frobenius norm
s1=sqrt(sumsqr(A)); % sum along rows of A
s2=sqrt(sumsqr(B)); % sum along rows of inv(A)
condition_number_Frobenius_sum_norm=s1*s2

  
%% b.
X1=NaiveGaussElim(A,b);
X2=GaussElim(A,b);
X3=Ludecomp(A,b);
X4=A\\b;
%% c.
X5=inverse_gauss(A)*b;
  
%% d.
X6=LU_inv(A)*b;
  
  
  
  
  
function x = NaiveGaussElim(A,b);
n = length(b); x = zeros(n,1);
for k=1:n-1 % forward elimination
for i=k+1:n
xmult = A(i,k)/A(k,k);
for j=k+1:n
A(i,j) = A(i,j)-xmult*A(k,j);
end
b(i) = b(i)-xmult*b(k);
end
end
% back substitution
x(n) = b(n)/A(n,n);
for i=n-1:-1:1
sum = b(i);
for j=i+1:n
sum = sum-A(i,j)*x(j);
end
x(i) = sum/A(i,i);
end
  
  
  
  
  
  
  
  
function x = GaussElim(A,b)

%GaussPP(A,b) solves the n-by-n linear system of equations using partial
%pivoting
%A is the coeficient matrix
%b the right-hand column vector

n = size(A,1); %getting n
A = [A,b]; %produces the augmented matrix

%elimination process starts
for i = 1:n-1
p = i;
%comparison to select the pivot
for j = i+1:n
if abs(A(j,i)) > abs(A(i,i))
U = A(i,:);
A(i,:) = A(j,:);
A(j,:) = U;
end
end
%cheking for nullity of the pivots
while A(p,i)== 0 && p <= n
p = p+1;
end
if p == n+1
disp(\'No unique solution\');
break
else
if p ~= i
T = A(i,:);
A(i,:) = A(p,:);
A(p,:) = T;
end
end
  
for j = i+1:n
m = A(j,i)/A(i,i);
for k = i+1:n+1
A(j,k) = A(j,k) - m*A(i,k);
end
end
end

%checking for nonzero of last entry
if A(n,n) == 0
disp(\'There is no unique solution\');
return
end

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

function x =Ludecomp(A,b)
% LU factorization of an n by n matrix A
n=size(A);
[L,U,P]=LU_factor(A);

y=L\\b;
x=U\\y;
end

function [L,U,P] = LU_factor(A)
% LU factorization.
%
%   
% [L,U,P] = Lu(A) returns unit lower triangular matrix L, upper
% triangular matrix U, and permutation matrix P so that P*A = L*U.
%
s=length(A);
U=A;
L=zeros(s,s);
PV=(0:s-1)\';
for j=1:s,
% Pivot Voting (Max value in this column first)
[~,ind]=max(abs(U(j:s,j)));
ind=ind+(j-1);
t=PV(j); PV(j)=PV(ind); PV(ind)=t;
t=L(j,1:j-1); L(j,1:j-1)=L(ind,1:j-1); L(ind,1:j-1)=t;
t=U(j,j:end); U(j,j:end)=U(ind,j:end); U(ind,j:end)=t;

% LU
L(j,j)=1;
for i=(1+j):size(U,1)
c= U(i,j)/U(j,j);
U(i,j:s)=U(i,j:s)-U(j,j:s)*c;
L(i,j)=c;
end
end
P=zeros(s,s);
P(PV(:)*s+(1:s)\')=1;

function [Inverse] = inverse_gauss(A)

[m,n] = size(A);
if (m ~= n)
display(\'Matrix must be square for this program\')
return
else
end
x = 1; % iterator for elimination matrix 3rd dimension
for j = 1:1:n-1,
for i = j+1:1:m,
E(:,:,x) = eye(n);
E(i,j,x) = -A(i,j) / A(j,j);
A = E(:,:,x) * A;
x = x + 1;
end
end
x = 1;
for j = n:-1:2,
for i = j-1:-1:1,
P(:,:,x) = eye(n);
P(i,j,x) = -A(i,j) / A(j,j);
A = P(:,:,x) * A;
x = x + 1;
end
end
for i = 1:1:n,
S(:,:,i) = eye(n);
S(i,i,i) = 1 / A(i,i);
A = S(:,:,i) * A;
end

ProdE = 1;
[~,~,c] = size(E);
for ii = c:-1:1,
ProdE = ProdE * E(:,:,ii);
end
ProdP = 1;
[~,~,c] = size(P);
for ii = c:-1:1,
ProdP = ProdP * P(:,:,ii);
end
ProdS = 1;
[~,~,d] = size(S);
for jj = d:-1:1,
ProdS = ProdS * S(:,:,jj);
end
Inverse = ProdS*ProdP*ProdE;
end

function Ainv = LU_inv(A)
% LU factorization of an n by n matrix A

[L,U,P] = LU_factor(A);
% Solve linear system for Identity matrix
I=eye(size(A));
s=size(A,1);
Ainv=zeros(size(A));
for i=1:s
b=I(:,i);
Ainv(:,i)=TriangleBackwardSub(U,TriangleForwardSub(L,P*b));
end

function C=TriangleForwardSub(L,b)
% Triangle Matrix Forward Substitution

s=length(b);
C=zeros(s,1);
C(1)=b(1)/L(1,1);
for j=2:s
C(j)=(b(j) -sum(L(j,1:j-1)\'.*C(1:j-1)))/L(j,j);
end

function C=TriangleBackwardSub(U,b)
% Triangle Matrix Backward Substitution

s=length(b);
C=zeros(s,1);
C(s)=b(s)/U(s,s);
for j=(s-1):-1:1
C(j)=(b(j) -sum(U(j,j+1:end)\'.*C(j+1:end)))/U(j,j);
end