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

原创 2011年09月13日 18:54:49
 

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

对中等规模的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 vector

N=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);

end

X(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 = B

N = length(B);

L=P;                   %L is a mediut

for 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

   end

end


 

分析误差:

这时候我们看到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附近的时候误差相对较小(有点郁闷,针对这个三对角矩阵时没能达到加速的目的)。总之,迭代法的舍入误差随着迭代次数的增加,能达到相当高的精度;而且收敛速度令人满意。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

相关文章推荐

噪声图像的超分辨重建

简介   本篇主要是对论文:Super-Resolving Noisy Images 的简单笔记记录。 该论文主要是亮点:对噪声图像做超分辨重建,能同时达到去噪和部分细节恢复效果。 实现原理   ...

投影矩阵与最小二乘(三)

其实这篇文章作为投影矩阵与最小二乘的完结篇,已经不完全是投影矩阵与最小二乘的关系了,更多的是在投影矩阵的基础上发展出来的一些理论 先说一下标准正交矩阵的概念: 对于一个矩阵A,如果A的列向量是标准...

非严格主对角矩阵的转化(matlab实现)

Question: Design/implement a method to transform a matrix into strict diagonal dominance, if possib...

Matlab实现——Jacobi Iteration and Gauss-Seidel Iteration

Use Jacobi Iteration and Gauss-Seidel Iteration to solve the following linear system with several di...

数值分析 Gauss-Seidel迭代法求解线性方程组 MATLAB程序实现

Gauss-Seidel迭代法 参考数值分析第四版 颜庆津著 P39 运行输入 运行结果 函数内容(保存为gauss.m文件,在MATLAB中运行) %Gauss-Seidel迭代法求解线性方程组。迭...

分别用雅可比(Jacobi)迭代法和高斯—塞德尔(Gauss—Seidel)迭代法求解线性方程组

算法介绍(迭代法介绍): 代码C语言实现; # include # include # define N 6 /* *使用雅可比迭代法和高斯-赛德尔迭代法 求解线性方程组 ...

高斯赛德尔迭代求解矩阵的解

  • 2016年01月20日 16:43
  • 121KB
  • 下载

MATLAB追赶法求解三对角方程

  • 2017年03月28日 22:13
  • 350B
  • 下载

Matlab的Gauss_Seidel迭代方法解线性方程组

用Matlab实现Gauss_Seidel迭代法解线性方程组今天中午看见代做群有个题目,就是做一个G-S迭代,本来想接下来,可是就慢了几分钟就被别人抢走。不过我反正打了一天的LOL累了,也没事干就把代...

带状对角矩阵的LU分解及回代求解算法实现

算法名称:带状对角矩阵的LU分解及回代求解   算法描述:        分解主要是使用笔者前面几篇文章提到过的Crout方法。因为不可能把一个带状对角矩阵A的LU分解也像其压缩形式本是一样紧凑...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:Matlab实现——严格对角占优三对角方程组求解(高斯赛尔德Gauss-Seidel迭代、超松弛)
举报原因:
原因补充:

(最多只允许输入30个字)