Matlab的Gauss_Seidel迭代方法解线性方程组

用Matlab实现Gauss_Seidel迭代法解线性方程组

今天中午看见代做群有个题目，就是做一个G-S迭代，本来想接下来，可是就慢了几分钟就被别人抢走。不过我反正也没事干就把代码敲了。

高斯－赛德尔迭代（Gauss–Seidel method）是数值线性代数中的一个迭代法，可用来求出线性方程组解的近似值。该方法以卡尔·弗里德里希·高斯和路德维希·赛德尔命名。同雅可比法一样，高斯－赛德尔迭代是基于矩阵分解原理。
OK简要的摘一下Wiki的内容，接下来来看看算法层是怎么实现的。

算法

对于一个含有n个未知量及n个等式的如下线性方程组：

$$\begin{aligned}a_{11}\cdot x_{1}+a_{12}\cdot x_{2}+\ldots +a_{1n}\cdot x_{n}&=b_{1},\\a_{21}\cdot x_{1}+a_{22}\cdot x_{2}+\ldots +a_{2n}\cdot x_{n}&=b_{2},\\\vdots \qquad \qquad \qquad &=\vdots \\a_{n1}\cdot x_{1}+a_{n2}\cdot x_{2}+\ldots +a_{nn}\cdot x_{n}&=b_{n}.\end{aligned}$$
高斯-赛德尔迭代法的基本步骤就是：
用k用来计量迭代步数，给定该方程组解的一个近似值 x0，然后对于每个元素可以使用如下公式进行迭代：
$$\begin{aligned}{\displaystyle x_{m}^{k+1}={\frac {1}{a_{mm}}}\left(b_{m}-\sum _{j=1}^{m-1}a_{mj}\cdot x_{j}^{k+1}-\sum _{j=m+1}^{n}a_{mj}\cdot x_{j}^{k}\right),\quad 1\leq m\leq n.} \end{aligned}$$

在求k+1步近似值时，我们利用第m个方程求解第m个未知量。在求解过程中，所有已解出的k+1步元素都被直接使用，重复上述的求解过程，可以得到一个线性方程组解的近似值数列，设置误差值error，如果方程最终有解，则该方法最终会收敛在一个稳定值。

举个小例子，可以看图：

注意：为了保证该方法可以进行，主对角线元素需非零。

一个DEMO

我之前还没有接触过这个数值算法，但是看到这题目，我抑制不住我内心的洪荒之力。于是顺手编了一个小家伙，看看效果怎么样。代码如下：

prompt= 'please input the  Coefficient matrix: ';
A=input(prompt);%%[m,n]=size(A);
prompt= ' please input the  Column Vector: ';
B=input(prompt);
B=B';
prompt= ' please input the  number of iterations: ';
k=input(prompt);
x1=1;x2=1;x3=1;
   for i=1:1:k
   x1=(B(1,1)-(A(1,2)*x2)-(A(1,3)*x3))/A(1,1);
   x2=(B(2,1)-(A(2,1)*x1)-(A(2,3)*x3))/A(2,2);
   x3=(B(3,1)-(A(3,2)*x2)-(A(3,1)*x1))/A(3,3);
   end
x1,x2,x3

迭代20次，结果出乎我的意料啊。
x1 =-1.9727e+27
x2 = 2.4828e+28
x3 =-1.6801e+28
这这这，是什么情况？28次方，你要上天啊？

干脆搞个一般程序

于是我怀恨在心，不怀好意地阅读了别人的代码，然后编了一个一般方程组都能用的G-S迭代方法程序，采用设置距离误差的方法停止迭代，试验结果还不错。

clear;clc
format compact
%% 输入线性方程组矩阵，A为系数矩阵，C为常数列向量。
A = [2 5 1 ;
     5 0.5 4 ;
     7 8 11 ];% 相关系数矩阵
C = [2;6;1;];% 常数列向量
n=length(A);
X = ones(n,1);
Error_eval = ones(n,1);

%% 检查矩阵是否是完全严格对角占优矩阵，如果不是则无解。
for i = 1:n
    j = 1:n;
    j(i) = [];
    B = abs(A(i,j));
    Check(i) = abs(A(i,i)) - sum(B);
    if Check(i) < 0
        fprintf('矩阵在第%3i行上不是严格的对角线\n\n',i)
    end
end

%% 迭代过程如下

k = 0;
while max(Error_eval) > 0.001
    k = k + 1;
    Z = X;  
    for i = 1:n
        j = 1:n; 
        j(i) = []; 
        Xtemp = X;  
        Xtemp(i) = []; 
        X(i) = (C(i) - sum(A(i,j) * Xtemp)) / A(i,i);
    end
    Xsolution(:,k) = X;
    Error_eval = sqrt((X - Z).^2);
end

%% 结果
GaussSeidelTable = [1:k;Xsolution]'   %输出迭代过程
Solution = [A X C]  %输出迭代结果

结果怎么说呢，这个方程就不能收敛。。不过用来解其他对角占优矩阵是没问题了。

---