Matlab实现——严格对角占优三对角方程组求解（高斯赛尔德Gauss-Seidel迭代、超松弛）

 

严格对角占优三对角方程组求解

对中等规模的n阶的(n<100)线性方程组，直接法的准确性和可靠性,所以常采用直接法

对于较高阶的方程组，特别是地于某些偏微分方程离散化后得到的大型稀疏方程组(系统矩

阵绝大多数为零元素)，由于直接解法的计算代价较高，使得迭代法更具有竞争力。

于是设计以下的2种算法：

                                       ——(1)

 其系数矩阵是对角的，且元素满足严格对角占优：

   

 

1)追赶法：

利用方程组(1)的特点，应用Gauss消元法求解时，每步只需消一个元素。其消元过程为：

 

    ——(2)

得到同解方程组(仍然严格对角占优)为：

 

                                            ——(3)

回代过程(对角占优，不必选主元)为：

 

                           ——(4)

 

用追赶法解方程组(2)仅需(5n-4)次乘除过程，(3n-3)次加减过程,算法时间复杂度O(n)。

程序：

function X=trisys(A,D,C,B)%Input- A is the subdiagonal of the coefficient matrix%     - D is the main diagonal of the coefficient matrix%     - C is the superdiagonal of the coefficient matrix %     - B is the constant vector of the linear system%Output - X is the solution vectorN=length(B);X=zeros(N,1);for k=2:N   mult=A(k-1)/D(k-1);   D(k)=D(k)-mult*C(k-1);   B(k)=B(k)-mult*B(k-1);endX(N)=B(N)/D(N);for k= N-1:-1:1   X(k)=(B(k)-C(k)*X(k+1))/D(k);end

这个方法的精度很高，和系统的内置函数linsolve(H,B)的求解结果一致

 

 

 

2)迭代法(采用改进的Gauss-Seidel迭代)（这个方法是看了超松弛迭代后，得出的类似方法）：

原理介绍： 

 

                        ——(1)

 

它的Gauss-Seidel迭代方法我们已经很熟悉了：

【1】      给一个初始列向量：

 

 

【2】利用迭代公式：

 

 

经过一定的迭代次数以后，就能得到近似解：

                

现在对上述的Gauss-Seidel迭代进行加速

得到迭代方法：

(注意：当ω=1时，就是我们所熟悉的Gauss-Seidel迭代)

其中：

ω是迭代加速的相关系数——松弛因子

上述方法可解释为第k+1次迭代近似解的各分量依次为用Gauss-Seidel方法求得的第k+1次迭代近似值和第次近似值的加权平均值。适当选取收敛因子ω(事实上叫做松弛因子)，可望该方法比Gauss-Seidel迭代法收敛得更快。

根据以上的原理分析，作出程序如下：

function X=acc(A,D,C,B,P,delta, max1,w)%Input- A is the subdiagonal of the coefficient matrix%     - D is the main diagonal of the coefficient matrix%     - C is the superdiagonal of the coefficient matrix %     - B is the constant vector of the linear system%     - P is an N x 1 matrix; the initial guess%     - w is the convergence multiplicate%     - delta is the tolerance for P%     - max1 is the maximum number of iterations% Output - X is an N x 1 matrix: the gauss-seidel approximation%           to the solution of AX = BN = length(B);L=P;                   %L is a mediutfor k=1:max1          %max1th iteration    X=L;               %initial the X=[x1;x2;…;xN]=L=[d01;d02;…;d0N]    % the kth iteration of valuing the X     for j=1:N      if j==1       X(1)=(1-w)*X(1)+w*(B(1)-C(1)*X(2))/D(1);      elseif j==N         X(N)=(1-w)*X(N)+w*(B(N)-A(N-1)*X(N-1))/D(N);      else        %X contains the kth approximations        X(j)=(1-w)*X(j)+w*(B(j)-A(j-1)*X(j-1)-C(j)*X(j+1))/D(j);      end   end   err=abs(norm(X-L));  %get the error              L=X;   relerr=err/(norm(X)+eps);   if (err<delta)|(relerr<delta) %fit the over condition of iteration       break   endend

 

分析误差：

这时候我们看到w=0.2时误差是e=0.01550147497154

 

经过类似的试验可以知道：

在w=0.98-w=1之间的时候存在最优松弛因子

 

我们看到迭代次数的增加，带来误差的显著减小，且迭代次数max1=20的时候精度达到1.0e-11

可见，该方法的求解精度还是令人满意的。

         在求解该问题的过程中，对于求解方程组的方法选择是一个很重要的因素，注意到这个系数矩阵是50阶严格对角占优三对角稀疏矩阵，查询了相关知识后，我个人认为，50阶的严格对角占优三对角稀疏矩阵，完全可以用高斯消去法，这是因为高斯消去后的上三角(或者下三角)仍然是严格对角占优，而对于这个稀疏矩阵，迭代法是一个非常不错的选择，而我采取的迭代法受限制的就是这个松弛因子w，注意到0<w≤1的时候，该方法是任何初始向量P都收敛，于是采取了w=0：0.2：1的选择方式，最后发现w=1附近的时候误差相对较小(有点郁闷，针对这个三对角矩阵时没能达到加速的目的)。总之，迭代法的舍入误差随着迭代次数的增加，能达到相当高的精度；而且收敛速度令人满意。