Jacobi迭代法
求解线性方程组 Ax=b\boldsymbol{Ax}=\boldsymbol{b}Ax=b,其中A\boldsymbol{A}A是n×nn\times nn×n维可逆系数矩阵,b\boldsymbol{b}b是nnn维列向量。
将系数矩阵A\boldsymbol{A}A分裂为 A=D+L+U,\boldsymbol{A}=\boldsymbol{D}+\boldsymbol{L}+\boldsymbol{U},A=D+L+U, 其中,D=diag(a11,a22,⋯ ,ann),\boldsymbol{D}=\text {diag}(a_{11},a_{22},\cdots,a_{nn}),D=diag(a11,a22,⋯,ann), 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}。L=0a21a31⋮an100a32⋮an2000⋱⋯⋯⋯⋯⋱an,n−1000⋮0,U=00⋮00a120⋮00a13a23⋱⋯⋯⋯⋯⋱00a1na2n⋮an−1,n0。 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},x(k+1)=−D−1(L+U)x(k)+D−1b,其中,−D−1(L+U)-\boldsymbol{D}^{-1}(\boldsymbol{L+U})−D−1(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)}})。xi(k+1)=aii1(bi−j=i∑aijxj(k))。 定理1:迭代法对任意初始向量x(0)\boldsymbol{x}^{(0)}x(0)都收敛的充分必要条件是迭代矩阵的谱半径小于1。
定理2:若A\boldsymbol{A}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);
}