-etapu eliminacji "wprzód", w którym doprowadzamy do macierzy trójkątnej górnej \(\displaystyle{ U, }\)
-etapu " podstawienia "od tyłu", w którym uzyskujemy macierz trójkątną dolną \(\displaystyle{ L. }\)
Mamy więc rozkład macierzy układu \(\displaystyle{ A }\) na iloczyn macierzy \(\displaystyle{ L\cdot U. }\)
Te dwuetapowe rozwiązywanie układu równań możemy przedstawić w postaci:
\(\displaystyle{ L\cdot U\cdot \vec{x} = \vec{b}: }\)
\(\displaystyle{ U\cdot \vec{x} = \vec{y}, }\)
\(\displaystyle{ L\cdot \vec{y} = \vec{b}. }\)
Program w OCTAVE
Kod: Zaznacz cały
function lufact(A,b)
% Rozwiązywanie układu Ax=b, używając rozkładu LU.
n=length(b);
y=zeros(n,1);
x=zeros(n,1);
fprintf('\n');
for i=1:n
U(i,i)=1;
end
L(1,1)=A(1,1)/U(1,1);
for j=2:n
L(j,1)=A(j,1)/U(1,1);
U(1,j)=A(1,j)/L(1,1);
end
for i=2:n-1
S=0;
for k=1:i-1
S=S+U(k,i)*L(i,k);
end
L(i,i)=(A(i,i)-S)/U(i,i);
for j=i+1:n
S=0;
for k=1:i-1
S=S+U(k,i)*L(j,k);
end
L(j,i)=(A(j,i)-S)/U(i,i);
S=0;
for k=1:i-1
S=S+U(k,j)*L(i,k);
end
U(i,j)=(A(i,j)-S)/L(i,i);
end
end
S=0;
for k=1:n-1
S=S+U(k,n)*L(n,k);
end
L(n,n)=(A(n,n)-S)/U(n,n)
% Podstawienie "wprzód"
y(1)=b(1)/L(1,1);
for i=2:n
S=b(i);
for j=1:i-1
S=S-L(i,j)*y(j);
end
y(i)=S/L(i,i);
end
%Podstawienie " od tyłu"
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
S=y(i);
for j=i+1:n
S=S-U(i,j)*x(j);
end
x(i)=S/U(i,i);
end
L
disp(' Podstawienie "wprzód" daje ')
y
U
disp(' Wektorem rozwiązań jest =')
x
Kod: Zaznacz cały
>> A=[1 1 1 1;2 3 1 5;-1 1 -5 3;3 1 7 -2];
>> b=[10 31 -2 18]';
L =
1 0 0 0
2 1 0 0
-1 2 -2 0
3 -2 2 -1
Podstawienie "wprzód" daje:
y =
10
11
7
4
U =
1 1 1 1
0 1 -1 3
0 0 1 1
0 0 0 1
Wektorem rozwiązań jest =
x =
1
2
3
4