线性方程组的直接解法
从这一章开始是数值代数课程的笔记，使用的教材是徐树方、高立、张平文编著的《数值线性代数》。
三角形方程组和三角分解
前代法
求解下三角形方程组
其中 \(b=(b_1,\cdots,b_n)^T\in\mathbb{R}^n\) 已知， \(y=(y_1,\cdots,y_n)^T\in\mathbb{R}^n\) 未知，而
是已知的非奇异下三角阵，且 \(l_{ii}\neq 0,\ i=1,\cdots,n\) .
前代法通过公式
从上到下求解方程组.
代码实现
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;
}
回代法
求解上三角形方程组
其中 \(y=(y_1,\cdots,y_n)^T\in\mathbb{R}^n\) 已知， \(x=(x_1,\cdots,x_n)^T\in\mathbb{R}^n\) 未知，而
是已知的非奇异上三角阵，且 \(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\) 进行行变换
这种类型的初等下三角阵称为 Gauss 变换，向量 \(l_k\) 称为 Gauss 向量.
高斯向量具有重要性质，因为 \(e_k^Tl_k = 0\) ，从而
即有
高斯变换的逆容易得到.
三角分解的计算
设 \(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_{11}^{(k-1)}\) 是 \(k-1\) 阶上三角阵，而
如果 \(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\) 有如下形状
此方法称为 Gauss 消去法；通常称 \(a_{kk}^{(k-1)}\) 为主元.
高斯消去法的存储
由于每次消元都与之前的列无关，并且下三角阵对角元为 \(1\) ，可以省略，因此可以将每一步高斯变换的结果直接存放在原先的列中，我们只展示最后的结果
高斯消去法是原地消去，不需要额外的存储空间.
代码实现
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\) .
选主元三角分解
由于计算机精度的限制，如果主元太小，则会导致高斯变换的误差增大，因此我们需要选取尽可能小的元素作为主元.
我们在每一步高斯变换前，进行初等变换 \(I_{pq}\) ，左乘 \(I_{pq}\) 会交换 \(p,q\) 行，右乘 \(I_{pq}\) 会交换 \(p,q\) 列。选取
如果 \(a_{pq}^{(k-1)} = 0\) ，说明 \(A\) 秩为 \(k-1\) ，结束过程；否则利用 \(I_{kp}\) 与 \(I_{kq}\) 交换 \(a_{kk}\) 与 \(a_{pq}\) ，最终得到
其中 \(P_k = I_{kp},\ Q_k={I_{kq}}\) ，此过程称为全主元 Gauss 消去法.
令
则有
称为全主元三角分解；注意到 \(P_k^{-1}=P_k\) ，有
如果我们以 \(L_1^{-1}\) 为中心，令
其中的每一步的 \(P_kL^{(k-1)}P_kL_k^{-1}\) 只改变 \(I_{n-k+1}\) 而不会改变下三角结构
因此结果 \(L\) 仍为单位下三角阵；并且由于每一步选取最大的主元，由高斯变换的构造得到 \(L\) 中元素的模不会超过 \(1\) .
定理 1.2.1 若 \(A\in\mathbb{R}^{n\times n}\) ，则有排列方阵 \(P,Q\in\mathbb{R}^{n\times n}\) 使得
其中 \(L\) 中所有元素的模不大于 \(1\) ， \(U\) 中非零对角元的个数恰为 \(A\) 的秩.
虽然全主元高斯消去法弥补了不选主元的高斯消去法的不足，但仍需要付出极大的代价，因为在 \(A\) 非奇异的情况下，选主元需进行
次比较判断。于是我们考虑列主元高斯消去法
这样只进行每一列的选主元，能够大量减少时间.
代码实现
// 我们同时进行对 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}\) ，使得
称为 Cholesky 分解，其中 \(L\) 称为 \(A\) 的 Cholesky 因子.
证明 由于对称正定，因此其全部主子阵均正定，则存在 \(\widetilde{L},U\) 使得 \(A = \widetilde{L}U\) ，令
则有
注意到左边是单位上三角阵，右边是下三角阵，因此两边均为单位矩阵，于是 \(\widetilde{U} = \widetilde{L}^T,\ A=\widetilde{L}D\widetilde{L}^T\) ；由于 \(D\) 对角元为正数，则取
有 \(A = LL^T\) ，此即为满足要求的分解.
平方根法通过比较元素来计算分解，对比 \(A\) 和 \(LL^T\) 两边的对应元素得到
从左向右，每次计算一列元素，先算出对角元素，然后推导其余元素.
代码实现
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];
}
}
}
为了避免开方运算，我们使用改进的平方根法
其中 \(L\) 是单位下三角阵， \(D\) 是对角元均为正的对角阵.
仍然通过对比元素计算分解
在实际计算时，将 \(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];
}
}
}
追赶法
三对角阵的分解
对比元素得到 \(\alpha_1 = a_1\) ，而
此方法称为追赶法.
