Jacobi迭代法求解线性方程组

Chenxr32

已于 2023-05-14 17:00:34 修改

阅读量1.9k

点赞数

CC 4.0 BY-SA版权

分类专栏： C/C++ 文章标签：线性代数 c语言

于 2023-05-14 16:31:11 首次发布

本文链接：https://blog.youkuaiyun.com/QDchenxr/article/details/130669361

C/C++ 专栏收录该内容

7 篇文章

订阅专栏

Jacobi迭代法是一种用于求解线性方程组Ax=b的技术，其中A是n×n的可逆矩阵。方法涉及到将矩阵A分解为D（对角部分）、L（下三角）和U（上三角）。迭代公式表示为x(k+1)=D^(-1)(b-Dx(k))，收敛条件是迭代矩阵的谱半径小于1。当A为严格对角占优或不可约对角占优时，该方法收敛。文章还给出了C语言实现Jacobi迭代法的代码示例。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Jacobi迭代法

求解线性方程组 $Ax=b\boldsymbol{Ax}=\boldsymbol{b}$ ，其中 $A\boldsymbol{A}$ 是 $n×nn\times n$ 维可逆系数矩阵， $b\boldsymbol{b}$ 是 $n$ 维列向量。
将系数矩阵 $A\boldsymbol{A}$ 分裂为 $A=D+L+U，\boldsymbol{A}=\boldsymbol{D}+\boldsymbol{L}+\boldsymbol{U}，$ 其中， $,ann)，\boldsymbol{D}=\text {diag}(a_{11},a_{22},\cdots,a_{nn})，$ $L=[000⋯0a2100⋯0a31a320⋯0⋮⋮⋱⋱⋮an1an2⋯an,n−10]，U=[0a12a13⋯a1n00a23⋯a2n⋮⋮⋱⋱⋮00⋯0an−1,n00⋯00]。\boldsymbol{L}= \begin{bmatrix} 0&0&0&\cdots&0 \\ a_{21}&0&0&\cdots&0 \\ a_{31}&a_{32}&0&\cdots&0\\ \vdots&\vdots&\ddots&\ddots&\vdots \\ a_{n1}&a_{n2}&\cdots&a_{n,n-1}&0 \end{bmatrix}， \boldsymbol{U}= \begin{bmatrix} 0&a_{12}&a_{13}&\cdots&a_{1n}\\ 0&0&a_{23}&\cdots&a_{2n} \\ \vdots&\vdots&\ddots&\ddots&\vdots \\ 0&0&\cdots&0&a_{n-1,n}\\ 0&0&\cdots&0&0 \end{bmatrix}。$ Jacobi迭代法的矩阵表示为 $x(k+1)=−D−1(L+U)x(k)+D−1b，\boldsymbol{x}^{(k+1)}=-\boldsymbol{D}^{-1}(\boldsymbol{L+U})\boldsymbol{x}^{(k)}+\boldsymbol{D}^{-1}\boldsymbol{b}，$ 其中， $−D−1(L+U)-\boldsymbol{D}^{-1}(\boldsymbol{L+U})$ 被称为迭代矩阵。上述矩阵表示可以展开成分量形式 $xi(k+1)=1aii(bi−∑j≠iaijxj(k))。x_i^{(k+1)}=\frac{1}{a_{ii}} (b_i - \sum_{j\neq i}{a_{ij}x_j^{(k)}})。$ 定理1：迭代法对任意初始向量 $x(0)\boldsymbol{x}^{(0)}$ 都收敛的充分必要条件是迭代矩阵的谱半径小于1。
定理2：若 $A\boldsymbol{A}$ 为严格对角占优，或不可约对角占优矩阵，则Jacobi迭代法收敛。

C语言实现Jacobi迭代法

void jacobi(const Matrix *A, const Matrix *b, Matrix *x, const int it)
{
	double tmp = 0;
	int r = 0, c = 0;
	double *x_tmp = (double *)malloc(sizeof(double) * x->row);
	for (int k = 0; k < it; k++)
	{
		for (r = 0; r < A->row; r++)
		{
			tmp = 0;
			for (c = 0; c < A->column; c++)
			{
				if (c != r)
					tmp += A->data[r * (A->column) + c] * x->data[c];
			}
			x_tmp[r] = (b->data[r] - tmp) / (A->data[r * (A->column) + r]);
		}
		memcpy(x->data, x_tmp, sizeof(double) * x->row);
	}
	free(x_tmp);
}