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

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

我们知道无论是有限差分方法还是有限元方法，最后都离不开解大规模方程组。除提到的这三种方法外，常用的方法还有共轭梯度法（CG）和多重网格方法以及他们的一些在不同条件下的改进，比如用$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++进行有限元和差分方法的编程工作。