MATLAB Code Help. Using Crout's Method, solve the system of linear equation $Mz=f$, where $M=\begin{pmatrix}I &A\\A^T&0\end{pmatrix}$

343 Views Asked by At

Using Crout's Method, solve the system of linear equation $Mz=f$, where $$M=\begin{pmatrix}I &A\\A^T&0\end{pmatrix}$$

I have implemented algorithm of Crout's method. But I don't have any idea how to create this $M$ function in MATLAB, and implement it in my algorithm. It is necessary to solve $Mz=f$ where $M$ is block-matrices. I need help with implementing this $M$ matrix. Please help me.

Code of Algorithm

function [L,U] = Crouts(A,B,n)

%n=input('Enter the size of the matrix\n');

%%%%%%%%%%EXAMPLES%%%%%%%%%%

%Please uncomment an example you prefer to use, or type it on command window

%example 0§
%A = [2 -6 6; 3 -7 13; -2 2 -11]
%B = [2; -13; 21]

%example 1
%A=[5 4 1; 10 9 4; 10 13 15]
%B=[3.4;8.8;19.2]

%example 2 
%A = [2 1 4; 8 -3 2; 4 11 -1]
%B = [12 ;20 ;33]

%example 3
%A = [4 -2 -3 6; -6 7 6.5 -6; 1 7.5 6.25 5.5; -12 22 15.5 -1]
%B = [ 12; -6.5; 16; 17 ]

%example 4
%A = [9 -4 -2 0; -4 17 -6 -3; -2 -6 14 -6; 0 -3 -6 11]
%B = [24; -16; 0; 18]

%example 4
%A=rand(n) , B= rand(n,1)

%example 5 (pdf)
%A = [10 3 4; 2 -10 3; 3 2 -10]
%B = [15 37 -10]

%example 6 (pdf)
%A = [9 3 3 3; 3 10 -2 -2; 3 -2 18 10; 3 -2 10 10]
%B = [24; 17; 45; 29]


%We should remeber that the n in both A and B has to be the same number so the method works

% Here we initialize L and U in order to prevent dynamic allocation during the process

%if one of the main diagonal's entries is 0, program doesnt calculate LU
for i=1:n
   if A(i,i) == 0
       L = 0;
       U=0;           
       det =0;     %if so then return 0's
       x=0;
       return;
   end
end

% initializing L and U,to prevent dynamic allocation during the process
L = zeros(n,n);
U = zeros(n,n);

%Substituting 1's im the diagonal of U
for i=1:n
   U(i,i)=1;
end

%Calculating the 1st column of L
%the first column of L in Crout's factorization is always equal to the
%first column of A

for i=1:n
   L(i,1) = A(i,1);
end


%Calculating the elements in the first row of U(Except U11 which already
%was already calculated
for i=2:n
   U(1,i) = A(1,i) / L(1,1);
end

%calculating the remaining  elements row after row.The elements of  L are calculated first
%because they are used for calculating the elements of U
for i = 2:n
   for j = 2:i
       L(i, j) = A(i, j) - L(i, 1:j - 1) * U(1:j - 1, j);%%%formula
   end
           
   for j = i + 1:n
       U(i, j) = (A(i, j) - L(i, 1:i - 1) * U(1:i - 1, j)) / L(i, i);%%%formula
   end  
end
L
U

%det(A)=det(L)*det(U). 
% As we know, det(U)=1.
%Calculate det(L), which is
%the product of diagonal elements
detA=1;
for i=1:n
   detA=detA*L(i,i);
end

detA

% AX=B
% let A=LU => LUX=B
% let UX=Y => LY=B

%forward substitution

Y = zeros(n,1); %Intializing Y
Y(1) = B(1)/L(1,1);
for m=2:n
    Y(m) = (B(m) - L(m,1:m-1)*Y(1:m-1))/L(m,m);
end
Y

%backward substitution 

X = zeros(n,1);
X(n) = Y(n)/U(n,n);
for m=n-1:-1:1
   X(m) = (Y(m) - U(m,m+1:n)*X(m+1:n))/U(m,m);
   
end
X

end