用Matlab實現Gauss_Seidel迭代法解線性方程組

今天中午看見代做群有個題目，就是做一個G-S迭代，本來想接下來，可是就慢了幾分鐘就被別人搶走。不過我反正也沒事幹就把程式碼敲了。

高斯－賽德爾迭代（Gauss–Seidel method）是數值線性代數中的一個迭代法，可用來求出線性方程組解的近似值。該方法以卡爾·弗里德里希·高斯和路德維希·賽德爾命名。同雅可比法一樣，高斯－賽德爾迭代是基於矩陣分解原理。
OK簡要的摘一下Wiki的內容，接下來來看看演算法層是怎麼實現的。

演算法

對於一個含有n個未知量及n個等式的如下線性方程組：

在求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]  %輸出迭代結果

結果怎麼說呢，這個方程就不能收斂。。不過用來解其他對角佔優矩陣是沒問題了。