Gauss-Seidel迭代法解线性方程组
迭代公式为:
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
i
−
1
a
i
j
x
j
(
k
+
1
)
−
∑
j
=
i
+
1
n
a
i
j
x
j
(
k
)
)
,
i
=
1
,
2
,
.
.
.
,
n
x_{i}^{(k+1)} = \frac{1}{a_{ii}}(b_{i} - \sum_{j=1}^{i-1}a_{ij}x_{j}^{(k+1)}-\sum_{j = i+1}^{n}a_{ij}x_{j}^{(k)}),i=1,2,...,n
xi(k+1)=aii1(bi−j=1∑i−1aijxj(k+1)−j=i+1∑naijxj(k)),i=1,2,...,n
function[x,err] = GaussSeidel(A,b,x0,tol,kmax)
%Input -A is n*n matrix
% -b is n*1 vector
% -x0 is initial solution row vector
% -tol is the tolerance
% -kmax is the maximum number of iteration
%Output -x is the solution
% -err is the error estimate
n = length(b);
err = 10; % 误差
k = 0; % 迭代次数
while(k < kmax)&(err>tol) % 小于最大迭代次数且误差大于容忍值时
for i = 1:n
t = b(i);
if i > 1
t = t - A(i,1:i-1)*x(1:i-1)';
end
if i < n
t = t - A(i,i+1:n)*x0(i+1:n)';
end
x(i) = t/A(i,i);
end
k = k+1;
err = max(abs(x-x0)); % 更新误差
x0 = x;
end
Tips:
与Jacobi迭代法的区别在于,使用了本次迭代已经求出的新的值来更新解,Jacobi迭代法的公式为:
x
i
(
k
+
1
)
=
1
a
i
i
(
b
i
−
∑
j
=
1
,
j
≠
i
n
a
i
j
x
j
(
k
)
)
,
i
=
1
,
2
,
.
.
.
,
n
x_{i}^{(k+1)} = \frac{1}{a_{ii}}(b_{i} - \sum_{j = 1,j\not=i}^{n}a_{ij}x_{j}^{(k)}),i=1,2,...,n
xi(k+1)=aii1(bi−j=1,j=i∑naijxj(k)),i=1,2,...,n