Matlab. LU decomposition Crout's Method. And usage of decomposition on accurate Block Matrix.

330 Views Asked by At

I have to implement such a program()Look at picture I attached

I mostly implemented everything: Crout's Algorithm, solving linear equations, I created this block matrix, but I don't know to use that all to solve the system of linear equations Mz = f INSIDE my program. Help me please.

I attached the code below:

% Explanation of the solution:

%This program decompose the square matrix of any size into a 
% product of a Lower-triangular matrix (L)
% and an Upper-triangular matrix (U). 


 
 function [L,U] = Crouts(A,B,n)
 
 
%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]
%n=3

%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]
% n = 4

%%%POSSIBLE ERRORS

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

%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

%If user input wrong size of matrix(n), program will return 0's
for i=1:n
    if n ~= size(A);
        return;
    end
end

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

%Starting finding L and U
%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
%DETERMINANT 
%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


%TRANSPOSE MATRIX
[m,n] = size(A); % m and n are the rows of matrix x
tran = zeros(n,m); % when we transpose, rows and columns are exchanged
for i= 1:n;
    for j=1:m ;
        tran(i,j) = A(j,i);
    end
end

 tran
 %%% IDENTITY MATRIX
 Identity_Matrix=zeros(n,n);
 for i=1:n
    Identity_Matrix(i,i)=1;
 end
 

E1 = Identity_Matrix;
E2 = A;
E3 = tran;
E4 = zeros(n);
M=[E1,E2;E3,E4]
f = rand(n*2,1)


 end