Jacobi迭代法
求解线性方程组
A
x
=
b
\boldsymbol{Ax}=\boldsymbol{b}
Ax=b,其中
A
\boldsymbol{A}
A是
n
×
n
n\times n
n×n维可逆系数矩阵,
b
\boldsymbol{b}
b是
n
n
n维列向量。
将系数矩阵
A
\boldsymbol{A}
A分裂为
A
=
D
+
L
+
U
,
\boldsymbol{A}=\boldsymbol{D}+\boldsymbol{L}+\boldsymbol{U},
A=D+L+U, 其中,
D
=
diag
(
a
11
,
a
22
,
⋯
,
a
n
n
)
,
\boldsymbol{D}=\text {diag}(a_{11},a_{22},\cdots,a_{nn}),
D=diag(a11,a22,⋯,ann),
L
=
[
0
0
0
⋯
0
a
21
0
0
⋯
0
a
31
a
32
0
⋯
0
⋮
⋮
⋱
⋱
⋮
a
n
1
a
n
2
⋯
a
n
,
n
−
1
0
]
,
U
=
[
0
a
12
a
13
⋯
a
1
n
0
0
a
23
⋯
a
2
n
⋮
⋮
⋱
⋱
⋮
0
0
⋯
0
a
n
−
1
,
n
0
0
⋯
0
0
]
。
\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
−
1
b
,
\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)被称为迭代矩阵。上述矩阵表示可以展开成分量形式
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
≠
i
a
i
j
x
j
(
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);
}