-
数值分析:Python实现高斯赛德尔迭代法(Gauss-Seidel)与雅可比迭代法(Jacobi)
2021-01-19 13:02:39Python实现高斯赛德尔迭代法(Gauss-Seidel)与雅可比迭代法(Jacobi) 数值分析:Python实现高斯赛德尔迭代法(Gauss-Seidel)与雅可比迭代法(Jacobi) 一、基本迭代法 二、雅可比迭代法(Jacobi) 三、高斯—...Python实现高斯赛德尔迭代法(Gauss-Seidel)与雅可比迭代法(Jacobi)
数值分析:Python实现高斯赛德尔迭代法(Gauss-Seidel)与雅可比迭代法(Jacobi)
一、基本迭代法
二、雅可比迭代法(Jacobi)
三、高斯—赛德尔迭代法(Gauss-Seidel)
题目:分别使用雅可比迭代法与高斯—赛德尔迭代法求解Ax=b线性方程组,其中A为10阶希尔伯特矩阵输出解向量与迭代次数。
希尔伯特矩阵 Hilbert matrix
希尔伯特矩阵(Hilbert matrix),矩阵的一种,其元素A(i,j)=1/(i+j-1),i,j分别为其行标和列标。
代码实现:
import numpy as np n=10 def Jacobi(A,b,x0,xstar): k=0 while True: for i in range(0, n): sum = 0.0 for j in range(0, i): sum = sum + A[i][j] * x0[j] for j in range(i + 1, n): sum = sum + A[i][j] * x0[j] xstar[i] = (b[i] - sum) / A[i][i] temp = np.fabs(xstar[0] - x0[0]) for j in range(1, n): if np.fabs(xstar[j] - x0[j]) > temp: temp = np.fabs(xstar[j] - x0[j]) for j in range(0, n): x0[j] = xstar[j] k = k + 1 if (temp < 1.0e-6 or k > 1000) : break print("Jacobi:",k) def Gauss(A,b,x0,xstar): k = 0 while True: for i in range(0, n): sum = 0.0 for j in range(0, i): sum = sum + A[i][j] * xstar[j] for j in range(i + 1, n): sum = sum + A[i][j] * x0[j] xstar[i] = (b[i] - sum) / A[i][i] temp = np.fabs(xstar[0] - x0[0]) for j in range(1, n): if np.fabs(xstar[j] - x0[j]) > temp: temp = np.fabs(xstar[j] - x0[j]) for j in range(0, n): x0[j] = xstar[j] k = k + 1 if (temp < 1.0e-6 or k > 1000) : break print("Gauss:",k) def main(): A=np.zeros([n,n],float) b=np.zeros(n,float) x0=np.zeros(n,float) Jex = np.zeros(n,float) GSex = np.zeros(n,float) for i in range(0,n): for j in range(0,n): A[i][j] = 1./(i+1+j) for i in range(0,n): A[i][i] = A[i][i] + 0.5 b[i] = 1 x0[i] = 0 Jacobi(A,b,x0,Jex) for i in range(0,n): x0[i] = 0 Gauss(A,b,x0,GSex) for i in range(0,n): print("Jex[",i,"]=",Jex[i]) print("\n") for i in range(0,n): print("GSex[",i,"]=",GSex[i]) if __name__ == '__main__': main()
通过编码计算,对比两种算法的收敛速度与结果。可以看出,高斯赛德
尔迭代法的收敛速度比雅可比迭代法收敛速度要快的多。 -
如何利用MATLAB求解线型方程组--雅可比迭代法、高斯赛德尔迭代法
2020-04-24 21:26:15文章目录前言1 直接法2 迭代法小结 前言 今天我们要说的就是数值微积分,赶紧看看他和高等数学中的微积分有什么区别吧。本文是科学计算与MATLAB语言专题六第2小节的学习笔记,如果大家有时间的话,可以去听听课,...前言
今天我们要说的就是数值微积分,赶紧看看他和高等数学中的微积分有什么区别吧。本文是科学计算与MATLAB语言专题六第2小节的学习笔记,如果大家有时间的话,可以去听听课,没有的话,可以看看我的笔记,还是很不错的。
1 直接法
1.线性方程组的直接解法
高斯(Gauss)消去法
列主元消去法
矩阵的三角分解法
高斯(Gauss)消去法是一个经典的直接法,由它改进得到的列主元消去法,是目前计算机上求解线性方程组的标准算法,其特点就是通过消元将一般线性方程组的求解问题转化为三角方程组的求解问题。此外,还有矩阵的三角分解法等许多直接求解算法。
(1)利用左除运算符的直接解法
MATLAB提供了一个左除运算符“\”用于求解线性方程组,它使用列主元消去法,使用起来十分方便。对于线性方程组,可以利用左除运算符反斜杠求解,b左除以A可获得线性方程组的数值解x。
→
注意,如果矩阵A是奇异的或接近奇异的,则MATLAB会给出警告信息。
例1 用左除运算符求解下列线性方程组。
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; x=A\b
(2)利用矩阵分解求解线性方程组矩阵分解是设计算法的重要技巧,是指将一个给定的矩阵分解成若干个特殊类型矩阵的乘积,从而将一个一般的矩阵计算问题转化为几个易求的特殊矩阵的计算问题。通过矩阵分解方法求解线性方程组的优点是运算速度快,可以节省存储空间。
LU分解
QR分解
Cholesky分解
①LU分解的基本思想
矩阵的LU分解就是将一个n阶矩阵表示为一个下三角矩阵和一个上三角矩阵的乘积。线性代数中已经证明,只要方阵是非奇异的,LU分解总是可以进行的。
→
对于三角方程很容易求解,于是可以首解向量y使Ly=b,再求解Ux=y,从而达到求解线性方程组Ax=b的目的。
②MATLAB的LU分解函数
LU分解函数是根据列主元LU分解算法定义的,具有较好的数据稳定性。lu函数有两种调用格式:
[L,U]=lu(A):
产生一个上三角阵U和一个变换形式的下三角阵L,使之满足A=LU。注意,这里的矩阵A必须是方阵。
[L,U,P]=lu(A):
产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P,使之满足PA=LU。同样,矩阵A必须是方阵。
当使用第一种格式时,矩阵L往往不是一个下三角阵,但可以通过行交换成为一个下三角阵。
③用LU分解求解线性方程组
→→
或
→→→
通过LU分解后可以大大提高运算速度。
例2 用LU分解求解例1中的线性方程组。A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; [L,U]=lu(A); x=U\(L\b)
x = -66.5556 25.6667 -18.7778 26.5556
2 迭代法
2.线性方程组的迭代解法
迭代法是一种不断用变量的原值推出它的新值的过程,是用计算机解决问题的一种基本方法。
(1)雅可比(Jacobi)迭代法
求解公式为:
与之对应的迭代公式为:
)
雅可比迭代法的函数文件jacobi.m:function [y,n]=jacobi(A,b,x0,ep) D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; y=B*x0+f; n=1; while norm(y-x0)>=ep x0=y; y=B*x0+f; n=n+1; end
(2)高斯-赛德尔(Gauss-Serdel)迭代法
Gauss-Serdel迭代法的函数文件gauseidel.mfunction [y,n]=gauseidel(A,b,x0,ep) D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=(D-L)\U; f=(D-L)\b; y=B*x0+f; n=1; while norm(y-x0)>=ep x0=y; y=B*x0+f; n=n+1; end
例3 分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0,迭代精度为10^(-6)。
=A=[4,-2,-1;-2,4,3;-1,-3,3]; b=[1,5,0]'; [x,n]=jacobi(A,b,[0,0,0]',1.0e-6) [x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)
x = 0.9706 0.8529 1.1765 n = 35 x = 0.9706 0.8529 1.1765 n = 16
例4 分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组,看是否收敛。
=A=[1,2,-2;1,1,1;2,2,1]; b=[9;7;6]; [x,n]=jacobi(A,b,[0;0;0],1.0e-6) [x,n]=gauseidel(A,b,[0;0;0],1.0e-6)
x = -27 26 8 n = 4 x = NaN NaN NaN n = 1012
小结
直接法:以矩阵初等变换为基础,可以求得方程组的精确解;占用的内存空间大、程序实现较为复杂;一般适合求解低阶稠密线性方程组。
迭代法:从给定初始值逐步逼近精确解的过程,求解过程占用存储空间小,程序设计简单;适用于求解大型稀疏矩阵线性方程组;要考虑算法的收敛性。
大家学会了吗?另外分享给大家一个Markdown的公式神器- ,为了编辑出优美的公式,大家可以多看看官方支持文档.最后没有一键三连,但欢迎大家
点赞👍,收藏⭐,转发🚀,
如有问题、建议,请您在评论区留言💬哦。 -
【雅可比迭代法&高斯—赛德尔迭代法】
2019-03-17 18:33:50雅可比迭代法欢迎使用Markdown编辑器新的改变功能快捷键合理的创建标题,有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右...雅可比迭代法
雅可比迭代法-高斯-赛德尔迭代法求线性方程组
- 求解方程组
- matlab程序实现
function [x,k] = Fjacobi(A,b,x0,nm,eps) %用雅可比迭代法求解方程组Ax=b %输入:A为方程组的系数矩阵,b为方程组右端的列向量, %输入:x0为迭代初值构成的列向量,nm为最大迭代次数,eps为误差精度 %输出:x为求得的方程组的解构成的列向量,k为迭代次数 D = diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); d=D\b; k=1; disp('雅可比迭代法'); while k<=nm x=B*x0+d; %雅可比迭代式 if norm(x-x0)<eps disp('迭代次数为');k disp('方程组的解为');x return; end x0=x; k=k+1; end %下面:如果达到最大迭代次数仍不收敛,输出警告语句及迭代的最终结果(并不是方程组的解) disp('在最大迭代次数内不收敛!'); disp('最大迭代次数后的结果为');x function [x,k] = Fgseid(A,b,x0,nm,eps) %用高斯-塞德尔迭代法法求解方程组Ax=b %输入:A为方程组的系数矩阵,b为方程组右端的列向量,X为迭代初值构成的列向量,nm为最大迭代次数,eps为误差精度 %输出: x为求得的方程组的解构成的列向量,k为迭代次数 D = diag(diag(A)); L = -tril(A,-1); U = -triu(A,1); G = (D-L)\U; d = (D-L)\b; k = 1; disp('高斯-塞德尔迭代法'); while k<=nm x = G * x0 + d;%计算迭代矩阵 if norm(x-x0)<eps disp('迭代次数为');k disp('方程组的解为');x return; %上面:达到精度要求就结束程序,输出迭代次数和方程组的解 end x0 = x; k = k + 1; end %下面:如果达到最大迭代次数仍不收敛,输出警告语句及迭代的最终结果(并不是方程组的解) disp('在最大迭代次数内不收敛!'); disp('最大迭代次数后的结果为');x
- 求解函数;
a = [5 -1 -1 -1; -1 10 -1 -1; -1 -1 5 -1; -1 -1 -1 10]; b = [ 4 13 8 34]'; x0 = [0 0 0 0]'; nm = 100; eps = 0.0001; [x,k] = Fjacobi(a,b,x0,nm,eps); [x,k] = Fgseid(a,b,x0,nm,eps);
结果:
-
高斯赛德尔迭代算法 C语言
2013-06-09 13:28:11迭代法是一种逐次逼近的方法,与直接法(高斯消元法)比较, 具有: 程序简单,存储量小的优点。特别适用于求解系数矩阵为大型稀疏矩阵的方程组。常用迭代方法:雅可比迭代,高斯-赛德尔迭代,松弛迭代等。 -
高斯赛德尔迭代c语言_超松弛迭代法(SOR)
2021-01-22 20:16:13本节,我们介绍超松弛迭代法。下面我们对该算法进行收敛性分析:这里我们解释一下什么叫严格对角占优矩阵:严格对角占优矩阵就是矩阵主对角线元素的的模大于与他同行的其他元素的模的总和。下面我们具体代码实现一下...一般而言,因雅可比迭代收敛速度不够快,所以在工程中应用不多。并且在雅可比迭代收敛速
度很慢的情况下,通常高斯-赛德尔方法也不会很快。我们可以对高斯-赛德尔方法做出一定的
修改,来提高收敛速度。本节,我们介绍超松弛迭代法。
下面我们对该算法进行收敛性分析:
这里我们解释一下什么叫严格对角占优矩阵:严格对角占优矩阵就是矩阵主对角线元素的的模
大于与他同行的其他元素的模的总和。
下面我们具体代码实现一下:
%超松弛(SOR)迭代法,计算线性方程组的解 function [x,k] = SOR_iteration(A,b,x0,w,tol) % tol为输入误差容限,x0为迭代初始值 % 默认最多迭代300次,超出300次会给出警示 max = 300; if(w<=0||w>=2) % MATLAB中error语句用于报错跳出,并可以给出相应提示 error('错啦!w的值不符合要求'); % 在MATLAB中遇到return的意思就是结束整个函数的执行, % break是仅仅跳出循环体 return; end % 取出X矩阵的对角元,然后构建一个以X对角元为对角的对角矩阵 D = diag(diag(A)); % 求A的下三角矩阵,对角线元素为0,再每个矩阵元素取负号 L = -tril(A,-1); % 求A的上三角矩阵,对角线元素为0,再每个矩阵元素取负号 U = -triu(A,1); % 在MATLAB中inv是求矩阵的逆矩阵的意思,同具有一样的功能 B = inv(D-L*w)*((1-w)*D+w*U); f = w*inv((D-L*w))*b; x = B*x0+f; k = 1; while norm(x-x0)>=tol x0 = x; x = B*x0+f; k = k+1; if(k>=max) disp('迭代次数超过',max1,'次,方程组可能不收敛'); return; end [k,x'] end
我们用超松弛迭代法求解一下下面的方程组:
在编辑界面输入如下命令:
A = [4,-2,-1;-2,4,-2;-1,-2,3]; % 注意必须加'表示转置 b = [0,-2,3]'; x0 = [1,1,1]'; [x,k] = SOR_iteration(A,b,x0,1.45,1e-6)
最后计算结果如下:
-
matlab 高斯迭代代码_超松弛迭代法(SOR)
2020-11-28 11:46:13本节,我们介绍超松弛迭代法。下面我们对该算法进行收敛性分析:这里我们解释一下什么叫严格对角占优矩阵:严格对角占优矩阵就是矩阵主对角线元素的的模大于与他同行的其他元素的模的总和。下面我们具体代码实现一下... -
数值中的一些方法概述
2015-11-22 15:05:34插值法: 拉格朗日插值牛顿插值多项式三次样条插值 微积分: 牛顿-柯斯特公式复合辛普森公式 解线性方程组的迭代法:雅可比迭代法与高斯-赛德尔迭代法 解非线性方程:牛顿法 -
常用计算方法(C语言代码)(计算方法课程)
2018-04-06 10:35:12目录 矩阵相乘 ...高斯赛德尔迭代法及超松弛迭代法 拉格朗日插值 最小二乘法求解拟合曲线 牛顿插值多项式 (均差) 变步长梯形求积法 复化梯形公式 与 复化辛普森公式求积 改进的欧拉法 求解 微分方... -
计算机数值方法
2013-10-17 21:25:44使用雅可比迭代法或高斯-赛德尔迭代法对下列方程组进行求解。 实验四 矩阵特征值与特征向量问题 使用幂法求A模为最大的特征值及其相应的特征向量。 实验五 代数插值 使用拉格朗日插值法或牛顿插值法求解:已知f(x)... -
《数值分析》作者: 钟尔杰,黄廷祝 出版时间: 2004年
2019-06-07 19:00:18解线性方程组的直接法§3.1 高斯消元法§3.2 列主元消元法与三角分解§3.3 直接三角分解法§3.4 向量和矩阵范数§3.5 方程组直接方法的误差估计应用:小行星轨道问题习题三第四章 线性方程组的迭代解法§4.1 ... -
C#版数值计算算法教材及实例代码
2011-05-13 11:17:434.12 高斯—赛德尔迭代法 4.13 求解对称正定方程组的共轭梯度法 4.14 求解线性最小二乘问题的豪斯荷尔德变换法 4.]5 求解线性最小二乘问题的广义逆法 4.16 病态方程组的求解 第5章 非线性方程与方程组的求解 5.1 非... -
常用算法程序集(C语言描述) 第三版 (PDF高清电子书+附书源码打包)
2017-11-02 18:39:146.10 高斯-赛德尔迭代法 6.11 求解对称正定方程组的共轭梯度法 6.12 求解线性最小二乘问题的豪斯荷尔德变换法 6.13 求解线性最小二乘问题的广义逆法 6.14 求解病态方程组 第7章 非线性方程与方程组的求解 7.1 求非... -
MATLAB语言常用算法程序集
2012-02-22 15:58:01gauseidel 高斯-赛德尔迭代法求线性方程组Ax=b的解 SOR 超松弛迭代法求线性方程组Ax=b的解 SSOR 对称逐次超松弛迭代法求线性方程组Ax=b的解 JOR 雅可比超松弛迭代法求线性方程组Ax=b的解 twostep 两步迭代法求线性... -
MATLAB常用算法
2010-04-05 10:34:28gauseidel 高斯-赛德尔迭代法求线性方程组Ax=b的解 SOR 超松弛迭代法求线性方程组Ax=b的解 SSOR 对称逐次超松弛迭代法求线性方程组Ax=b的解 JOR 雅可比超松弛迭代法求线性方程组Ax=b的解 twostep 两步迭代法求线性... -
C#数值计算算法编程 代码
2009-12-11 11:31:584.12 高斯—赛德尔迭代法 4.13 求解对称正定方程组的共轭梯度法 4.14 求解线性最小二乘问题的豪斯荷尔德变换法 4.]5 求解线性最小二乘问题的广义逆法 4.16 病态方程组的求解 第5章 非线性方程与方程组的求解 5.1 非... -
C++常用算法程序集
2011-03-24 09:15:323.10 高斯\|赛德尔迭代法161 3.11 求解对称正定方程组的共轭梯度法165 3.12 求解线性最小二乘问题的豪斯荷尔德变换法169 3.13 求解线性最小二乘问题的广义逆法175 3.14 求解病态方程组189 第4章 非线性方程与方程组... -
C开发金典随书源码:含数据结构 数值计算分析 图形图像处理 目录和文件操作 系统调用方面的范例
2013-10-25 13:12:12范例1-32 头插法建立单链表 75 ∷相关函数:createlist函数 1.3.2 限制链表长度建立单链表 77 范例1-33 限制链表长度建立长单链表 77 ∷相关函数:createlist函数 1.3.3 尾插法建立单链表 79 范例1-34 尾插法... -
C语言通用范例开发金典.part2.rar
2012-08-31 14:18:18范例1-32 头插法建立单链表 75 ∷相关函数:createlist函数 1.3.2 限制链表长度建立单链表 77 范例1-33 限制链表长度建立长单链表 77 ∷相关函数:createlist函数 1.3.3 尾插法建立单链表 79 范例1-34 尾插法... -
C 开发金典
2013-06-20 16:20:03范例1-32 头插法建立单链表 75 ∷相关函数:createlist函数 1.3.2 限制链表长度建立单链表 77 范例1-33 限制链表长度建立长单链表 77 ∷相关函数:createlist函数 1.3.3 尾插法建立单链表 79 范例1-34 尾插法...