Java教程

线性方程组的直接解法

本文主要是介绍线性方程组的直接解法,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!

三角形方程组和三角分解

前代法

求解下三角形方程组

\[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\) .
 
回代法通过公式

\[x_i = \left(y_i-\sum_{j=i+1}^{n}u_{ij}x_j\right)\bigg{/}u_{ii} \]

从下到上求解方程组.

 

代码实现

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 \]

此方法称为追赶法.

这篇关于线性方程组的直接解法的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!