前言
将高斯消元法改写为紧凑形式,可以直接从矩阵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版) 李庆扬,王能超,易大义 *编