1、待定系数法LU分解
设
A
=
L
U
A=LU
A=LU,即
[
a
11
a
12
a
13
⋯
a
1
n
a
21
a
22
a
23
⋯
a
2
n
a
31
a
32
a
33
⋯
a
3
n
⋮
⋱
a
n
1
a
n
2
a
n
3
⋯
a
n
n
]
\begin{bmatrix} a_{11}&a_{12}&a_{13}& \cdots&a_{1n} \\ a_{21}&a_{22}&a_{23}& \cdots&a_{2n} \\ a_{31}&a_{32}&a_{33}& \cdots&a_{3n} \\ \vdots&& \ddots&\\a_{n1}&a_{n2}&a_{n3}& \cdots&a_{nn} \\ \end{bmatrix}
⎣
⎡a11a21a31⋮an1a12a22a32an2a13a23a33⋱an3⋯⋯⋯⋯a1na2na3nann⎦
⎤=
[
1
l
21
1
l
31
l
32
1
⋮
⋱
⋮
l
n
1
l
n
2
⋯
l
n
,
n
−
1
1
]
\begin{bmatrix} 1&&&& \\ l_{21}&1&& & \\ l_{31}&l_{32}&1& & \\ \vdots&& \ddots&\vdots\\l_{n1}&l_{n2}& \cdots&l_{n,n-1}&1 \\ \end{bmatrix}
⎣
⎡1l21l31⋮ln11l32ln21⋱⋯⋮ln,n−11⎦
⎤
[
u
11
u
12
u
13
⋯
u
1
n
u
22
u
23
⋯
u
2
n
u
33
⋯
u
3
n
⋮
⋱
u
n
n
]
\begin{bmatrix} u_{11}&u_{12}&u_{13}&\cdots&u_{1n} \\ &u_{22}&u_{23}&\cdots & u_{2n}\\ &&u_{33}& \cdots&u_{3n} \\&&&\vdots&\ddots\\&&&&u_{nn} \\ \end{bmatrix}
⎣
⎡u11u12u22u13u23u33⋯⋯⋯⋮u1nu2nu3n⋱unn⎦
⎤
1、比较等式两边的第一行,可得
u
1
j
=
a
1
j
,
j
=
1
,
2
,
.
.
.
.
,
n
.
u_{1j}=a_{1j}, j=1,2,....,n.
u1j=a1j,j=1,2,....,n.
再比较等式两边的第一列,可得
a
i
1
=
l
i
1
u
11
⇒
l
i
1
=
a
i
1
/
u
11
,
i
=
2
,
3
,
.
.
.
,
n
a_{i1}=l_{i1}u_{11} \Rightarrow l_{i1}=a_{i1}/u_{11}, i=2,3,...,n
ai1=li1u11⇒li1=ai1/u11,i=2,3,...,n
2、比较等式两边的第二行,可得
a
2
j
=
l
21
u
1
j
+
u
2
j
⇒
u
2
j
=
a
2
j
−
l
21
u
1
j
,
j
=
2
,
3
,
.
.
.
,
n
a_{2j}=l_{21}u_{1j}+u_{2j}\Rightarrow u_{2j}=a_{2j}-l_{21}u_{1j}, j=2,3,...,n
a2j=l21u1j+u2j⇒u2j=a2j−l21u1j,j=2,3,...,n
再比较等式两边的第二列,可得
a
i
2
=
l
i
1
u
12
+
l
i
2
+
u
22
⇒
l
i
1
=
(
a
i
2
−
l
i
1
u
12
)
/
u
22
,
i
=
2
,
3
,
.
.
.
,
n
a_{i2}=l_{i1}u_{12}+l_{i2}+u_{22}\Rightarrow l_{i1}=(a_{i2}-l_{i1}u_{12})/u_{22}, i=2,3,...,n
ai2=li1u12+li2+u22⇒li1=(ai2−li1u12)/u22,i=2,3,...,n
3、依次类推比较等式两边的第k行,可得
u
k
j
j
=
a
k
j
−
(
l
k
1
u
1
j
+
⋯
+
l
k
,
k
−
1
u
k
−
1
,
j
)
,
j
=
k
,
k
+
1
,
.
.
.
,
n
u_{kjj}=a_{kj}-(l_{k1}u_{1j}+\cdots+l_{k,k-1}u_{k-1,j}), j=k,k+1,...,n
ukjj=akj−(lk1u1j+⋯+lk,k−1uk−1,j),j=k,k+1,...,n
依次类推比较等式两边的第k列,可得
l
i
k
=
(
a
i
k
−
l
i
1
u
1
k
−
⋯
−
l
i
,
k
−
1
u
k
−
1
,
k
)
/
u
k
k
,
i
=
k
,
k
+
1
,
.
.
.
,
n
l_{ik}=(a_{ik}-l_{i1}u_{1k}-\cdots-l_{i,k-1}u_{k-1,k})/u_{kk}, i=k,k+1,...,n
lik=(aik−li1u1k−⋯−li,k−1uk−1,k)/ukk,i=k,k+1,...,n
直到第n步算出L和U的所有元素。
2、总结LU分解公式
u
1
i
=
a
1
i
(
i
=
1
,
2
,
3
,
⋯
,
n
)
u_{1i}=a_{1i}(i=1,2,3,⋯,n)
u1i=a1i(i=1,2,3,⋯,n)
l
i
1
=
a
i
1
/
u
11
(
i
=
2
,
3
,
⋯
,
n
)
l_{i1}=a_{i1}/u_{11}(i=2,3,⋯,n)
li1=ai1/u11(i=2,3,⋯,n)
u
r
i
=
a
r
i
−
∑
k
=
1
r
−
1
l
r
k
u
k
i
(
i
=
r
,
r
+
1
,
r
+
2
,
⋯
,
n
)
u_{ri}=a_{ri}−\sum_{k=1} ^{r−1}l_{rk}u_{ki}(i=r,r+1,r+2,⋯,n)
uri=ari−k=1∑r−1lrkuki(i=r,r+1,r+2,⋯,n)
l
i
r
=
(
a
i
r
−
∑
k
=
1
r
−
1
l
i
k
u
k
r
)
/
u
r
r
(
i
=
r
+
1
,
⋯
,
n
;
且
r
≠
n
)
l_{ir}=(a_{ir}−\sum_{k=1}^{r−1}l_{ik}u_{kr})/u_{rr}(i=r+1,⋯,n;且r≠n)
lir=(air−k=1∑r−1likukr)/urr(i=r+1,⋯,n;且r=n)
C++实现方面,首先LU分解的函数传入两个参数,方阵的一阶数组和方阵的阶数(方阵用一维数组的行优先表示)。
计算步骤:
1. 初始化LU矩阵,L矩阵上三角为0,对角线为1,U矩阵下三角为0.
2. 计算U矩阵的第一行和L矩阵的第一列
3. 循环计算U矩阵和L矩阵,第一层循环表示U计算到第几行,同时也表示L计算到第几列,先计算U后计算L。第二层循环分别表示U矩阵改行的第几列元素。第三层循环就是公式中的求和符号部分。
下面展示一些 内联代码片
。
#include <iostream>
// 参数:一个order阶矩阵,和矩阵的阶数
void lowerUpperFactor(double *matrix, int order) {
printf("--------原矩阵:--------\n");
printMatrix(matrix, order,order);
// 结果变量 L矩阵和U矩阵都是order阶矩阵
double *L = new double[order*order];
double *U = new double[order*order];
// 初始化全为0
for (int i = 0; i < order; i++) {
// 初始化U下三角为0
for (int j = 0; j < i; j++) {
*(U + i * order + j) = 0;
}
//初始化L对角线为1,上三角为0
*(L + i * order + i) = 1;
for (int j = i + 1; j < order; j++) {
*(L + i * order + j) = 0;
}
}
// 计算U的第一行和L的第一列
int i = 0;
for (i = 0; i < order; i++) {
*(U + i) = *(matrix + i);
}
for (i = 1; i < order; i++) {
*(L + i * order) = *(matrix + i * order) / *U;
}
// 计算其余行列
int temp;
for (int i = 1; i < order; i++) {
// 计算矩阵U
for (int j = i; j < order; j++) {
temp = 0;
for (int k = 0; k < i; k++) {
temp+= (*(U + k * order + j) * (*(L + i * order + k)));
}
*(U + i * order + j) = *(matrix + i * order + j) - temp;
}
// 计算矩阵L
for (int j = i+1; j < order; j++) {
temp = 0;
for (int k = 0; k < i; k++) {
temp += *(U + k * order + i) * (*(L + j * order + k));
}
*(L + j * order + i) = (*(matrix +j * order + i) - temp) / (*(U+i* order + i));
}
}
printf("------矩阵U------\n");
printMatrix(U, order,order);
printf("------矩阵L------\n");
printMatrix(L, order, order);
if (L) {
delete[] L;
}
if (U) {
delete[] U;
}
}
void printMatrix(double *matrix, int row, int column) { for (int i = 0; i < row; i++) { for (int j = 0; j < column; j++) { printf("%6.2lf ", *(matrix + i * column + j)); } printf("\n"); }}
int main() {
//double matrix[] = { 1,2,3,2,5,2,3,1,5 };
//int order = 3;
double matrix1[] = { 1,1,-1,2,1,2,0,2,-1,-1,2,0,0,0,-1,1 };
int order1 = 4;
lowerUpperFactor(matrix1, order1);
}
3、方程求解
现在已经知道矩阵 A x = b Ax=b Ax=b的LU分解了,就要去求解方程的解了,已知方程 A x = b Ax=b Ax=b
矩阵A可做如下分解
A
=
L
U
A=LU
A=LU
则
A
x
=
b
⇔
L
U
x
=
b
Ax=b\Leftrightarrow LUx=b
Ax=b⇔LUx=b
令
U
x
=
y
Ux=y
Ux=y
则
A
x
=
b
⇔
L
y
=
b
⇒
y
=
L
−
1
b
U
x
=
y
⇒
x
=
U
−
1
y
Ax=b\Leftrightarrow Ly=b\Rightarrow y=L^{-1}b Ux=y\Rightarrow x=U^{-1}y
Ax=b⇔Ly=b⇒y=L−1bUx=y⇒x=U−1y
算法总共分为两步:
向前回代求解
L
U
x
=
b
LUx=b
LUx=b
向后回代求解
U
x
=
y
Ux=y
Ux=y
如果矩阵 A 不进行分解直接通过
A
−
1
b
A^{-1}b
A−1b 来求解方程运算量会比较大,因为通常来说 是不容易求的。而将 A 分解为 LU 的形式,单位下三角矩阵L 的逆和上三角矩阵 U 的逆是容易求的,因此很容易可以求出 y和x ,运算量将非常小。