前言
将高斯消元法改写为紧凑形式，可以直接从矩阵A的元素的得到计算L，U元素的递推公式，而不需要任何中间步骤，这就是直接三角分解法。这样我们就可以把求解$Ax = b$的问题转化成求解两个三角形方程组的问题：
（1）Ly = b，求y；
（2）Ux = y，求x。
算法分析
对于一个非奇异矩阵A， 有分解式
$$
A = LU
$$
，其中L为单位下三角矩阵，U为上三角矩阵，即
$$
A=
\begin{bmatrix}
1 \
l_{21} & 1 \
\vdots & \vdots & \ddots \
l_{n1} & l_{n2} & \cdots & 1\
\end{bmatrix}
\begin{bmatrix}
u_{11} & u_{12} & \cdots & u_{1n} \
& u_{22} & \cdots & u_{2n} \
& & \ddots & \vdots \
& & & u_{nn}\
\end{bmatrix}
$$
上面只是我们的猜想，现在我们就来说明L，U的元素可以由$n$步直接计算定出，其中第$r$步定出U的第$r$行和L的第$r$列元素：
$$
a_{1i} = u_{1i}, i=1,2,…,n
$$
得U的第1行元素；
$$
a_{i1} = l_{i1}u_{11}, l_{i1} = a_{i1} / u_{i1}, i=2,3,…,n
$$
得L的第1列元素。
设已经定出U的第1行到第$r-1$行元素与L的第$1$列到第$r-1$列元素，由分解式得到：
$$
a_{ri} = \sum_{k=1}^nl_{rk}u_{ki} = \sum_{k=1}^{r-1}l_{rk}u_{ki} + u_{ri},
$$
故
$$
u_{ri} = a_{ri} - \sum_{k=1}^{r-1}, i=r,r+1,…,n
$$
所以，
$$
a_{ir} = \sum_{k=1}^nl_{ik}u_{kr} = \sum_{k=1}^{k-1}l_{ik}x_{kr} + l_{ir}u_{rr}
$$
综上所述，分解的计算公式为：
$$① u_{1i} = a_{1i}, l_{i1} = a_{i1} / u_{11}, i=2,3,…,n$$
$$② u_{ri} = a_{ri} - \sum_{k=1}^{r-1}l_{rk}u_{ki}, i=r,r+1,…,n$$
$$③ l_{ir} = (a_{ir} - \sum_{k=1}^{r-1}l_{ik}u_{kr} / u_{rr}, i=r+1,…,n, 且r≠n
.$$
所以求解Ly=b和Ux = y的计算公式为：
$$
\begin{cases}
y_1 = b_1\
y_i = b_i - \sum_{k=1}^{i-1}l_{ik}y_k, i=2,3,…,n;\
\end{cases}
$$
$$
\begin{cases}
x_n = y_n / u_{nn}\
x_i = (y_i - \sum_{k=i+1}^{n}l_{ik}x_k) / u_{ii}, i=n-1, n-2,…,1\
\end{cases}
$$
代码分析
变量定义
1
2
3
4
5
6
7
| int n;
double A[maxn][maxn] = {0};
double L[maxn][maxn] = {0};
double U[maxn][maxn] = {0};
double b[maxn] = {0};
double y[maxn] = {0};
double x[maxn] = {0};
|
使用二维数组存储矩阵，使用一维数组存储向量。
输入与输出
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
|
void read() {
cout << "Please input the scale of matrix n: ";
cin >> n;
cout << "|--------------------\n";
cout << "|Please input the data: " << endl;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cin >> A[i][j];
}
}
for (int i = 1; i <= n; i++) {
cin >> b[i];
}
}
void printLU() {
cout << "L matrix: " << endl;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cout << setw(10) << L[i][j] << " ";
}
cout << endl;
}
cout << "U matrix: " << endl;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cout << setw(10) << U[i][j] << " ";
}
cout << endl;
}
}
void print() {
cout << "|-------------------" << endl;
for (int i = 1; i <= n; i++) {
cout << "y[" << i << "] = " << y[i] << endl;
}
cout << "|-------------------" << endl;
for (int i = 1; i <= n; i++) {
cout << "x[" << i << "] = " << x[i] << endl;
}
}
|
分别有处理输入、输出分解后结果和输出最终解的三个方法。
LU分解
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
| void LUDecomposition() {
int i, r, k;
for (i = 1; i <= n; i++) {
U[1][i] = A[1][i];
}
for (i = 2; i <= n; i++) {
L[i][1] = A[i][1] / U[1][1];
}
for (r = 2; r <= n; r++) {
for (i = r; i <= n; i++) {
double sum1 = 0;
for (k = 1; k < r; k++) {
sum1 += L[r][k] * U[k][i];
}
U[r][i] = A[r][i] - sum1;
}
if (r != n) {
for (i = r + 1; i <= n; i++) {
double sum2 = 0;
for (k = 1; k < r; k++) {
sum2 += L[i][k] * U[k][r];
}
L[i][r] = (A[i][r] - sum2) / U[r][r];
}
}
}
for (int i = 1; i <= n; i++) {
L[i][i] = 1;
}
y[1] = b[1];
for (i = 2; i <= n; i++) {
double sum3 = 0;
for (k = 1; k < i; k++) {
sum3 += L[i][k] * y[k];
}
y[i] = b[i] - sum3;
}
x[n] = y[n] / U[n][n];
for (i = n - 1; i >= 1; i--) {
double sum4 = 0;
for (k = i + 1; k <= n; k++) {
sum4 += U[i][k] * x[k];
}
x[i] = (y[i] - sum4) / U[i][i];
}
printLU();
print();
}
|
按照上述的算法进行编程，最后输出LU矩阵和x。
完整代码
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
| #include <iostream>
#include <iomanip>
using namespace std;
#define maxn 50
int n;
double A[maxn][maxn] = {0};
double L[maxn][maxn] = {0};
double U[maxn][maxn] = {0};
double b[maxn] = {0};
double y[maxn] = {0};
double x[maxn] = {0};
void read() {
cout << "Please input the scale of matrix n: ";
cin >> n;
cout << "|--------------------\n";
cout << "|Please input the data: " << endl;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cin >> A[i][j];
}
}
for (int i = 1; i <= n; i++) {
cin >> b[i];
}
}
void printLU() {
cout << "L matrix: " << endl;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cout << setw(10) << L[i][j] << " ";
}
cout << endl;
}
cout << "U matrix: " << endl;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cout << setw(10) << U[i][j] << " ";
}
cout << endl;
}
}
void print() {
cout << "|-------------------" << endl;
for (int i = 1; i <= n; i++) {
cout << "y[" << i << "] = " << y[i] << endl;
}
cout << "|-------------------" << endl;
for (int i = 1; i <= n; i++) {
cout << "x[" << i << "] = " << x[i] << endl;
}
}
void LUDecomposition() {
int i, r, k;
for (i = 1; i <= n; i++) {
U[1][i] = A[1][i];
}
for (i = 2; i <= n; i++) {
L[i][1] = A[i][1] / U[1][1];
}
for (r = 2; r <= n; r++) {
for (i = r; i <= n; i++) {
double sum1 = 0;
for (k = 1; k < r; k++) {
sum1 += L[r][k] * U[k][i];
}
U[r][i] = A[r][i] - sum1;
}
if (r != n) {
for (i = r + 1; i <= n; i++) {
double sum2 = 0;
for (k = 1; k < r; k++) {
sum2 += L[i][k] * U[k][r];
}
L[i][r] = (A[i][r] - sum2) / U[r][r];
}
}
}
for (int i = 1; i <= n; i++) {
L[i][i] = 1;
}
y[1] = b[1];
for (i = 2; i <= n; i++) {
double sum3 = 0;
for (k = 1; k < i; k++) {
sum3 += L[i][k] * y[k];
}
y[i] = b[i] - sum3;
}
x[n] = y[n] / U[n][n];
for (i = n - 1; i >= 1; i--) {
double sum4 = 0;
for (k = i + 1; k <= n; k++) {
sum4 += U[i][k] * x[k];
}
x[i] = (y[i] - sum4) / U[i][i];
}
printLU();
print();
}
int main() {
while (1) {
read();
LUDecomposition();
}
}
|
LU分解法求解$Ax = b$的算法复杂度大约是O($\frac{1}{3} n^3$)**，但是如果对于不同b来说，一旦完成了分解，之后的每一次求解复杂度仅为O($n^2$)**，这是LU分解相比起高斯消元法有优势的地方。
至此，矩阵的LU分解算法及代码的讲解已经结束了，谢谢！
参考资料：
1.数值分析（第5版） 李庆扬，王能超，易大义 *编
