Descompunere LUP

-eliminare gaussiana

Contents

Codul pentru sisteme triunghiulare

function [L,U,P]=lup(A)
%LUP find LUP decomposition of matrix A
%call [L,U,P]=lup(A)
%permute effectively lines
[m,n]=size(A);
P=zeros(m,n);
piv=(1:m)';
for i=1:m-1
   %pivoting
   [pm,kp]=max(abs(A(i:m,i)));
   kp=kp+i-1;
   %line interchange
   if i~=kp
       A([i,kp],:)=A([kp,i],:);
       piv([i,kp])=piv([kp,i]);
   end
   %Schur complement
   lin=i+1:m;
   A(lin,i)=A(lin,i)/A(i,i);
   A(lin,lin)=A(lin,lin)-A(lin,i)*A(i,lin);
end;
for i=1:m
   P(i,piv(i))=1;
end;
U=triu(A);
L=tril(A,-1)+eye(m);

Codul pentru sisteme triunghiulare

function x=backsubst(U,b)
%BACKSUBST - backward substitution
%U - upper triangular matrix
%b - right hand side vector
n=length(b);
x=zeros(size(b));
for k=n:-1:1
   x(k)=(b(k)-U(k,k+1:n)*x(k+1:n))/U(k,k);
end
function x=forwardsubst(L,b)
%FORWARDSUBST - forward substitution
%L - lower triangular matrix
%b - right hand side vector
x=zeros(size(b));
n=length(b);
for k=1:n
   x(k)=(b(k)-L(k,1:k-1)*x(1:k-1))/L(k,k);
end
%
A = [10,7,0; -3, 2, 6; 5, -1, 5];
b = [7; 4; 6];

Exemplu de utilizare

Descompunere

[L,U,P] = lup(A)
norm(L*U-P*A)
L =

   1.0000e+00            0            0
   5.0000e-01   1.0000e+00            0
  -3.0000e-01  -9.1111e-01   1.0000e+00


U =

   1.0000e+01   7.0000e+00            0
            0  -4.5000e+00   5.0000e+00
            0            0   1.0556e+01


P =

     1     0     0
     0     0     1
     0     1     0


ans =

   6.6613e-16

function x=lupsolve(A,b)
%LUPSOLVE - solution of an algebraic system by LUP decomposition
[L,U,P]=lup(A);
y=forwardsubst(L,P*b);
x=backsubst(U,y);

rezolvare efectiva

x = lupsolve(A,b)
x =

   4.7158e-01
   3.2632e-01
   7.9368e-01

Exemplul de pe slide-uri

A = [1, 1, 1; 1, 1, 2; 2, 4, 2];
b = [3,4,8]';
[L,U,P] = lup(A)
y = forwardsubst(L,P*b)
x= backsubst(U,y)
L =

   1.0000e+00            0            0
   5.0000e-01   1.0000e+00            0
   5.0000e-01   1.0000e+00   1.0000e+00


U =

     2     4     2
     0    -1     1
     0     0    -1


P =

     0     0     1
     0     1     0
     1     0     0


y =

     8
     0
    -1


x =

     1
     1
     1