欢迎前往个人博客 驽马点滴 和视频空间 哔哩哔哩-《挨踢日志》

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

对中等规模的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:Nmult=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:1X(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 iterationX=L;               %initial the X=[x1;x2;…;xN]=L=[d01;d02;…;d0N]% the kth iteration of valuing the X for j=1:Nif j==1X(1)=(1-w)*X(1)+w*(B(1)-C(1)*X(2))/D(1);elseif j==NX(N)=(1-w)*X(N)+w*(B(N)-A(N-1)*X(N-1))/D(N);else%X contains the kth approximationsX(j)=(1-w)*X(j)+w*(B(j)-A(j-1)*X(j-1)-C(j)*X(j+1))/D(j);endenderr=abs(norm(X-L));  %get the error           L=X;relerr=err/(norm(X)+eps);if (err<delta)|(relerr<delta) %fit the over condition of iterationbreakendend

分析误差:

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

欢迎前往个人博客 驽马点滴 和视频空间 哔哩哔哩-《挨踢日志》

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

  1. matlab trisys,Matlab实现——严格对角占优三对角方程组求解(高斯赛尔德Gauss-Seidel迭代、超松弛) | 学步园...

    严格对角占优三对角方程组求解 对中等规模的n阶的(n<100)线性方程组,直接法的准确性和可靠性,所以常采用直接法 对于较高阶的方程组,特别是地于某些偏微分方程离散化后得到的大型稀疏方程组(系统 ...

  2. 数值计算大作业:Jacobi与Gauss -Seidel迭代求解线性方程组(Matlab实现)

    作为研究生的入门课,数值计算的大作业算是所有研究生开学的重要编程作业. 我把Jacobi与Gauss -Seidel迭代求解线性方程组的数值计算作业在MATLAB中编程实现.具体的程序详细标注后放在文 ...

  3. 高斯赛尔德、牛顿拉尔逊matlab潮流计算

    个人电气博文目录传送门 :学好电气全靠它,个人电气博文目录(持续更新中-) 本文讲解复杂电力系统的潮流计算.具体文档程序见下载链接: 文章目录 一.复杂电力系统潮流计算 二.程序 三.下载链接 一.复 ...

  4. 追赶法求解三对角方程组

    1. 来源和背景 对于一个(主)三对角方程组,我们常用"追赶法"来进行求解. 而三对角方程组常常出现于微分方程的数值求解,例如热传导方程的边值问题 {y′′(x)=f(x,y,y′ ...

  5. 运用雅可比和高斯赛德尔迭代公式求解方程组,并尝试将矩阵变为主对角占优矩阵

    程序描述 首先要求用户输入矩阵的大小n(默认不超过10),然后再提示用户输入大小为n的方阵.因为输入的方阵可能含有较多的0元素,因此用了数据结构上的矩阵的压缩方法来存储稀疏矩阵.矩阵的每一个非零元用一 ...

  6. 10种基于MATLAB的方程组求解方法

    线性方程组的求解包括直接法和迭代法,其中迭代法包括传统的高斯消元法,最速下降法,牛顿法,雅克比迭代法,共轭梯度法,以及智能启发式算法求解法和神经网络学习算法,传统算法可以相互组合改进,智能仿生启发式算 ...

  7. MATLAB数值分析学习笔记:线性代数方程组的求解和高斯消元法

    工程和科学计算的许多基本方程都是建立在守恒定律的基础之上的,比如质量守恒等,在数学上,可以建立起形如 [A]{x}={b} 的平衡方程.其中{x}表示各个分量在平衡时的取值,它们表示系统的状态或响应: ...

  8. 龙格库塔法解微分方程组的matlab程序,MATLAB实例源码教程:龙格库塔法求解微分方程组源代码实例.doc...

    MATLAB实例源码教程:龙格库塔法求解微分方程组源代码实例.doc MATLAB实例源码教程龙格库塔法求解微分方程组源代码实例题目用经典 Runge-Kutta方法求下列一阶微分方程组的近似解y1 ...

  9. Python解线性方程组的直接法(6)————求解三对角方程组的追赶法

    求解三对角方程组的追赶法 import numpy as npdef zuiganfa(A, d):n = A.shape[0]l = np.mat(np.zeros(n, dtype=float

最新文章

  1. 财务大数据比赛有python吗-【教改实验班简介】财务大数据分析班
  2. C语言程序设计 | 操作符介绍与使用方法
  3. 测量仪图片_南昌高度仪价格,大行程非标影像测量仪组装
  4. uniapp防抖操作
  5. 快来一起玩转LiteOS组件:Curl
  6. USACO_1_2_Dual Palindromes
  7. linux NFS 配置步骤
  8. keras 初步学习
  9. MySQL游标(cursor) 定义及使用
  10. 【DS3231 RTC实时时钟模块与Arduino接口构建数字时钟】
  11. 东线报接口 全网一手线报全网(京东,淘宝,天猫)最全优惠信息
  12. 华为主题包hwt下载_hwtTool下载-华为主题开发工具下载 v9.0.2.301 官方版[百度网盘资源] - 安下载...
  13. 别再用手机管家了!华为手机删除这几个文件夹,能瞬间释放大量内存
  14. Chrome打不开baidu的解决方法
  15. [译] UX 设计实践:如何设计可扫描的 Web 界面
  16. 部署测试fabric1.0及源码解析
  17. PHP对接国际验证码接口DEMO示例
  18. 借助 Material Design,帮助您打造更好的无障碍应用 (中篇)
  19. 30岁改行学python_我30岁了,转行学编程可以吗? 排除法告诉你答案
  20. 飘在上海,心路历程【2】2011年11月 反思与拷问

热门文章

  1. EMOTIV Epoc X 无线便携式脑电仪
  2. PyTorch 1.0 中文文档:torch.utils.model_zoo
  3. 最小系统板 STM32入门,点亮 LED 灯(STM32F103C6T6)
  4. nag在逆向中是什么意思_NAG在医学是什么意思
  5. IO端口、IO内存、IO空间、内存空间的含义和联系
  6. github.io 公共博客
  7. 自定义SeekBar实现实现进度提示随thum移动
  8. 常垒·视频:股权投资的终极思维
  9. 各大OJ刷题平台汇总
  10. emplace_back()