设线性方程组为AX=B,则Gauss-Seidel迭代法的迭代公式为
\[{x}_{m}^{k+1}=\frac{1}{{a}_{mm}}(-\sum_{i=1}^{m-1} {{a}_{mi}{x}_{i}^{k+1}-\sum_{i=m+1}^{n} {{a}_{mi}{x}_{i}^{k}}+{b}_{m}})\]
对应的矩阵表达式为
\[{x}^{\{k+1\}}={B}_{0}{x}^{k}+f\]
其中,B0为Gauss-Seidel迭代矩阵,B0=(D+L) 1U,f=(D+L)1b,D为对角矩阵,L和U分别为严格下三角矩阵和严格上三角矩阵。据此,可以编写实现Gauss-Seidel迭代法的M文件Gauss.m,如下所示。
code.matlab
function s=Gauss(a,b,x0,eps)
% 用Gauss-Seidel迭代法解线性方程组
% a为系数矩阵,b为方程组ax=b的右端项,x0为初值
if nargin==3
eps=1.0e-6;
elseif nargin<3
error
return
end
D=diag(diag(a)); %求对角矩阵
L=tril(a,-1); %求严格下三角矩阵
U=triu(a,1); %求严格上三角矩阵
C=inv(D+L);
B=-C*U;
f=C*b;
s=B*x0+f;
while norm(s-x0)>=eps
x0=s;
s=B*x0+f;
end
return
【例4】用上面编写的Gauss函数求解下列方程组
\[\left\{ \begin{matrix} 10{x}_{1}-2{x}_{2}-{x}_{3}=6 \\ -2{x}_{1}+2{x}_{2}-{x}_{3}=10 \\ -{x}_{1}-2{x}_{2}+5{x}_{3}=10 \end{matrix} \right.\]
下面用MATLAB实现,在命令窗口中输入
code.matlab
>> a=[10 -2 -1;-2 2 -1;-1 -2 5];
>> b=[6 10 10]';
>> x0=[0 0 0]';
>> gauss(a,b,x0)
ans =
4.0000
13.0000
8.0000