迭代精化下求解三对角Toeplitz线性方程组的快速算法
Matlab追赶法和迭代法解线性方程组
Matlab追赶法和迭代法解线性⽅程组实验⽬的:1)追赶法解三对⾓阵;2)掌握解线性⽅程组的迭代法;3)⽤Matlab实现Jacobi及超松弛迭代法实验要求:1)掌握追赶法解三对⾓阵2)掌握解线性⽅程组的迭代法3)提交追赶法、Jacobi及超松弛迭代法的m⽂件实验内容:1)追赶法解三对⾓矩阵⽅程(m⽂件)习题1. ⽤追赶法的m⽂件求解2)Jacobi迭代法解线性⽅程组(m⽂件)对不同初值⽤Jacobi迭代法解习题1并⽐较结果。
3)超松弛迭代法解线性⽅程组(m⽂件)对不同松弛因⼦解习题1并⽐较结果。
实验步骤: 代码:1 %追赶法2 %输⼊:系数矩阵A和因变量d;3 %输出:⾃变量x4 function z=zuigan(A,d)5 n=length(d);6 %取三对⾓元素a,b,c7for i=1:n-18 a(i)=A(i,i);9 b(i)=A(i+1,i);10 c(i)=A(i,i+1);11 end12 a(n)=A(n,n);13 %分解系数矩阵14 u(1)=a(1);15 l(1)=c(1)/a(1);16for i=2:n-117 u(i)=a(i)-b(i-1)*l(i-1);18 l(i)=c(i)/u(i);19 end20 u(n)=a(n)-c(n-1)*l(n-1);21 %解y22 y(1)=d(1)/u(1);23for k=2:n24 y(k)=d(k)-c(k-1)*y(k-1)/u(k);25 end26 %解x27 x(n)=y(n);28for k=n-1:-1:129 x(k)=y(k)-l(k)*x(k+1);30 end31 z=x;32 endzuigan 运⾏: 所得结果,较为粗糙。
代码:1 %雅克⽐迭代法2 %输⼊系数矩阵A,因变量b,初始向量x0,容许误差eps,最⼤迭代次数t3 %输出⾃变量x和迭代数n4 function [z,k]=jacobi(A,b,x0,e,t)5 %默认eps和最⼤迭代次数m6if nargin==37 e=1e-6;8 m=200;9 elseif nargin<310 error('输⼊的参数不⾜');11return;12 elseif nargin==513 m=t;14 end15 n=length(b);16 x(1,:)=x0;17 z(1,:)=x0;18for k=2:m19 sum=0;20for i=1:n21 w=0;22 u=0;23for j=i+1:n24 w=w+A(i,j)*x(k-1,j);25 end26for j=1:i-127 u=u+A(i,j)*x(k-1,j);28 end29 x(k,i)=(-1/A(i,i))*(u+w-b(i));30if sum<abs(x(k,i)-x(k-1,i))31 sum=abs(x(k,i)-x(k-1,i));32 end33 end34if sum<e35 z(k,:)=x(k,:);36return;37 end38 z(k,:)=x(k,:);39 end40 endjacobi 运⾏⽰例,初始向量x0=[0 0 0 0 0 0];和初始向量x0=[1 1 1 1 1 1]; 初始值不同,迭代次数可能不同。
线性方程组的直接法和迭代法
线性方程组的直接法直接法就是经过有限步算术运算,无需迭代可直接求得方程组精确解的方法。
线性方程组迭代法迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法.该方法具有对计算机的存贮单元需求少,程序设计简单、原始系数矩阵在计算过程中不变等优点,是求解大型稀疏矩阵方程组的重要方法.迭代法不是用有限步运算求精确解,而是通过迭代产生近似解逼近精确解.如Jacobi 迭代、Gauss — Seidel 迭代、SOR 迭代法等。
1. 线性方程组的直接法直接法就是经过有限步算术运算,无需迭代可直接求得方程组精确解的方法。
1.1 Cramer 法则Cramer 法则用于判断具有n 个未知数的n 个线性方程的方程组解的情况。
当方程组的系数行列式不等于零时,方程组有解且解唯一。
如果方程组无解或者有两个不同的解时,则系数行列式必为零。
如果齐次线性方程组的系数行列式不等于零,则没有非零解。
如果齐次线性方程组有非零解,则系数行列式必为零。
定理1如果方程组Ax b =中0D A =≠,则Ax b =有解,且解事唯一的,解为1212,,...,n n D D Dx x x D D D===i D 是D 中第i 列换成向量b 所得的行列式。
Cramer 法则解n 元方程组有两个前提条件: 1、未知数的个数等于方程的个数。
2、系数行列式不等于零 例1 a 取何值时,线性方程组12312312311x x x a ax x x x x ax ++=⎧⎪++=⎨⎪++=⎩有唯一解。
解:211111111011(1)11001A a a a a a a ==--=--- 所以当1a ≠时,方程组有唯一解。
定理2当齐次线性方程组0Ax =,0A ≠时该方程组有唯一的零解。
定理3 齐次线性方程组0Ax =有非零解0A <=>=。
1.2 Gauss 消元法Gauss 消元法是线性代数中的一个算法,可用来为线性方程组求解,求出矩阵的秩,以及求出可逆方阵的逆矩阵。
三对角方程组的追赶法
2013-2014(1)专业课程实践论文题目:三对角方程组的追赶法一、算法理论在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对角线方程组11112222211111n n n n n n n n n b c x f a b c x f a b c x f a b x f -----⎛⎫⎛⎫⎛⎫⎪⎪ ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪= ⎪⎪⎪ ⎪⎪⎪ ⎪⎪ ⎪⎝⎭⎝⎭⎝⎭, 简记为Ax f =。
求解Ax f =:等价于解两个三角形方程组,Ly f y =求;,Ux y x =求.从而得到解三对角线方程组的追赶法公式:(1)计算{}i β的递推公式()111/,/,2,3,,1;i i i i i c b c b a i n βββ==-=- (2) 解Ly f =()()11111/,/,2,3,,;i i i i i i i y f b y f a y b a i n β--==--=(3)解Ux y =1,,1,2,2,1.n n i i i i x y x y x i n n β+==-=--我们将计算系数12112n n y y y βββ-→→→→→→ 及的过程称为追的过程, 将计算方程组的解11n n x x x -→→→ 的过程称为赶的过程。
#include <stdio.h>#include <math.h>#include<stdlib.h>#define N 20double a[N], b[N], c[N-1], f[N], r[N];int n;int i;void LUDecompose(); // LU分解void backSubs(); // 回代void main(){printf("请输入方程的维数n=");scanf("%d",&n);getchar();if(n>N||n<=0){printf("由于该维数过于犀利, 导致程序退出!");return;}printf("\n输入下三角元素\n");printf("输入%d个a值: ", n-1);for (i=1; i<n; i++)scanf("%lf", &a[i]);getchar();printf("\n输入主对角线元素\n");printf("输入%d个b值: ", n);for (i=0; i<n; i++)scanf("%lf", &b[i]);getchar();printf("\n输入上三角元素\n");printf("输入%d个c值: ", n-1);for (i=0; i<n-1; i++)scanf("%lf", &c[i]);getchar();printf("\n输入%d个方程组右端项: \n", n);for (i=0; i<n; i++)scanf("%lf", &f[i]);getchar();LUDecompose();backSubs();printf("\n线性方程组的解为: \n");for (i=0; i<n; i++)printf("x%d=%lf\n", i+1, f[i]);}void LUDecompose(){ //α被b取代, β被c取代, 以节省存储空间c[0]=c[0]/b[0];for(i=1;i<n-1;i++){r[i]=a[i];b[i]=b[i]-r[i]*c[i-1];c[i]=c[i]/b[i];}r[i]=a[i];b[i]=b[i]-r[i]*c[i-1];}void backSubs(){ // y被f取代, x也被f取代, 以节省存储空间f[0]=f[0]/b[0];for(i=1; i<n; i++)f[i]=(f[i]-r[i]*f[i-1])/b[i];f[n-1]=f[n-1];for(i=n-2;i>=0;i--)f[i]=f[i]-c[i]*f[i+1];}四、 算法实现例1.用该程序计算三对角线方程组2100012100A 012100012100012⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭---=-----, 10000b ⎛⎫⎪ ⎪⎪= ⎪ ⎪ ⎪⎝⎭计算其方程组的解。
关于Toeplitz矩阵的计算_1_5
摘要本文的研究涉及三个方面的内容:Toeplitz矩阵的逆的求法与分解式、特殊Toeplitz矩阵的逆的求法与分解式、Toeplitz线性方程组的解法。
全文共分四章三个部分,创新成果着重体现在第二章、第三章及第四章。
第一部分(第一章)介绍有关Toeplitz矩阵和特殊Toeplitz矩阵的定义以及一些的简单性质。
第二部分(第二章和第三章)主要研究Toeplitz矩阵的逆的求法与分解式,特殊Toeplitz矩阵的逆的求法。
在第二章,本文先给出了Toeplitz矩阵的逆的求法和Toeplitz矩阵的逆的分解式。
本文得到一个新的结论,Toeplitz矩阵的逆矩阵显式表示为循环矩阵与下三角Toeplitz矩阵的乘积之和,并讨论了此分解式的稳定性。
在第三章,本文先给出循环Toeplitz矩阵和三对角Toeplitz矩阵的逆的求法,并且给出了显式逆。
还对五对角Toeplitz矩阵的逆进行了研究,给出了新的结论,得到五对角Toeplitz矩阵的逆的求法,而且显式地表示了五对角Toeplitz矩阵的逆。
第三部分(第四章)主要讨论了Toeplitz线性方程组的解法,介绍了用Zohar 算法、Akaike算法、Bareiss变换法、Gohberg-Kailath-Koltracht算法以及快速傅立叶方法来求解一般Toeplitz线性方程组。
同时,本文给出了新的算法来求解五对角Toeplitz线性方程组和循环五对角Toeplitz线性方程组。
关键词:Toeplitz矩阵,三对角矩阵,循环Toeplitz矩阵,五对角矩阵,逆IABSTRACTThis thesis presents a systematic research on Toeplitz matrices such as computing the inversion of Toeplitz matrices and solving the Toeplitz linear systems. The thesis consists three parts with four chapers.In part one (chapter one), we give the definitions of the Toeplitz matrix and the special Toeplitz matrix. The simple properties of the Toeplitz matrix are presented.The algorithms and the expressions for the inversion of Toeplitz matrix and special Toeplitz matrix are given in part two (chapter two and three). In chapter two, we introduce the methods computing the inversion of Toeplitz matrix and the factorization for the inversion of Toeplitz matrix. A new conclusion is obtained: the inversion of a Toeplitz matrix can be denoted as a sum of products of circulant matrices and upper triangular Toeplitz matrices. In chapter three, we give a new fast algorithm to computing the inversion of the five-diagonal Toeplitz matrix. The inversion of the tridiagonal Toeplitz matrix is also considered.In the last part (chapter four), the algorithms for solving the Toeplitz linear systems, circulant Toeplitz linear systems and the band Toeplitz linear systems are presented. At first, we introduce some classical methods for solving the Toeplitz linear equations. Then, two algorithms for solving the five-diagonal Toeplitz matrix linear equations are given.Keywords: Toeplitz matrix, tridiagonal matrix, circulant Toeplitz matrix, five-diagonal matrix, inversionII。
三对角toeplitz矩阵 python 迭代法 -回复
三对角toeplitz矩阵python 迭代法-回复三对角Toeplitz矩阵是一种特殊类型的方阵,其中除了对角线和相邻的两个对角线上的元素外,其余的元素都为零。
这种矩阵在数学、物理学和工程学等领域中具有重要的应用。
本文将介绍如何利用Python语言中的迭代法来解决三对角Toeplitz矩阵的问题。
在开始之前,我们首先需要了解一下三对角Toeplitz矩阵的定义和特点。
一个n\times n的三对角Toeplitz矩阵可以表示为:\[\begin{bmatrix}b_1 & c_1 & 0 & \ldots & \ldots & \ldots & \ldots & 0 \\a_1 & b_2 & c_2 & 0 & \ldots & \ldots & \ldots & 0 \\0 & a_2 & b_3 & c_3 & 0 & \ldots & \ldots & 0 \\\ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \\0 & \ldots & \ldots & \ldots & \ldots & a_{n-2} & b_{n-1} &c_{n-1} \\0 & \ldots & \ldots & \ldots & \ldots & 0 & a_{n-1} & b_n \\\end{bmatrix}\]其中,a_i, b_i, c_i是已知的实数。
特别地,当a_i = c_i时,该矩阵称为对称三对角Toeplitz矩阵。
MATLAB-追赶法求解三对角方程组的算法原理例题与程序
3)三对角形线性方程组123456789104100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014x x x x x x x x x x -⎡⎤⎡⎤⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎢⎥--⎢⎢⎥⎢⎢⎥-⎣⎦⎣⎦7513261214455⎡⎤⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥-⎢⎥⎢⎥⎢⎥-⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥-⎣⎦*(2,1,3,0,1,2,3,0,1,1)T x =--- 二、数学原理设系数矩阵为三对角矩阵112223311100000000000000000n n n n n b c a b c a b A a b c a b ---⎛⎫ ⎪⎪ ⎪=⎪ ⎪ ⎪⎪⎪⎝⎭L L L M M MM M M L L则方程组Ax=f 称为三对角方程组。
设矩阵A 非奇异,A 有Crout 分解A=LU ,其中L 为下三角矩阵,U 为单位上三角矩阵,记11222331100001000000010000000100,00000000000001n n nn b L U γαβγββγβ--⎛⎫⎛⎫⎪ ⎪⎪ ⎪ ⎪ ⎪∂==⎪⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪⎪ ⎪∂⎝⎭⎝⎭L L L LL L M M M MM M M MM M L L LLL可先依次求出L ,U 中的元素后,令Ux=y ,先求解下三角方程组Ly=f 得出y ,再求解上三角方程组Ux=y 。
事实上,求解三对角方程组的2追赶法将矩阵三角分解的计算与求解两个三角方程组的计算放在一起,使算法更为紧凑。
其计算公式为:1111,1111,111,2,3,,,1,2,,1ii i i i i i i ii i i i i n ni i i i c f b y i n c a b a f y y x y i n n x y x βγββαβγγβαβγ--+⎧===⎪⎪=⎪⎪⎪==-=⎪⎪⎨-⎪=⎪⎪=⎪⎪=--⎪=-⎪⎩L L 对对(*)三、程序设计function x=chase(a,b,c,f)%求解线性方程组Ax=f,其中A 是三对角阵 %a 是矩阵A 的下对角线元素a(1)=0 %b 是矩阵A 的对角线元素%c 是矩阵A 的上对角线元素c(n)=0 %f 是方程组的右端向量 n=length(f);x=zeros(1,n);y=zeros(1,n); d=zeros(1,n);u= zeros(1,n); %预处理 d(1)=b(1); for i=1:n-1 u(i)=c(i)/d(i);d(i+1)=b(i+1)-a(i+1)*u(i); end%追的过程y(1)=f(1)/d(1); for i=2:ny(i)=(f(i)-a(i)*y(i-1))/d(i); end%赶的过程 x(n)=y(n); for i=n-1:-1:1x(i)=y(i)-u(i)*x(i+1); end>> a=[0,-1,-1,-1,-1,-1,-1,-1,-1,-1];>> b=[4,4,4,4,4,4,4,4,4,4];>> c=[-1,-1,-1,-1,-1,-1,-1,-1,-1,0];>> f=[7,5,-13,2,6,-12,14,-4,5,-5];>> x=chase(a,b,c,f)x =2.00001.0000-3.00000.00001.0000-2.00003.0000-0.00001.0000-1.0000四、结果分析和讨论追赶法求解的结果为x=(2,1,-3,0,1,-2,3,0,1,-1)T。
最新三对角线性方程组的解法知识讲解
Jacobi迭代法
1.迭代过程: 对于Ax b ,设A可逆,且对于aii0,i1,2, ,n 令D为A的对角元素部分,则 A(AD )D, 继而对于方程组有 D x(D A )xb,所以有 解 ,令 xD 1(D A )xD 1b B JD 1(D A )D 1(LU ) 为迭代矩阵,fJ D1b。则Jacobi迭代法的迭 代格式为: 。 x(k1) Bjx(k) fJ
二、实验内容
考虑线性方程组:
A xb , A R n n,b R n
其中A为三对角矩阵,编制程序求解该方程组
三、实验要求
(1)考虑不同的数值解法,对方程组进行编程,编写其通用代码; (2)对不同的程序代码用实例进行验证,取矩阵:
2 1 0 0 0 1
1
2
1
0
0
0
A 0 1 2 1 0 ,b 1
0
0 1 2 1
0
0 0 0 1 2 0
x 方程初值 :
(0,0,0,0,0,0)T
0
。
(3)比较不同方法对求解三对角线性方程组的优异性。
四、实验过程
解法
追赶法
Jacobi迭代法
Gauss迭代法
SOR迭代法
(一)追赶法
• 编程思想:
矩阵的三对角线性方程组是由矩阵的直接三角分解法来推到的 来的,求解பைடு நூலகம்价于解两个三角形方程组:
B w(D w) 1 L (1 (w )D w)U
其中 x(k1) Bwx(k)fw 迭代法
且 fw wb 。当w=1时,SOR法就是Gauss
• SOR方法是Gauss迭代法的一种修正,主要思想是:设已知 x (k )
(1)Ly f ,求y;(2)Ux y ,求x
Toeplitz系数矩阵方程组的迭代解法
Toeplitz系数矩阵方程组的迭代解法方雅敏【摘要】主要讨论系数矩阵为非对称正定的Toeplitz的迭代求解,运用以系数矩阵的一个对称、反对称分裂为基础的SSS迭代方法.特别地分裂是一个中心对称分裂,可以利用中心对称矩阵的可约性来减少计算量和存储量.再通过几个数值例子验证了此方法的有效性.【期刊名称】《丽水学院学报》【年(卷),期】2008(030)002【总页数】4页(P25-28)【关键词】Toeplitz矩阵;迭代方法;矩阵分裂;谱半径;中心对称矩阵【作者】方雅敏【作者单位】丽水学院,数理学院,浙江,丽水,323000【正文语种】中文【中图分类】O241.80 引言考虑线性方程组Ax=b(1)的迭代求解,其中A∈Rn×n是非对称正定的Toeplitz矩阵,b∈Rn已知。
构造求解(1)的迭代方法需要系数矩阵A的有效的矩阵分袭。
例如Jacobi’s方法[1]和Gauss-seidel迭代方法[1],在这几种方法中A=D-L-U,其中D为对角矩阵,-L 为A的严格下三角部分,-U为A的严格上三角部分。
在广义的CG方法[2]和广义的Lanczos方法[3]中矩阵A的分裂为A=H+S,其中是矩阵A的共轭转置,而SSS迭代方法用的矩阵分裂是A=F-G,(2)这里表示矩阵A的转置。
定义1[1] (1)方阵A=M-N称为A的一个分裂,如果|M|≠0。
(2)如方阵的各对角线上的元素分别都相等,则称此矩阵为Toeplitz矩阵。
定义2[2] 方阵A如满足JAJ=A,称其为中心对称矩阵;若JAJ=-A,则称其为反中心对称矩阵。
这里J仅次对角线元为1的其余元为0的方阵。
在分裂(2)中,由Toeplitz矩阵、中心对称矩阵及反中心对称矩阵的定义可知:F 是一个对称的Toeplitz矩阵,即其为一个中心对称矩阵;G是一个反对称的Toeplitz矩阵,则其为一个反中心对称矩阵。
对中心对称矩阵和反中心对称矩阵的可约性质,为了简单起见只考虑n=2m的情况。
3线性方程组的迭代解法
三、逐次超松弛法(SOR方法)
逐次超松弛法(Successive Over Relaxation Method)可 看成是Gauss-Seidel方法的加速,Seidel迭代法是SOR方法的 特例。将Seidel方法的迭代公式
改写为
x(k1) i
1 aii
(bi
i 1
a x(k 1) ij j
k
0
1
2
3
4
5
6
x1
0
2.5000 2.9773 3.0098 2.9998 2.9999 3.0000
x2
0
2.0909 2.0289 1.9968 1.9997 2.0001 2.0000
x3
0
1.2273 1.0041 0.9959 1.0002 1.0001 1.0000
可见Gauss-Seidel迭代法比Jacobi迭代法收敛要快一些。
x(k 1) BJ x(k ) f J
0
其中
a21
a22
BJ D1(L U )
an1 ann
a12 a11 0
an2 ann
a13 a11
a23 a22
7
a1n1 a11
a2n1 a22
ann1 ann
a1n
a11
a2 n
a22 , fJ D1b
0
二、 Gauss-Seidel 迭代法
x(k ) i
xi(k )
x(k ) i
1 aii
bi
i 1
a x(k 1) ij j
j 1
n
aij
x(jk
)
j i
为加快收敛,在增量 xi(k ) 前加一个因子
三对角toeplitz矩阵 python 迭代法
三对角toeplitz矩阵python 迭代法关于三对角Toeplitz矩阵的迭代法,我们需要先了解什么是三对角矩阵和Toeplitz矩阵,然后介绍迭代法的基本原理以及如何应用于解三对角Toeplitz矩阵。
接下来,我们将详细讨论这个算法的步骤,并提供一个Python代码实现示例。
一、什么是三对角Toeplitz矩阵和近似解?1. 三对角矩阵:三对角矩阵是指除了对角线上的元素外,只有三个带状元素不为零的矩阵。
具体来说,如果一个矩阵A的元素aij满足以下条件之一,则该矩阵为三对角矩阵:(1) 当i-j >1时,元素aij为零;(2) 当i-j =1时,元素aij不为零;2. Toeplitz矩阵:Toeplitz矩阵是指具有下列形式的矩阵:A = [c0, c1, c2, ..., cd,b0, c0, c1, ..., cd,b0, b1, c0, ..., cd,...b0, b1, ..., bn-2, c0];其中,c0,c1,c2,...,cd,b0,b1,...,bn-2均为已知常数。
3. 解的近似:对于三对角Toeplitz矩阵A,我们想要找到一个向量x,使得Ax=b,其中b是一个已知向量。
由于矩阵A是特殊结构的,我们可以使用迭代法来近似求解x。
迭代法的基本思想是通过不断迭代更新当前解的值,直到满足一定的收敛条件。
二、迭代法的基本原理及应用于解三对角Toeplitz矩阵的步骤1. 基本原理:迭代法是一种通过多次迭代逼近解的方法。
在解三对角Toeplitz矩阵的问题中,我们可以使用一种叫做Thomas算法的迭代法。
该算法基于高斯消元的思想,通过将矩阵A转化为上三角矩阵,然后反向求解得到解向量x。
2. 迭代法步骤:针对解三对角Toeplitz矩阵的问题,我们可以按照以下步骤进行迭代法的求解:(1) 初始化:令a、b、c、d分别表示矩阵A的上、中、下对角线和常数向量b 的元素。
初始化变量n为矩阵A的维度。
关于Toeplitz矩阵的计算_21_25
15
第三章 特殊 Toeplitz 矩阵的逆矩阵的研究
特殊的 Toeplitz 矩阵有广泛的应用前景,特别在生物学、物理学、数学、社会 科学中的许多问题都和 Toeplitz 矩阵的理论有着密切的关系。在数理统计、石油、 地震物探及其它应用科学中,尤其是图象、数字信号处理中,常会遇到循环矩阵 这类特殊矩阵。带状 Toeplitz 方程组广泛应用于科学计算和工程计算中,尤其是在 双曲偏微分方程的数值求解中,随着科学技术的不断发展,问题的规模越来越大, 适应于大规模计算的并行算法也被陆续提出来。因此,对于这类特殊的 Toeplitz 矩 阵的研究是有必要的。
分析:经过简单计算可知,
λ k = ∑ ξ j ω kj , ( k = 1,2, " , n − 1 ) 。
j =0
n −1
通过这个定理,可以知道循环 Toeplitz 矩阵的逆矩阵仍然为循环 Toeplitz 矩阵,而 且
−1 −1 1 −1 C −1 = Fdiag (λ 0 , λ1 , " λ− , n −1 ) F
(4).循环 Toeplitz 矩阵 C = ∑ ξ j K J ;
j =0
n −1
(5). CK = KC 。 定理 3.1.1 矩阵 C 为 n 阶循环 Toeplitz 矩阵的充要条件是存在数 λ 0 , λ1 , " λ n −1
使
F −1CF = diag (λ 0 , λ1 , " λ n −1 ) 。
j =0
n −1
如果第一步和第三步借助于离散富氏变换,则整个算法的计算量为
17
O(n log n) 。
下面简单地讨论一下 r-循环 Toeplitz 矩阵的逆矩阵。 设形如 n 阶矩阵称为 r-循环 Toeplitz 矩阵为 ξ 1 " ξ n −1 ⎤ ⎡ ξ0 ⎢ rξ ξ 0 " ξ n−2 ⎥ ⎥ C (r ) = ⎢ n −1 ⎥ ⎢ # # % ⎥ ⎢ ⎣ r ξ 1 rξ 2 " ξ 0 ⎦ 特别地,当 r = 1 时,即得循环 Toeplitz 矩阵;当 r = −1 时,称之为斜循环 Toeplitz 矩阵(或反循环 Toeplitz 矩阵) ;当 r = 0 时,即得上三角 Toeplitz 矩阵。由 r-循环
解三对交线方程组的追赶法
VS
矩阵元素的微小变化
在三对交线方程组中,矩阵元素的微小变 化可能会导致解的巨大变化。这种敏感性 使得追赶法在面对某些问题时表现出数值 不稳定性。
提高数值稳定性和减小误差方法
选择合适的算法参数
在追赶法中,可以通过选择合适的算法参数来提高数值稳定性。例如,可以采用部分选主元策略来避免矩阵元素的微 小变化对解的影响。
优缺点分析
优点
追赶法具有计算量小、存储量低、易于编程实现等优点。对于大规模的三对角 线性方程组,追赶法通常比其他方法更加高效。
缺点
追赶法的适用范围有限,仅适用于系数矩阵为三对角矩阵的线性方程组。此外, 当系数矩阵不满足对角占优等条件时,追赶法可能无法收敛或收敛速度较慢。
Part
02
三对交线方程组数学模型建立
问题描述与定义
三对交线方程组
在二维平面上,给定三对直线,每对直线相交于一个点,这三对交线构成的方程组称为三对交线方程 组。
求解目标
通过给定的三对交线信息,求解出这三对直线的交点坐标。
数学模型构建方法
直线方程表示
在二维平面上,一条直线可以用一般式方程 $Ax + By + C = 0$ 表示,其中 $A, B$ 不同时为0。
THANKS
感谢您的观看
回代过程
从最后一个方程开始,依次将已知量代入方程求 解,得到未知量的值。此过程称为回代过程。
关键算法实现技巧
存储优化
追赶法中的系数矩阵是三对角 的,因此可以采用一维数组进 行存储,节省存储空间。
消元技巧
在消元过程中,需要注意消元 顺序和消元系数的选择,以确 保消元过程的稳定性和效率。
回代技巧
在回代过程中,需要按照正 确的顺序将已知量代入方程 求解,避免计算错误。
三对角toeplitz矩阵 python 迭代法 -回复
三对角toeplitz矩阵python 迭代法-回复什么是三对角Toeplitz矩阵?三对角Toeplitz矩阵是一种特殊的方阵,它具有特定的元素分布规律。
其每一行和每一列的元素满足同一个常数平移关系。
具体地说,给定一个n \times n的矩阵A,当且仅当A_{ij} = a_{i-j}时,矩阵A为三对角Toeplitz矩阵,其中a_{-n+1}, a_{-n+2}, ..., a_{0}, ..., a_{n-2}, a_{n-1}是已知的实数。
三对角Toeplitz矩阵具有诸多特性,例如可以通过简单的存储方式来节省空间、可以通过显式的公式计算逆等。
这些特性使得三对角Toeplitz矩阵在科学计算和工程应用中具有广泛的应用。
如何使用迭代法求解三对角Toeplitz矩阵?迭代法是一种逐步逼近解的方法,它通过反复迭代更新近似解,直到达到预设的收敛条件。
对于求解三对角Toeplitz矩阵的问题,我们可以借助迭代法来逐步逼近其解。
下面将给出一种常用的递推迭代方法,即Thomas算法。
首先,我们可以将三对角Toeplitz矩阵表示为:Ax = d其中,A是n \times n的三对角Toeplitz矩阵,x是未知向量,d是已知向量。
然后,我们可以将矩阵A分解为:A = LDU其中,L是下三角矩阵,D是对角矩阵,U是上三角矩阵。
根据上述分解,我们可以将原方程组重写为:LDUx = d我们可以利用前向替代和后向替代的方法来迭代求解y = Ux和z = Dy。
具体地,我们可以按照以下步骤进行迭代:1. 初始化:将原方程组中的A分解为L,D和U;2. 前向替代:计算y,满足L y = d;3. 后向替代:计算z,满足D z = y;4. 后向替代:计算x,满足U x = z;5. 返回结果:得到方程组的解x。
需要注意的是,为了保证迭代的收敛性,我们需要对原矩阵做一些预处理,例如通过对角线优化等。
总结三对角Toeplitz矩阵是一种具有特定元素分布规律的方阵,它在科学计算和工程应用中广泛使用。
数学中的线性方程组求解方法
数学中的线性方程组求解方法线性方程组是数学中的基础概念之一,广泛应用于各个领域。
求解线性方程组是数学中的一项重要任务,而线性方程组的求解方法也是数学中的研究热点之一。
本文将介绍几种常见的线性方程组求解方法,并探讨它们的优缺点。
一、高斯消元法高斯消元法是一种常见且简单的线性方程组求解方法。
其基本思想是通过初等行变换将线性方程组转化为简化的行阶梯形式,从而得到方程组的解。
高斯消元法的优点是简单易懂,适用于各种规模的线性方程组。
然而,高斯消元法的缺点也是显而易见的,当方程组的规模较大时,计算量会变得非常庞大,从而导致求解时间过长。
二、LU分解法LU分解法是一种常用的线性方程组求解方法,它将系数矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积,从而简化了方程组的求解过程。
LU分解法的优点是可以将求解过程分解为两个独立的步骤,从而提高了计算效率。
然而,LU分解法的缺点是需要消耗大量的存储空间,尤其是在方程组规模较大时,存储开销会变得非常庞大。
三、迭代法迭代法是一种逐步逼近的线性方程组求解方法,其基本思想是通过不断迭代计算,逐渐逼近方程组的解。
迭代法的优点是可以在较短的时间内得到一个近似解,并且可以通过控制迭代次数来控制解的精度。
迭代法的缺点是收敛性不易保证,有时候可能会产生发散的情况。
此外,迭代法的计算过程也比较复杂,需要进行大量的计算。
四、矩阵分解法矩阵分解法是一种将线性方程组转化为矩阵的乘法形式进行求解的方法。
常见的矩阵分解方法有QR分解、Cholesky分解等。
矩阵分解法的优点是可以将原始的线性方程组转化为一系列简单的矩阵乘法运算,从而提高了计算效率。
然而,矩阵分解法的缺点是需要对系数矩阵进行分解,这一过程可能会比较复杂且耗时。
五、特殊结构法对于一些特殊结构的线性方程组,可以使用特殊结构法进行求解。
例如,对于对称正定矩阵的线性方程组,可以使用共轭梯度法进行求解;对于稀疏矩阵的线性方程组,可以使用稀疏矩阵的存储和计算技巧进行求解。
数学中的线性方程组求解方法探究
数学中的线性方程组求解方法探究线性方程组是数学中一个非常基础且重要的概念,广泛应用于各个领域。
在实际问题中,我们常常需要求解线性方程组,以得到问题的解析解或者数值解。
本文将探究数学中的线性方程组求解方法,包括直接法和迭代法两种主要方法。
一、直接法直接法是求解线性方程组的一种常用方法,其思想是通过一系列数学运算,将线性方程组转化为一个简化形式的方程组,从而得到方程组的解。
直接法的优点是求解过程简单明了,计算量相对较小。
下面我们将介绍两种常用的直接法:高斯消元法和LU分解法。
高斯消元法是一种基本的直接法,它通过一系列的行变换将线性方程组转化为阶梯形或行最简形,从而得到方程组的解。
具体步骤如下:1. 将线性方程组写成增广矩阵的形式;2. 选取主元,即选取一个非零元素作为当前行的主元;3. 通过行变换,将主元所在列的其他元素变为零;4. 重复步骤2和步骤3,直到得到阶梯形或行最简形。
高斯消元法的缺点是在某些情况下可能会出现主元为零的情况,导致无法继续进行下去。
为了解决这个问题,可以使用列主元高斯消元法,即每次选取主元时选择当前列中绝对值最大的元素作为主元。
LU分解法是另一种常用的直接法,它将线性方程组的系数矩阵分解为一个下三角矩阵L和一个上三角矩阵U的乘积,即A=LU。
通过LU分解,可以将线性方程组的求解问题转化为两个较为简单的方程组的求解问题。
具体步骤如下:1. 对系数矩阵A进行LU分解,得到L和U;2. 将原方程组Ax=b转化为LUx=b;3. 先解Ly=b,得到向量y;4. 再解Ux=y,得到向量x,即为方程组的解。
LU分解法的优点是可以在多次求解同一个系数矩阵的线性方程组时节省计算量,因为LU分解只需要进行一次,之后只需要解两个三角方程组即可。
二、迭代法迭代法是求解线性方程组的另一种常用方法,其思想是通过迭代的方式逐步逼近方程组的解。
迭代法的优点是可以求解大规模的线性方程组,但其收敛速度较慢。
下面我们将介绍两种常用的迭代法:雅可比迭代法和高斯-赛德尔迭代法。
蒙伟-拟Toeplitz带状矩阵线性方程组的并行算法(最终版)
咸阳师范学院2010届本科毕业论文(设计)拟Toeplitz带状矩阵线性方程组的并行算法蒙伟(咸阳师范学院数学与信息科学学院陕西咸阳712000 )摘要本文主要研究了Toeplitz矩阵线性方程组的并行算法。
首先介绍了Toeplitz矩阵的定义、分类及半正定性等,然后分别给出了三对角Toeplitz线性方程组、拟对称七对角Toeplitz线性方程组的并行算法的误差分析。
最后,用Matlab语言编写了算法程序,通过实例计算,结果表明本文算法是可行的,并且验证了理论分析与实际计算的结果是一致的。
关键词:Toeplitz矩阵;Toeplitz线性方程组;误差分析;并行算法;并行性分析;拟对称七对角Toeplitz线性方程组拟Toeplitz带状矩阵线性方程组的并行算法Parallel Algorithm for Solving Banded quasi-ToeplitzMatrix Linear EquationsMeng wei(School of Mathematics and Information Science of Xianyang NormalUniversity Xianyang 712000)AbstractThis paper studies parallel algorithm of the Toeplitz matrix linear equations. Firstly, I introduce the definition, classification, and semi-positive definite of Toeplitz matrix and so on, and then give the error analysis of parallel algorithm for solving the three-diagonal Toeplitz matrix linear equations and quasi-symmetry seven diagonal Toeplitz matrix linear equations. Finally, I program the algorithm procedure using Matlab software. The results show that the algorithm is feasible, and the theoretical analysis is consistent with the result of practical calculation.KEY WORDS:Toeplitz matrix; Toeplitz matrix linear equations; Error analysis; Parellel algorithm; Quasi-symmetry seven diagonal Toeplitz matrix linear equations目录摘要 (I)Abstract (II)前言 (1)1 Toeplitz矩阵及其性质 (3)1.1 Toeplitz的定义 (3)1.2 Toeplitz的半正定性 (4)2 三对角Toeplitz矩阵线性方程组的并行算法 (5)2.1 并行算法 (5)2.2 误差分析 (8)2.3 并行性分析 (9)2.4 算例 (10)3 拟对称七对角Toeplitz矩阵线性方程组的并行算法 (11)3.1并行算法 (11)3.2 误差分析 (17)4 结论 (17)参考文献 (19)附录 (20)谢辞 (25)咸阳师范学院2010届本科毕业论文(设计)前 言Toeplitz 矩阵是数学和工程中应用广泛的特殊矩阵之一。
求解Toeplitz线性系统的快速CSCS分裂方法
求解Toeplitz线性系统的快速CSCS分裂方法史红芳【摘要】基于循环和反循环分裂迭代法(CSCS),提出一种系数矩阵为Toeplitz矩阵线性系统的新的分裂迭代法,该方法在一定的条件下收敛于Toeplitz线性系统的唯一解.数值实验表明新方法是有效的和可行的.【期刊名称】《太原师范学院学报(自然科学版)》【年(卷),期】2016(015)004【总页数】5页(P4-8)【关键词】Toeplitz矩阵;分裂;迭代方法;循环矩阵;反循环矩阵【作者】史红芳【作者单位】太原师范学院数学系,山西晋中030619【正文语种】中文【中图分类】O241考虑如下大型稀疏线性系统式中,A∈Cn×n是一个Toeplitz矩阵,且b,x∈Cn.Toeplitz线性系统被广泛应用于诸多领域,尤其在信号处理和控制理论中,见文献[1].求解线性方程有直接法和迭代法.文献[2]中介绍了求解Toeplitz线性系统的一些直接法.如Zohar算法,Akaike算法,Bareiss变换法等.在[3]中,借助于Toeplitz矩阵的特殊结构,Ng Michael K给出了一种基于循环和反循环的分裂迭代法(CSCS).在该方法中,收敛因子的上界仅依赖于循环矩阵和反循环矩阵的谱.在文献[4]的基础之上,一些学者对非Hermitian正定线性系统的分裂方式做了改进[5],该方法明显优于原先的方法.为了方便表达,这里给出一些简单的记号.通常情况下,用Cn×n表示n×n复矩阵的集合,Cn表示n维复向量空间.X*代表矩阵或向量X的共轭转置.对于复数x,Re(x)和Im(x)分别表示x的实部和虚部.ρ(A)表示矩阵A的谱半径.circ(a0,a1,…,an-1)表示以a0,a1,…,an-1为第一行元素的循环矩阵,icirc(a0,a1,…,an-1)表示以a0,a1,…,an-1为第一行元素的反循环矩阵.定义1([2])一个n×n矩阵A称为是一个Toeplitz矩阵,如果其满足ai,j=ai-j,A的对角线上的元素相同且平行于主对角线的各对角线上的元素也彼此相等.两个n阶Toeplitz矩阵相加或数乘Toeplitz矩阵,其结果仍是Toeplitz矩阵. 定义2([2])形如的n阶方阵称为r-循环矩阵.特别地,当r=1时,称之为循环矩阵,当r=-1时,称之为斜循环矩阵(或反循环矩阵).循环矩阵可以被离散傅立叶矩阵对角化,反循环矩阵可以被尺度傅里叶矩阵对角化[6].因此,我们得到这里ΛC和ΛS分别是包含C和S特征值的对角矩阵,循环矩阵和反循环矩阵的精确解可以使用快速傅里叶变换求解(FFTs).定理1C为形如(3)中r=1时的循环矩阵,且那么C的特征值为f(ωk),k=0,1,…,n-1,其中ωk为方程xn-1=0,0<k<n的解,并且满足Toeplitz矩阵A的循环反循环分裂为A=C+S:这里C是一个循环矩阵,S为反循环矩阵.CSCS迭代法给定一个初始值x(0),对于k=0,1,2,…直到x(k)收敛,计算这里α是一个正的常数.理论分析表明,如果C和S都是正定的,那么CSCS迭代法收敛于线性系统(1)的唯一解.本文其余部分安排如下.第二节介绍新的分裂迭代法,第三节讨论新算法的收敛性,最后,通过具体的数值例子证明它的有效性.借助于线性系统(1)系数矩阵A∈Cn×n的特殊结构,在本节中我们介绍一种新的迭代算法.本节中,我们讨论迭代矩阵的谱半径性质,并分析前面章节中算法1的收敛性.在上一节中,我们在理论上证明了新算法的收敛性,本节通过具体的数值例子来证明该算法的收敛性,同时从谱半径ρ、迭代次数(IT)和运行时间(CPU,单位s)这三个维度上与原来的CSCS方法做对比.下面的测试中,在MATLAB上执行且从零向量开始,迭代取Toeplitz方程的右端项为一个随机的生成向量,迭代误差ε=10-7,CSCS迭代的最优参数[3]如下:这里γmin和γmax分别为矩阵C和S的特征值的实部的上界和下界,βmin和βmax分别为矩阵C和S的特征值虚部的绝对值的上界和下界.例1 考虑线性系统(1),矩阵A由以下生成函数给出:aj=(1+|j|)-1,j≥0,aj=i(1+|j|)-1,j<0.并设D=d1*diag(A),这里取d1=2.8.数值实验结果见表1.由表1,我们可以看出,与CSCS迭代方法相比,新的分裂迭代算法1)在谱半径和迭代次数上没有原来的方法好,但是,随着矩阵维数的逐渐增大,这种差异逐渐减少.2)在程序运行时间上占明显的优势,并且随着矩阵维数的逐渐增大,这种差异更大. [1]CHAN R H,NG M K.Conjugate Gradient Methods for Toeplitz Systems[J].Siam Review,1996,38(3):427-482[2]徐仲,张凯院,陆全.Toeplitz矩阵类的快速算法[M].西安:西北工业大学出版社,1999:29-37[3]NG MK.Circulant and Skew-circulant Splitting Methods for Toeplitz Systems[J].Journal of Computational and Applied Mathematics,2003,159(1):101-108[4]BAI Z Z,GOLUB G H,Ng M K.Hermitian and Skew-Hermitian Splitting Methods for Non Positive Definite Linear Systems[J].Siam J Matrix Anal Appl,2001,24(3):603-626[5]WANG C L,WEN R P,BAI Y H.A New Splitting and Preconditioner for Iteratively Solving Non-Hermitian Positive Definite Systems[J].Comput.Math.Appl,2013,65:1047-1058[6]P J.Davis.Circulant Matrices[M].New York:Wiley,1979[7]ZHANG F Z.Matrix Theory[M].NewYork:Springer,2011:138-140 [8]Ng M K.Iterative Methods for Toeplitz Systems[M].New York:Oxford University Press,2004[9]夏长富.矩阵正定性的进一步推广[J].数学研究与评论,1988,8(4):499-504[10]金小庆.Toeplitz系统预处理方法[M].北京:高等教育出版社,2013:60-64[11]徐成贤,徐宗本.矩阵分析[M].西安:西北工业大学出版社,1991。
线性方程组的四种数值解法
线性方程组的四种数值解法(电子科技大学物理电子学院,四川 成都 610054)摘要:本文介绍了四种求解线性方程组的数值解法: 雅克比迭代法、高斯赛德尔迭代法、高斯消去法和改进的平方根法的基本原理和算法流程,通过求解具体方程,对四种求解方法进行了对比。
对于雅克比迭代法和高斯赛德尔迭代法,研究了两种算法对求解同一方程组的迭代效率差异,结果表明高斯赛德尔迭代法达到同样精度所需迭代次数较少。
对于高斯消去法,通过选择列主元的方法提高算法的准确度,计算结果表明高斯消去法计算精确,且运算复杂度也不是很高。
对于改进的平方根法,其运算复杂度低,但对于给定的方程组有着严苛的要求。
关键词:雅克比迭代法;高斯赛德尔迭代法;高斯消去法;改进的平方根法;线性方程组引言线性方程组的求解在日常生活和科研中有着极其重要的应用,但在实际运算中,当矩阵的维数较高时,用初等方法求解的计算复杂度随维数的增长非常快,因此,用数值方法求解线性方程组的重要性便显现出来。
经典的求解线性方程组的方法一般分为两类:直接法和迭代法。
前者例如高斯消去法,改进的平方根法等,后者的例子包括雅克比迭代法,高斯赛德尔迭代法等。
这些方法的计算复杂度在可以接受的范围内,因此被广泛采用。
一般来说,直接法对于阶数比较低的方程组比较有效;而后者对于比较大的方程组更有效。
在实际计算中,几十万甚至几百万个未知数的方程组并不少见。
在这些情况下,迭代法有无可比拟的优势。
另外,使用迭代法可以根据不同的精度要求选择终止时间,因此比较灵活。
在问题特别大的时候,计算机内存可能无法容纳被操作的矩阵,这给直接法带来很大的挑战。
而对于迭代法,则可以将矩阵的某一部分读入内存进行操作,然后再操作另外部分。
本文使用上述四种算法求解对应的方程组,验证各种算法的精确度和计算速度。
1 算法介绍1.1 雅克比迭代法 1.1.1 算法理论设线性方程组(1)b Ax的系数矩阵A 可逆且主对角元素 均不为零,令并将A 分解成 (2)从而(1)可写成令其中. (3)以B 1为迭代矩阵的迭代法(公式)(4)称为雅克比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为(5)其中为初始向量.1.1.2 算法描述 1给定迭代初始向量X 0以及误差要求delta 2根据雅克比迭代公式计算出下一组向量 3判断X 是否满足误差要求,即||X k+1 – X k || < delta4若误差满足要求,则停止迭代返回结果;若否,则返回第二步进行下一轮迭代1.2 高斯赛德尔迭代法nna ,...,a ,a 2211()nna ,...,a ,a diag D 2211=()D D A A +-=()b x A D Dx +-=11f x B x +=b D f ,A D I B 1111--=-=()()111f x B x k k +=+⎩⎨⎧[],...,,k ,n ,...,i x a ba xnij j )k (j j i iii)k (i21021111==∑-=≠=+()()()()()Tn x ,...x ,x x 002010=1.2.1 算法理论由雅克比迭代公式可知,在迭代的每一步计算过程中是用的全部分量来计算的所有分量,显然在计算第i 个分量时,已经计算出的最新分量没有被利用,从直观上看,最新计算出的分量可能比旧的分量要好些.因此,对这些最新计算出来的第次近似的分量加以利用,就得到所谓解方程组的高斯—塞德尔(Gauss-Seidel )迭代法.把矩阵A 分解成(6)其中,分别为的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成即其中(7)以为迭代矩阵构成的迭代法(公式)(8)称为高斯—塞德尔迭代法(公式),用变量表示的形式为(9)1.2.2 算法描述 1给定迭代初始向量X 0以及误差要求delta2根据高斯赛德尔迭代公式计算出下一组向量()k x ()1+k x ()1+k ix ()()1111+-+k i k x ,...,x 1+k()1+k x()1+k jx U L D A --=()nna ,...,a ,a diag D 2211=U ,L --A ()b Ux x L D +=-22f x B x +=()()b L D f ,U L D B 1212---=-=2B ()()221f x B x k k +=+⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++3判断X 是否满足误差要求,即||X k+1 – X k || < delta4若误差满足要求,则停止迭代返回结果;若否,则返回第二步进行下一轮迭代1.3 高斯消去法 1.3.1 算法理论下面三种变换称为初等行变换:1.对调两行;2.以数k ≠0乘某一行中的所有元素;3.把某一行所有元素的k 倍加到另一行对应的元素上去。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
DOI: 10.12677/pm.2020.105052
427
李姗 等
(1.2) (1.3) (1.4) (1.5) (1.6) (1.7)
理论数学
李姗 等
我们注意到,为了得到(1.6)的解 x1 和 xn ,我们需要解(1.7)中的两个线性方程组,其中 A11 是一个上三 对角矩阵。显然,向量 u 和 v 可以通过向后代入法求解。此外,由于 A11 是对角占优的,计算出的向量 u 和 v 都是稳定可靠的。
文章引用: 李姗, 刘仲云, 张育林. 迭代精化下求解三对角 Toeplitz 线性方程组的快速算法[J]. 理论数学, 2020, 10(5): 425-432. DOI: 10.12677/pm.2020.105052
李姗 等
摘要
本文主要讨论如何对三对角Toeplitz线性方程组 Ax = b 进行高精度数值求解。由于系数矩阵A这种比较 特殊的结构,使得我们可以设计出快速求解 Ax = b 的直接算法。我们将该算法应用到实际例子的计算过 程中,发现大部分例子计算效果显著,但部分例子的计算精度还不能达到计算机机器精度。针对这类达 不到计算机机器精度的例子,本文将在快速求解三对角Toeplitz线性方程组 Ax = b 的直接算法基础上, 进一步进行迭代精化,从而提高这类例子的计算精度。数值实验表明通过迭代精化,我们算法计算精度 可以达到计算机机器精度[1]。
Received: Apr. 16th, 2020; accepted: May 5th, 2020; published: May 12th, 2020
Abstract
This paper mainly discusses how to numerically solve tridiagonal Toeplitz linear systems Ax = b efficiently. Since the coefficient matrix A has a special structure, we can design a direct algorithm to quickly solve Ax = b . We will apply the above algorithm to the calculation of practical examples and find that the calculation precision of some examples is not as high as that of computer's machine accuracy. In order to improve the precision of algorithm, this paper further carries out iterative refinement to quickly solve the tridiagonal Toeplitz linear equations of Ax = b . Numerical experiments show that the computational accuracy of our algorithm can reach computer's machine accuracy by iterative refinement [1].
Open Access
1. 引言
本文考虑迭代精化下求解三对角 Toeplitz 线性方程组
Ax = b,
(1.1)
α γ
β A=
γ
β
α
其中 A = Tritoep (β ,α ,γ ) 为三对角 Toeplitz 矩阵。
关于三对角 Toeplitz 线性方程组的求解方法主要分为两类:一类是直接法如高斯消元法、循环约简 法和特殊 LU 分解法[2] [3] [4]等等。另一类就是迭代法如古典迭代法(Smith 迭代法,ADI 迭代法[5] [6] [7] [8]等)和投影迭代法(Krylov 子空间法等)。在不考虑舍入误差的情况下,对于小型稀疏求解线性方程组的 问题,我们通常使用直接法求解。因为直接法可以通过改善置换策略,引入尽量少的填充,充分利用硬 件的性能设计算法。但当遇到大型稀疏求解线性方程组的问题时,直接解法计算量太大且数值不稳定, 而迭代解法数值稳定且易于并行计算,所以这时我们通常使用迭代法求解。
我们首先从系数矩阵为次对角占优情况开始考虑。
DOI: 10.12677/pm.2020.105052
426
理论数学
0 1
01
C=
0
1
1
0
我们容易得到 Aˆ = CA 并且具有如下 2× 2 的块结构
β α γ
=Aˆ
β
α
β
γ α
≡
A11
w
T
p
0
.
α γ
0
我们对 Aˆ 作如下的 2× 2 的块 LU 分解
4:输出
x
=
xT 1
xn
T
。
( ) 算法 2 的稳定性取决于第三步的求解= xn wTv − b1 wTu 。如果 wTu 不是足够小,那么 xn 的就是相
当精确的。因此,我们可以得出结论,我们求解这类线性方程组的算法在数值上是稳定的,计算出的解 是可靠的,如果 A 是下次对角占优的并且 wTu 不是足够小。
关键词
Toeplitz矩阵,迭代精化,快速算法
Copyright © 2020 by author(s) and Hans Publishers Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). /licenses/by/4.0/
=x
= xxn1 , bˆ
b2
b1
.
用块 LU 分解法分解(1.3),我们就可以通过求下面方程的解来获得方程(1.1)的解
A11 x1 + xn p = b2 ,
w
T
x1
=
b1.
( ) = xn wTv − b1 wTu ,
x1= v − xnu.
为了获得 x = x1T xn T 的解,我们首先求解如下方程中的 u 和 v
线性方程组(1.1)转化为以下新的线性方程组
A x = b
(1.9)
x = Jx , b = Jb 。因此,结合算法 2 得到下面的算法 3 用于求解方程(1.1)。
算法 3 一个求解三对角 Toeplitz 线性方程组 Ax = b 的算法 1:输入 β ,α,γ , b ; 2:调用算法 2 去得到方程(1.9)的解 x ; 3:输出 x = Jx 。
A11u A11v
= =
p, b2 .
通过如下算法,我们得到式(1.7)中方程组的解 u 和 v ,具体算法如下:
算法 1 向后代入法求解 A11z = q
1:输入 β ,α ,γ , q ;
2:计算= α1
αβ= , γ 1
γ ;
β
3:计算 zn−=1 α1; zn−=2 γ1 − α z1 n−1 ; = 4: i 2,, n − 2 计算 zn−i−1 = −α1zn−i γ − z1 n−i+1 ;
A Fast Algorithm for Solving Tridiagonal Toeplitz Linear Systems with Iterative Refinement
Shan Li1, Zhongyun Liu1, Yulin Zhang2
1School of Mathematics and Statistics, Changsha University of Science and Technology, Changsha Hunan 2Centro de Matemática, Universidade do Minho, Braga, Portugal
算法 4 一个求解三对角 Toeplitz 线性方程组 Ax = b 的算法
1:输入 β ,α ,γ , b (若 β > γ )或 b (若 β < γ ); 2:调用算法 2 (若 β > γ )或调用算法 3 (若 β < γ )去求解方程(1.1)的解 x ; 3:输出 x 。
Pure Mathematics 理论数学, 2020, 10(5), 425-432 Published Online May 2020 in Hans. /journal/pm https:///10.12677/pm.2020.105052
Aˆ
=
I
w
T
A1−11
A11
1
p −wT A1−11 p .
−wT A1−11 p 是叫做 Aˆ 的 Schur 补。 同样地,要求问题(1.1)的解也就变成了求如下线性方程组的解
Aˆ x = bˆ
其中=bˆ C=b (b2 ,b3 ,,bn ,b1 )T 。对 x 和 bˆ 做如下格式的划分
现在,我们考虑上次对角占优的情况。设 J 为交换对角矩阵,在交叉对角上为 1 (从左下到右上),其 他地方为 0,即
0 0 0 1 0 0 1 0
J
=
,
(1.8)
0 1 0 0
1 0 0 0
因为 A 是上次对角占优的,=A J= AJ Tritoep (γ ,α , β ) 是下次对角占优的。因此,我们可以将原来的线性
对于计算复杂度,我们的算法需要大约 12n (flops)的浮点运算,选主元的 LU 分解法需要 13n (flops) 浮点运算,我们的方法与选主元的 LU 分解法相比需要更少的浮点运算量。对于内存存储空间,我们的 算法只需要存储 2 个大小为 n 向量,少于选主元的 LU 分解法需要存储 5 个大小为 n 向量。特别是,我 们的算法需要较少的数据传输。在数据传输过程中它只需要读取 1 个向量即右边的向量并写入一个向量 (即方程的解)。但是选主元的 LU 分解法需要读取个向量。正如我们所知,现代计算机有多层的内存结构, 有不同的存储级别,如较小的高速缓存和较大的低速磁盘存储。在计算过程中,数据在不同级别的高速 缓存中传输。因此,较少数据传输的算法可能会显示出更好的计算性能。这使得该算法比选主元的 LU 分解方法更有效。