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附近的時候誤差相對較小(有點鬱悶，針對這個三對角矩陣時沒能達到加速的目的)。總之，迭代法的舍入誤差隨著迭代次數的增加，能達到相當高的精度；而且收斂速度令人滿意。