Jacobi迭代、Gauss_Seidel迭代和最佳因子SOR迭代的比較

zzu小陸發表於2017-04-01
IDE

Jacobi迭代、Gauss_Seidel迭代和最佳因子SOR迭代的比較

我們知道無論是有限差分方法還是有限元方法，最後都離不開解大規模方程組。除提到的這三種方法外，常用的方法還有共軛梯度法（CG）和多重網格方法以及他們的一些在不同條件下的改進，比如用LDLT

LDL^T
去改進的共軛梯度演算法。這些工作分為幾次來做。

這裡編寫了Jacobi、Gauss_Seidel和SOR三種迭代演算法，對於一些簡單的東西，直接利用的matlab內建的演算法。可以很方便地line by line地改寫成C或者C++的演算法。

  • jacobi迭代Gauss_Seidel和SOR方法沒有更多優勢，特別是在速度方面。但是良好的並行性質使得它在解方程組當中有了一席之地。 jacobi.m如下：
function [x,out] = jacobi(A,b,x0,eps, iterNum)
%輸入引數eps表示精度，iterNum表示迭代步數
%當輸入引數為五個時，eps失效，以迭代步數作為終止條件
%% 
n = length(b);
if det(A) == 0 %借用了內建函式，行列式可用迴圈求解
    error('No Unique Solution Or No Solution!');
end
%% 當主對角線上元素不全為零時，可以通過交換兩行來實現
for i = 1 : n
    j = 1;
    while A(i,i) == 0
        if A(j,i) ~= 0 && A(i,j) ~= 0
            tempV = A(j,:);
            A(j,:) = A(i,:);
            A(i,:) = tempV;
            tempS = b(j);
            b(j) = b(i);
            b(i) = tempS;
        else
            j = j + 1;
        end
        if j > n
            error('No Exchange Exist!');
            %若無論如何都無法通過交換兩行使得對角線上元素為0時，原則上可以通過交換
            %兩列實現，此時求得的結果對應元素也需交換。這裡若通過交換行不能使得主對
            %角線非0，產生錯誤資訊。
        end
    end
end
%% 針對不同的引數輸入，來選擇精度控制模式或者步數控制模式。
D = diag(diag(A));
invD = diag(1./diag(A));
J = invD * (D - A);
if max(abs(eig(J))) >= 1
    error('Jacobi Algorithm Cannot Convergence！')
end
%這裡利用了Jacobi迭代收斂的充要條件判斷是否收斂，求譜半徑內建函式。事實上，可以
%通過矩陣分解、反冪法等手段求特徵值。此判斷並非必要，只是給A的選取是否合適，作為
%一個參考。收斂與否可通過最後x的變化來判斷。
invDb = invD * b;
f = @(x)(J*x + invDb);
x = x0;
if nargin == 5
    for k = 1:iterNum-1
        x = f(x);
    end
    x_0 = x;
    x = f(x);
    out = norm(x-x_0,2);
else
    if nargin == 3
        eps = 1.0e-6;%預設精度
    end
    out = 0;
    while 1
        x_0 = x;
        x = f(x);
        out = out + 1 ;
        if norm( x - x_0 ,2 ) < eps
            %這裡以最後向量之間的距離來作為迭代精度，事實上不大精準。
            break;
        end
    end
end
end
  • Gauss_Seidel迭代其實是，jacobi迭代的一種改進，在計算x的每一分量時，有新值不用舊值，使得迭代的速度有了一個質的飛躍。Gauss_Seddel.m如下：
function [x,out] = Gauss_Seidel(A,b,x0,eps,iterNum)
%% 同Jacobi迭代，只是迭代函式f有所改變。
n = length(b);
if det(A) == 0
    error('No Unique Solution Or No Solution!');
end
for i = 1 : n
    j = 1;
    while A(i,i) == 0
        if A(j,i) ~= 0 && A(i,j) ~= 0
            tempV = A(j,:);
            A(j,:) = A(i,:);
            A(i,:) = tempV;
            tempS = b(j);
            b(j) = b(i);
            b(i) = tempS;
        else
            j = j + 1;
        end
        if j > n
            error('No Exchange Exist!');
        end
    end
end
D = diag(diag(A));
invD = diag(1./diag(A));
J = invD * (D - A);
invDb = invD * b;
H = eye(n) - tril(A)^-1*A;
if max(abs(eig(H))) >= 1
    error('Gauss_Seidel Algorithm Cannot Convergence！')
end

    function x = f(x)
        for l = 1:n
            temp_x = J(l,:)*x + invDb(l);
            x(l) = temp_x;
        end
    end
x = x0;
if nargin == 5
    for k = 1:iterNum-1
        x = f(x);
    end
    x_0 = x;
    x = f(x);
    out = norm(x-x_0,2);
else
    if nargin == 3
        eps = 1.0e-6;
    end
    out = 0;
    while 1
        x_0 = x;
        x = f(x);
        out = out + 1 ;
        if norm( x - x_0 ,2 ) < eps
            break;
        end
    end
end
end
  • SOR超鬆弛迭代方法，和Gauss_Seidel迭代相比，在有新值的時候不用舊值，但也不用新值，去二者的加權平均，當權重取得比較合適時，在三者當中是最優的。SOR.m如下：
function [x,out] = SOR(A,b,x0,w,eps,iterNum)
%% 同Gauss_Seidel迭代，只是迭代函式f有所改變。增加了鬆弛因子w作為引數。
n = length(b);
if det(A) == 0
    error('No Unique Solution Or No Solution!');
end
for i = 1 : n
    j = 1;
    while A(i,i) == 0
        if A(j,i) ~= 0 && A(i,j) ~= 0
            tempV = A(j,:);
            A(j,:) = A(i,:);
            A(i,:) = tempV;
            tempS = b(j);
            b(j) = b(i);
            b(i) = tempS;
        else
            j = j + 1;
        end
        if j > n
            error('No Exchange Exist!');
        end
    end
end
D = diag(diag(A));
invD = diag(1./diag(A));
J = invD * (D - A);
% if strcmp(w ,'best') %最佳鬆弛因子，對矩陣的要求較高。不在程式中考慮。
%     w = 2/(   1+sqrt(   1-max(abs(eig(J)))^2   )   );
% end
w_ = 1 - w;
invDb = invD * b;
L = invD * (D - tril(A));
B_w = D*(eye(n) - w*L)/w;
H_w = eye(n) - B_w^-1*A;
if max(abs(eig(H_w))) >= 1
    error('SOR Algorithm Cannot Convergence！')
end

    function x = f(x)
        for l = 1:n
            temp_x = J(l,:)*x + invDb(l);
            x(l) = w_*x(l) + w*temp_x;
        end
    end
x = x0;
if nargin == 6
    for k = 1:iterNum-1
        x = f(x);
    end
    x_0 = x;
    x = f(x);
    out = norm(x-x_0,2);
else
    if nargin == 4
        eps = 1.0e-6;
    end
    out = 0;
    while 1
        x_0 = x;
        x = f(x);
        out = out + 1 ;
        if norm( x - x_0 ,2 ) < eps
            break;
        end
    end
end
end
  • 有了這三種迭代函式，就可以做一些比較簡單的比較工作。
clc
clear
%% 報錯檢驗

%當矩陣A奇異時迭代產生錯誤（這裡只以Jacobi迭代為例）
A = [1 1 ; 2 2];
b = [1 1]';
x0 = [1 1]';
jacobi(A,b,x0)

%當不能通過行交換使對角線上元素都變成非零時，產生錯誤。
A = [1 1 1;
     1 1 0;
     0 1 0];
b = [1 1 1]';
x0 = [1 1 1]';
jacobi(A,b,x0)

%當迭代不收斂時，該迭代方法失效
A = [1 2 ; 3 4];
b = [1 1]';
x0 = [1 1]';
jacobi(A,b,x0)

%% 造一個次序相容的A，且滿足Jacobi迭代矩陣J的特徵值都為實數，J譜半徑小於1.
%在這樣的條件下，我們看看SOR方法最佳鬆弛因子實驗值與理論值是否相符。
% 實驗中，可以控制相同步數，選精度最高因子，也可以選定精度，看那個因子對應迭代
%的步數最少。實驗選到的最優因子和理論最優因子有一定差距，原因在於：1、所採用的精度
%標準，並非真正的精度，這裡的精度和解的逼近只有一個粗糙關係。2、不同因子對應的迭代
%步數可能會相同，相同步數對應的精度也有可能對相同。
clc
clear
A = [-3 1 0 1;1 -3 0 0; 0 0 -3 1; 1 0 1 -3];%A次序相容
n = size(A,2);
b = ones(n,1);
x0 = ones(n,1);

D = diag(diag(A));
J = eye(n) - D^-1*A;
if isreal(eig(J)) && (max(abs(eig(J))) < 1)
sep = 0:0.001:2;
w0 = sep(2:end-1);
m = size(w0,2);
Ts = zeros(1,m);
step_sets = zeros(1,m);
h = 0;
for w = w0 %採用不同因子用迭代法嘗試，在相同迭代步數下選取精度最高。
[~,out] = SOR(A,b,x0,w,1.0e-4,15);
h = h+1;
Ts(h) = out;
[~,out] = SOR(A,b,x0,w,1.0e-4);
step_sets(h) = out;
end
inds = find(Ts==min(Ts));
ind = round((max(inds) + min(inds))/2);
%當多個精度相同時，取中間一個。
wb = w0(ind);%w 均勻分割迭代產生的最優鬆弛因子
wb0 = 2/( 1+sqrt( 1-max(abs(eig(J)))^2 ) );% 理論最優鬆弛因子
wb_error = wb - wb0;%理論與實際差距
[~,ind1]=min(abs(w0 - wb0));%最接近理論最優因子的實驗因子的位置
step_error = ind - ind1;%位置步數差距
p = 500;
plot(w0(ind1-p:ind1+p),step_sets(ind1-p:ind1+p),'-mo',...
                'MarkerSize',3)%繪圖直觀
            title('同精度所需迭代步數隨鬆弛因子變化圖')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%這裡提供一個求次序相容矩陣的方法
% k1 = 2;
% k2 = 2;
% n = k1 + k2;
% M1 = rand(k1,k2);
% D1 = diag(rand(1,k1)) + eye(k1)*n;
% D2 = diag(rand(1,k2)) + eye(k2)*n;
% M2 = M1';
% A = [D1 M1 ; M2 D2];
% D = diag(diag(A));
% while ~isreal(eig(eye(n)-D^-1*A)) || max(abs(eig(eye(n)-D^-1*A))) >= 1
% D1 = diag(rand(1,k1)) + eye(k1)*n;
% D2 = diag(rand(1,k2)) + eye(k2)*n;
% M1 = rand(k1,k2);
% M2 = M1';
% A = [D1 M1 ; M2 D2];
% D = diag(diag(A));
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% 我們來比較Jacobi迭代、Gauss_Seidel迭代和最優因子SOR迭代的速度。由於三
%%個演算法每一步迭代所需的運算量相差不大，而程式中使用的內建函式不一而同，所以我們不
%%能通過比較時間，可通過比較步數來看看他們之間差距。程式中，我們通過不同迭代步數下
%%達到的實時精度情況，來直觀地看三種方法的收斂效能。
clc
clear
close all
%%造一個主對角佔優的正定矩陣，使之三種方法都收斂。
n = 10;
A = rand(n,n);
A = A'+A + eye(n)*2*n;
x0 = ones(n,1);
true_x = rand(n,1)*10;%真解
b = A*true_x;

% 通過上面介紹的方法，求SOR方法的最優鬆弛因子。
sep = 0:0.001:2;
w0 = sep(2:end-1);
m = size(w0,2);
Ts = zeros(1,m);
step_sets = zeros(1,m);
h = 0;
for w = w0
[~,out] = SOR(A,b,x0,w,1.0e-4,15);
h = h+1;
Ts(h) = out;
[~,out] = SOR(A,b,x0,w,1.0e-4);
step_sets(h) = out;
end
inds = find(Ts==min(Ts));
ind = round((max(inds) + min(inds))/2);
%當多個精度相同時，取中間一個。
wb = w0(ind);%w 均勻分割迭代產生的最優鬆弛因子

%不同迭代方式的一個比較
iterNum_Down = 1; %設定迭代步上下限
iterNum_Up = 12;
eps = 0.0001;
Ts_j = [];
Ts_G = [];
Ts_S = [];
Xs_j = [];
Xs_G = [];
Xs_S = [];
iterNums = iterNum_Down : iterNum_Up;
for iterNum = iterNums
[x_j,out_j] = jacobi(A,b,x0,'',iterNum);
Ts_j(end+1) = out_j;
Xs_j(:,end+1) = x_j;
[x_G,out_G] = Gauss_Seidel(A,b,x0,'',iterNum);
Ts_G(end+1) = out_G;
Xs_G(:,end+1) = x_G;
[x_S,out_S] = SOR(A,b,x0,wb,'',iterNum);
Ts_S(end+1) = out_S;
Xs_S(:,end+1) = x_S;
end
plot(iterNums,Ts_j,iterNums,Ts_G,'--mo',iterNums,Ts_S,'-.r*');
legend('jacobi迭代','Gauss迭代','最優因子SOR迭代');
title('不同迭代方式的收斂速度比較');
% 三者真正精度和matlab內建函式的一個比較。
figure;
for k = 1:iterNum_Up - iterNum_Down + 1
error_j(k)= norm((Xs_j(:,k) - true_x), 2);
error_G(k)= norm((Xs_G(:,k) - true_x), 2);
error_S(k)= norm((Xs_S(:,k) - true_x), 2);
end
error_m = norm(A\b - true_x);
plot(iterNums,error_j,'b--o',iterNums,error_G,'--c*',iterNums,error_S,'r',iterNums,error_m,'-mo');
legend('jacobi迭代精度變化','Gauss迭代精度變化','最優因子SOR迭代精度變化','matlab內建函式精度');
title('不同迭代方式的收斂精度比較');
%由於收斂速度都是呈高數量級的且不同方法收斂速度的巨大差異（特別是Jacobi迭代在圖上
%降低了解析度），所以我們很難從圖上直觀看出我們這三種方法隨著迭代的進行，細緻
%地看出收斂效能的不同。但是，我們可以從數值上簡單地看到，通過十幾步的迭代，
%三種方法的計算精度就已經達到了matlab資料型別（顯示上的）精度（e-15左右）的上限。

%去掉jacobi迭代來考察Gauss和SOR收斂速度比較。
figure;
plot(iterNums,error_G,'--c*',iterNums,error_S,'r',iterNums,error_m,'-mo');
legend('Gauss迭代精度變化','最優因子SOR迭代精度變化','matlab內建函式精度');
title('不同迭代方式的收斂精度比較');

%看看Gauss_Seidel迭代幾步就能夠達到matlab內建演算法的精度。
steps = 0;
x = x0 ;
while 1
    x = Gauss_Seidel(A,b,x,'',1);
    steps = steps + 1;
    if norm(x-true_x,2) < error_m
        break;
    end
end
%發現一個神奇的現象是，當某種演算法迭代到一定的步數之後，和真解之間就不再靠近了，總
%是保持同一個距離，這是因為，數量達到了機器內部所能表示的精度極限了（e-15），所以
%繼續做迭代，解不會發生變化，總保持和真解之間有一個（e-15）數量級的距離。
% error_G(19) - error_G(20)
% vpa(error_G(16),50)
% vpa(error_G(18),50)
% vpa(error_G(19),50)
% vpa(error_G(20),50)

這裡寫圖片描述
這裡寫圖片描述
這裡寫圖片描述
這裡寫圖片描述
程式設計過程中對於一些必要的部分儘量避開了使用MATLAB內建函式，這樣在程式改寫為C或者C++程式時變得特別方便。但是為了工作的快速和高效以及小矩陣運算的速度追求，這裡暫時先採用MATLAB編寫。

MATLAB採用了微軟的MKL函式庫，所以在矩陣運算方面特別有優勢。但是，當矩陣的規模比較大時，一個變數存不下一個矩陣，那麼此時只能另外在申請記憶體，MATLAB也就失去了它得天獨厚的MATLAB矩陣運算能力，變得沒有太大意義了。另一方面，在計算速度上，由於MATLAB採用解釋性編譯，變數記憶體也是動態分配的，在速度上也是被C等甩下了不只一條街。故後面可能會採用C或者C++進行有限元和差分方法的程式設計工作。

相關文章