标签：right matrix int sum 线性方程组 cdots 直接 解法 left
三角形方程组和三角分解
前代法
求解下三角形方程组
\[Ly = b \]其中 \(b=(b_1,\cdots,b_n)^T\in\mathbb{R}^n\) 已知， \(y=(y_1,\cdots,y_n)^T\in\mathbb{R}^n\) 未知，而
\[L = \left( \begin{matrix} l_{11}\\ l_{21} & l_{22}\\ \vdots & \vdots & \ddots\\ l_{n1} & l_{n2} & \cdots & l_{nn} \end{matrix} \right) \]是已知的非奇异下三角阵，且 \(l_{ii}\neq 0,\ i=1,\cdots,n\) .
前代法通过公式
\[y_i = \left(b_i-\sum_{j=1}^{i-1}l_{ij}y_j\right)\bigg{/}l_{ii} \]从上到下求解方程组.
代码实现
vector<double> lower(vector<vector<double>> &A, vector<double> &b)
{
int n = b.size();
vector<double> x(n);
double sum;
x[0] = b[0];
for (int i = 1; i < n; i++)
{
sum = 0;
for (int j = 0; j < i; j++)
{
sum += A[i][j] * x[j];
}
x[i] = (b[i] - sum) / A[i][i];
}
return x;
}
回代法
求解上三角形方程组
\[Ux = y \]其中 \(y=(y_1,\cdots,y_n)^T\in\mathbb{R}^n\) 已知， \(x=(x_1,\cdots,x_n)^T\in\mathbb{R}^n\) 未知，而
\[U = \left( \begin{matrix} u_{11} & u_{12} & \cdots & u_{1n}\\ & u_{22} & \cdots & u_{2n}\\ & & \ddots & \vdots\\ & & & u_{nn}\\ \end{matrix} \right) \]是已知的非奇异上三角阵，且 \(u_{ii}\neq 0,\ i=1,\cdots,n\) .
回代法通过公式
从下到上求解方程组.
代码实现
vector<double> upper(vector<vector<double>> &A, vector<double> &b)
{
int n = b.size();
vector<double> x(n);
double sum;
for (int i = n - 1; i >= 0; i--)
{
sum = 0;
for (int j = i + 1; j < n; j++)
{
sum += A[i][j] * x[j];
}
x[i] = (b[i] - sum) / A[i][i];
}
return x;
}
Gauss 变换
通过一系列初等变换，逐步用下三角阵对 \(A\) 进行约化
\[L_k = I - l_ke_k^T,\quad l_k = (0,\cdots,0,l_{k+1,k},\cdots, l_{nk})^T \]左乘 \(L_k\) 进行行变换
\[L_k = \left( \begin{matrix} 1\\ & \ddots\\ & & 1\\ & & -l_{k+1,k} & 1\\ & & \vdots & & \ddots\\ & & -l_{nk}& & & 1 \end{matrix} \right) \]这种类型的初等下三角阵称为 Gauss 变换，向量 \(l_k\) 称为 Gauss 向量.
高斯向量具有重要性质，因为 \(e_k^Tl_k = 0\) ，从而
\[(I-l_ke_k^T)(I+l_ke_k^T) = I - l_ke_k^Tl_ke_k^T = I \]即有
\[L_k^{-1} = I + l_ke_k^T \]高斯变换的逆容易得到.
三角分解的计算
设 \(A\in\mathbb{R}^{n\times n}\) ，则 \(A\) 的三角分解指分解 \(A=LU\) ，其中 \(L\in\mathbb{R}^{n\times n}\) 为下三角阵， \(U\in\mathbb{R}^{n\times n}\) 为上三角阵，此分解也称为 LU 分解.
我们通过左乘一系列高斯变换来消去 \(A\) 的下三角部分
\[A^{(k-1)} = L_{k-1}\cdots L_1A = \left( \begin{matrix} A_{11}^{(k-1)} & A_{12}^{(k-1)}\\ 0 & A_{22}^{(k-1)} \end{matrix} \right) \]其中 \(A_{11}^{(k-1)}\) 是 \(k-1\) 阶上三角阵，而
\[A_{22}^{(k-1)} = \left( \begin{matrix} a_{kk}^{(k-1)} & \cdots & a_{kn}^{(k-1)}\\ \vdots & \ddots & \vdots\\ a_{nk}^{(k-1)} & \cdots & a_{nn}^{(k-1)} \end{matrix} \right) \]如果 \(a_{kk}^{(k-1)}\neq 0\) ，则可以再确定一个高斯变换 \(L_k = I-l_ke_k^T,\ l_{ki} = a_{ki}^{(k-1)}/a_{kk}^{(k-1)},\ i=k+1,\cdots,n\) .
最终我们得到
\[L= (L_{n-1}L_{n-2}\cdots L_1)^{-1},\quad U = A^{(n-1)},\quad A = LU \]利用高斯变换的特点
\[L = (L_{n-1}L_{n-2}\cdots L_1)^{-1} = L_1^{-1}\cdots L_{n-1}^{-1} = (I+l_1e_1^T)\cdots (I+l_{n-1}e_{n-1}^T) = I + l_1e_1^T + \cdots l_{n-1}e_{n-1}^T \]即 \(L\) 有如下形状
\[L = \left( \begin{matrix} 1\\ l_{21} & 1\\ \vdots & \vdots & \ddots\\ l_{n1} & l_{n2} & \cdots & 1 \end{matrix} \right) \]此方法称为 Gauss 消去法；通常称 \(a_{kk}^{(k-1)}\) 为主元.
高斯消去法的存储
由于每次消元都与之前的列无关，并且下三角阵对角元为 \(1\) ，可以省略，因此可以将每一步高斯变换的结果直接存放在原先的列中，我们只展示最后的结果
\[A^{(n-1)} = \left( \begin{matrix} u_{11} & u_{12} & \cdots & u_{1n}\\ l_{21} & u_{22} & \cdots & u_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ l_{n1} & l_{n2} & \cdots & u_{nn} \end{matrix} \right) \]高斯消去法是原地消去，不需要额外的存储空间.
代码实现
void LU(vector<vector<double>> &A)
{
int n = A.size();
for (int i = 0; i < n; i++)
{
double key = A[i][i];
// 无法进行高斯变换
if (key == 0)
{
return;
}
// 第 i 列高斯变换
for (int j = i + 1; j < n; j++)
{
A[j][i] = A[j][i] / key;
for (int k = i + 1; k < n; k++)
{
A[j][k] -= A[i][k] * A[j][i];
}
}
}
}
定理 1.1.1 主元 \(a_{ii}^{(i-1)},\ i=1,\cdots,k\) 均不为零的充要条件是 \(A\) 的 \(i\) 阶顺序主子式 \(A_i\) 非奇异.
证明 容易应用归纳法证明.
定理 1.1.2 若 \(A\in\mathbb{R}^{n\times n}\) 的顺序主子式 \(A_k\in\mathbb{R}^{k\times k},\ k=1,\cdots,n-1\) 非奇异，则存在唯一分解 \(A=LU\) .
选主元三角分解
由于计算机精度的限制，如果主元太小，则会导致高斯变换的误差增大，因此我们需要选取尽可能小的元素作为主元.
\[A^{(k-1)} = \left( \begin{matrix} A_{11}^{(k-1)} & A_{12}^{(k-1)}\\ 0 & A_{22}^{(k-1)} \end{matrix} \right) \]我们在每一步高斯变换前，进行初等变换 \(I_{pq}\) ，左乘 \(I_{pq}\) 会交换 \(p,q\) 行，右乘 \(I_{pq}\) 会交换 \(p,q\) 列。选取
\[\left|a_{pq}^{(k-1)}\right| = \max\left\{\left|a_{ij}^{(k-1)}\right|:k\le i,j\le n\right\} \]如果 \(a_{pq}^{(k-1)} = 0\) ，说明 \(A\) 秩为 \(k-1\) ，结束过程；否则利用 \(I_{kp}\) 与 \(I_{kq}\) 交换 \(a_{kk}\) 与 \(a_{pq}\) ，最终得到
\[L_rP_r\cdots L_1P_1AQ_1\cdots Q_r = U \]其中 \(P_k = I_{kp},\ Q_k={I_{kq}}\) ，此过程称为全主元 Gauss 消去法.
令
\[\begin{aligned} Q &= Q_1\cdots Q_r\\ P &= P_r\cdots P_1\\ L &= P(L_rP_r\cdots L_1P_1)^{-1} \end{aligned} \]则有
\[PAQ = LU \]称为全主元三角分解；注意到 \(P_k^{-1}=P_k\) ，有
\[L = P_r\cdots P_1\cdot (L_rP_r\cdots L_1P_1)^{-1} = P_r\cdots P_2L_1^{-1}P_2L_2^{-1}\cdots P_rL_r^{-1} \]如果我们以 \(L_1^{-1}\) 为中心，令
\[L^{(1)} = L_1^{-1},\quad L^{(k)} = P_kL^{(k-1)}P_kL_k^{-1} \]其中的每一步的 \(P_kL^{(k-1)}P_kL_k^{-1}\) 只改变 \(I_{n-k+1}\) 而不会改变下三角结构
\[L^{(k-1)} = \left( \begin{matrix} L_{11}^{(k-1)} & 0\\ L_{21}^{(k-1)} & I_{n-k+1} \end{matrix} \right) \]因此结果 \(L\) 仍为单位下三角阵；并且由于每一步选取最大的主元，由高斯变换的构造得到 \(L\) 中元素的模不会超过 \(1\) .
定理 1.2.1 若 \(A\in\mathbb{R}^{n\times n}\) ，则有排列方阵 \(P,Q\in\mathbb{R}^{n\times n}\) 使得
\[PAQ = LU \]其中 \(L\) 中所有元素的模不大于 \(1\) ， \(U\) 中非零对角元的个数恰为 \(A\) 的秩.
虽然全主元高斯消去法弥补了不选主元的高斯消去法的不足，但仍需要付出极大的代价，因为在 \(A\) 非奇异的情况下，选主元需进行
\[\sum_{k=1}^{n-1}(n-k+1)^2 = \dfrac{1}{3}n^3+O(n^2) \]次比较判断。于是我们考虑列主元高斯消去法
\[\left|a_{pk}^{(k-1)}\right| = \max\left\{\left|a_{ik}^{(k-1)}\right|:k\le i\le n\right\} \]这样只进行每一列的选主元，能够大量减少时间.
代码实现
// 我们同时进行对 A, b 的变换
void column_lu(vector<vector<double>> &A, vector<double> &b)
{
int n = A.size();
double max, tmp; // 记录最大值
int index; // 记录要交换的主元位置
for (int i = 0; i < n; i++)
{
max = 0;
index = i;
// 选取列主元
for (int p = i; p < n; p++)
{
tmp = fabs(A[p][i]);
if (tmp > max)
{
max = tmp;
index = p;
}
}
// 在这里奇异，说明 A 的秩为 i
if (max == 0)
{
cout << "Singular Matrix!" << endl;
return;
}
// 交换主元
for (int q = 0; q < n; q++)
{
tmp = A[i][q];
A[i][q] = A[index][q];
A[index][q] = tmp;
}
// 同时交换 b
tmp = b[i];
b[i] = b[index];
b[index] = tmp;
// 正常的高斯消去法
for (int j = i + 1; j < n; j++)
{
A[j][i] = A[j][i] / A[i][i];
for (int k = i + 1; k < n; k++)
{
A[j][k] = A[j][k] - A[i][k] * A[j][i];
}
}
}
}
平方根法
平方根法又叫做 Cholesky 分解法，是求解对称正定线性方程组最常用的方法之一.
定理 1.3.1 (Cholesky 分解定理) 若 \(A\in\mathbb{R}^{n\times n}\) 对称正定，则存在对角元均为正数的下三角阵 \(L\in\mathbb{R}^{n\times n}\) ，使得
\[A = LL^T \]称为 Cholesky 分解，其中 \(L\) 称为 \(A\) 的 Cholesky 因子.
证明 由于对称正定，因此其全部主子阵均正定，则存在 \(\widetilde{L},U\) 使得 \(A = \widetilde{L}U\) ，令
\[D = \mathrm{diag}(u_{11},\cdots,u_{nn}),\quad \widetilde{U} = D^{-1}U \]则有
\[\widetilde{U}^TD\widetilde{L}^T = A^T = A = \widetilde{L}D\widetilde{U}\ \Rightarrow\ \widetilde{L}^T\widetilde{U}^{-1} = D^{-1}\widetilde{U}^{-T}\widetilde{L}D \]注意到左边是单位上三角阵，右边是下三角阵，因此两边均为单位矩阵，于是 \(\widetilde{U} = \widetilde{L}^T,\ A=\widetilde{L}D\widetilde{L}^T\) ；由于 \(D\) 对角元为正数，则取
\[L = \widetilde{L}\mathrm{diag}(\sqrt{u_{11}},\cdots,\sqrt{u_{nn}}) \]有 \(A = LL^T\) ，此即为满足要求的分解.
平方根法通过比较元素来计算分解，对比 \(A\) 和 \(LL^T\) 两边的对应元素得到
\[l_{kk} = \left(a_{ik}-\sum_{p=1}^{k-1}l_{kp}^2\right)^{\frac{1}{2}},\quad l_{ik} = \left(a_{ik}-\sum_{p=1}^{k-1}l_{ip}l_{kp}\right)\bigg{/}l_{kk} \]从左向右，每次计算一列元素，先算出对角元素，然后推导其余元素.
代码实现
void cholesky(vector<vector<double>> &A)
{
double sum = 0;
int n = A.size();
for (int i = 0; i < n; i++)
{
sum = 0;
for (int p = 0; p < i; p++)
{
sum += A[i][p] * A[i][p];
}
A[i][i] = sqrt(A[i][i] - sum);
// 奇异矩阵
if (A[i][i] == 0)
{
return;
}
for (int j = i + 1; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
{
sum += A[j][k] * A[i][k];
}
A[j][i] = (A[j][i] - sum) / A[i][i];
A[i][j] = A[j][i];
}
}
}
为了避免开方运算，我们使用改进的平方根法
\[A = LDL^T \]其中 \(L\) 是单位下三角阵， \(D\) 是对角元均为正的对角阵.
仍然通过对比元素计算分解
\[d_j = a_{jj} - \sum_{k=1}^{j-1}l_{jk}d_kl_{jk},\quad l_{ij} = \left(a_{ij}-\sum_{k=1}^{j-1}l_{ik}d_{k}l_{jk}\right)\bigg{/}d_{j} \]在实际计算时，将 \(L\) 的严格下三角元素存储在 \(A\) 的对应位置， \(D\) 的对角元存储在 \(A\) 的对应对角位置.
代码实现
void cholesky_d(vector<vector<double>> &A)
{
double sum = 0;
int n = A.size();
for (int i = 0; i < n; i++)
{
sum = 0;
for (int p = 0; p < i; p++)
{
sum += A[i][p] * A[i][p] * A[p][p];
}
A[i][i] -= sum;
// 奇异矩阵
if (A[i][i] == 0)
{
return;
}
for (int j = i + 1; j < n; j++)
{
sum = 0;
for (int k = 0; k < i; k++)
{
sum += A[j][k] * A[i][k] * A[k][k];
}
A[i][j] = A[j][i] - sum;
A[j][i] = A[i][j] / A[i][i];
}
}
}
追赶法
三对角阵的分解
\[\left( \begin{matrix} a_{1} & b_2\\ c_{2} & a_{2} & b_3\\ & \ddots & \ddots & \ddots\\ & & c_{n-1} & a_{n-1} & b_{n}\\ & & & c_{n} & a_n \end{matrix} \right) = \left( \begin{matrix} \alpha_{1}\\ \gamma_{2} & \alpha_{2}\\ & \ddots & \ddots\\ & & \gamma_{n-1} & \alpha_{n-1} \\ & & & \gamma_{n} & \alpha_n \end{matrix} \right) \left( \begin{matrix} 1 & \beta_2\\ & 1 & \beta_3\\ & & \ddots & \ddots\\ & & & 1 & \beta_{n}\\ & & & & 1 \end{matrix} \right) \]对比元素得到 \(\alpha_1 = a_1\) ，而
\[\gamma_i = c_i,\quad \beta_i = b_i/\alpha_{i-1},\quad \alpha_i = a_i-\gamma_i\beta_i,\quad i=2,\cdots,n \]此方法称为追赶法.
标签：right,matrix,int,sum,线性方程组,cdots,直接,解法,left 来源： https://www.cnblogs.com/Bluemultipl/p/15913455.html
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。