Crout 방법‐> LU분해‐> 해구하기

function [L, U] = LUdecompCrout(A)

% The function decomposes the matrix A into a lower triangular matrix L

% and an upper triangular matrix U, using Crout's method such that A=LU.

% Input variables:

% A The matrix of coefficients.

% Output variable:

% L Lower triangular matrix.

% U Upper triangular matrix.

 

[R, C] = size(A);

for i = 1:R

L(i,1) = A(i,1);

U(i,i) = 1;

end

for j = 2:R

U(1,j)= A(1,j)/L(1,1);

end

for i = 2:R

for j = 2:i

L(i,j)=A(i,j)-L(i,1:j-1)*U(1:j-1,j);

end

for j = i+1:R

U(i,j)=(A(i,j)-L(i,1:i-1)*U(1:i-1,j))/L(i,i);

end

end

 

function [x,lu,piv] = GEpivot(A,b)

%

% function [x,lu,piv] = GEpivot(A,b)

%

% This program employs the Gaussian elimination method with partial

% pivoting to solve the linear system Ax=b.

%

% Input

% A: coefficient square matrix

% b: right side column vector

%

% Output

% x: solution vector

% lu: a matrix whose upper triangular part is the upper

% triangular matrix resulting from the Gaussian elimination

% with partial pivoting, and whose strictly lower triangular

% part stores the multipliers. In other words, it stores

% the LU factorization of the matrix A with row permutations

% determined by the pivoting vector piv.

% piv: a pivoting vector recording the row permutations.

 

% check the order of the matrix and the size of the vector

[m,n] = size(A);

if m ~= n

error('The matrix is not square.')

end

m = length(b);

if m ~= n

error('The matrix and the vector do not match in size.')

end

 

% initialization of the pivoting vector

piv = (1:n)';

 

% elimination step

for k = 1:n-1

% Find the maximal element in the pivot column, below the

% pivot position, along with the index of that maximal element.

 

[col_max index] = max(abs(A(k:n,k)));

index = index + k-1;

 

if index ~= k

% Switch rows k and index, in columns k thru n. Do similarly

% for the right-hand side b.

 

tempA = A(k,k:n);

A(k,k:n) = A(index,k:n);

A(index,k:n) = tempA;

 

tempb = b(k);

b(k) = b(index);

b(index) = tempb;

 

temp = piv(k);

piv(k) = piv(index);

piv(index) = temp;

end

 

% Form the needed multipliers and store them into the pivot

% column, below the diagonal.

A(k+1:n,k) = A(k+1:n,k)/A(k,k);

 

% Carry out the elimination step, first modifying the matrix,

% and then modifying the right-hand side.

for i = k+1:n

A(i,k+1:n) = A(i,k+1:n) - A(i,k)*A(k,k+1:n);

end

b(k+1:n) = b(k+1:n) - A(k+1:n,k)*b(k);

end

 

% Solve the upper triangular linear system.

x = zeros(n,1);

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

 

% Record the LU factorization with row permutation.

lu = A;

 

 

>> A=[4 3 2 1;3 4 3 2;2 3 4 3 ; 1 2 3 4];

>> b=[1 ;1 ;-1; -1];

>> [L, U]=LUdecompCrout(A)

L =

4 0 0 0

3 7/4 0 0

2 3/2 12/7 0

1 5/4 10/7 5/3

 

U =

1 3/4 1/2 1/4

0 1 6/7 5/7

0 0 1 5/6

0 0 0 1

 

>> y=GEpivot(L, b)

y =

1/4

1/7

-1

-1/7505999378950826

 

>> GEpivot(U, y)

ans =

1/9007199254740992

1

-1

-1/7505999378950826