LU分解(直接三角分解,Doolittle分解)
A
x
=
b
,
A
=
L
U
Ax=b \,\,,\,\, A=LU
Ax=b,A=LU
{
L
y
=
b
U
x
=
y
\begin{cases} Ly=b \\ Ux=y \end{cases}
{Ly=bUx=y
矩阵
L
{L}
L 的对角元素为
1
{1}
1 ,矩阵
U
{U}
U 的第一行和
A
{A}
A 相同。
步骤:
1.
矩阵
L
的对角元素为
1
,矩阵
U
的第一行和
A
相同。
2.
迭代
,
j
=
1
,
2
,
⋯
n
−
1
算
L
的第
j
列
,
L
i
,
j
=
A
i
,
j
−
∑
r
=
1
j
−
1
L
i
,
r
U
r
,
j
U
j
,
j
,
i
=
j
+
1
,
j
+
2
,
⋯
,
n
算
U
的第
j
+
1
行
,
U
j
+
1
,
k
=
A
j
+
1
,
k
−
∑
r
=
1
j
L
j
+
1
,
r
U
r
,
k
L
j
+
1
,
j
+
1
,
k
=
j
+
1
,
j
+
2
,
⋯
,
n
3.
回代
,
y
i
=
b
i
−
∑
j
=
1
i
−
1
L
i
,
j
y
j
,
i
=
1
,
2
,
⋯
,
n
x
i
=
y
i
−
∑
j
=
i
+
1
n
x
j
⋅
U
i
,
j
U
i
,
i
,
i
=
n
,
n
−
1
,
⋯
,
1
\begin{align*} 1.& 矩阵 L 的对角元素为 1 ,矩阵U 的第一行和A相同。 \\ \\ 2. & 迭代 \,\,,\,\, j=1,2, \cdots n-1 \\ \\ &算L的第j列 \,\,,\,\, L_{i,j}= \frac{A_{i,j}- \sum_{r=1}^{j-1}L_{i,r}U_{r,j}}{U_{j,j}},i=j+1,j+2,\cdots ,n \\ \\ &算U的第j+1行 \,\,,\,\, U_{j+1,k}= \frac{A_{j+1,k}- \sum_{r=1}^{ j}L_{j+1,r}U_{r,k}}{L_{j+1,j+1}} ,k=j+1,j+2,\cdots ,n \\ \\ 3.& 回代 \,\,,\,\, \\ \\ & y_i= b_i- \sum_{j=1}^{ i-1}L_{i,j}y_j,i=1,2,\cdots ,n \\ \\ &x_i= \frac{y_i- \sum_{j=i+1}^{ n}x_j \cdot U_{i,j}}{U_{i,i}} \,\,,\,\, i=n,n-1, \cdots ,1 \end{align*}
1.2.3.矩阵L的对角元素为1,矩阵U的第一行和A相同。迭代,j=1,2,⋯n−1算L的第j列,Li,j=Uj,jAi,j−∑r=1j−1Li,rUr,j,i=j+1,j+2,⋯,n算U的第j+1行,Uj+1,k=Lj+1,j+1Aj+1,k−∑r=1jLj+1,rUr,k,k=j+1,j+2,⋯,n回代,yi=bi−j=1∑i−1Li,jyj,i=1,2,⋯,nxi=Ui,iyi−∑j=i+1nxj⋅Ui,j,i=n,n−1,⋯,1
matlab实现
%% Ax=b例子
A = [16 -12 2 4;
12 -8 6 10;
3 -13 9 23;
-6 14 1 -28];
b = [17 36 -49 -54]';
[x,L,U] = LUsolve(A,b)
%% LU分解解Ax=b
% 输入方阵A,向量b
% 输出解x,L、U矩阵
function [x,L,U] = LUsolve(A,b)
n = size(A);
L = eye(n);
U(1,[1:n]) = A(1,[1:end]);
for j = 1:n-1 % 对U是行号,对L是列号
for i = j+1:n % 算L第i行j列
L(i,j) = A(i,j);
for r = 1:j-1
L(i,j) = L(i,j)- L(i,r)*U(r,j);
end
L(i,j) = L(i,j)/U(j,j);
end
for k = j+1:n % 算U第j+1行k列
U(j+1,k) = A(j+1,k);
for r = 1:j
U(j+1,k) = U(j+1,k)-L(j+1,r)*U(r,k);
end
U(j+1,k) = U(j+1,k)/L(j+1,j+1);
end
end
% 回代
for i = 1:n
y(i) = b(i);
for j = 1:i-1
y(i) = y(i)-L(i,j)*y(j);
end
end
for i=n:-1:1
x(i) = y(i);
for j=n:-1:i+1
x(i) = x(i)-U(i,j)*x(j);
end
x(i) = x(i)/U(i,i);
end
x = x';
end