数值线性代数最重要和最常见的应用之一是线性系统的解可以用这种形式表示A*x=b
.当一个
是一个大型稀疏矩阵,可以使用迭代方法求解线性系统,这使您能够在计算的运行时间和解的精度之间进行权衡。本主题介绍MATLAB中可用的迭代方法®解方程A*x=b
.
求解线性方程组有两种方法A*x=b
:
直接的方法是高斯消去法的变体。这些方法通过矩阵运算(如LU、QR或Cholesky分解)直接使用单个矩阵元素。您可以使用直接方法以高精度求解线性方程组,但在处理大型稀疏矩阵时,这些方法可能会变慢。求解直线的速度采用直接法的ar系统强烈依赖于系数矩阵的密度和填充模式。
例如,此代码求解一个小型线性系统。
A=幻数(5);b=和(A,2);x=A\b;范数(A*x-b)
ans = 1.4211 e-14
迭代的方法在有限步数后,生成线性系统的近似解。这些方法对于大型方程组非常有用,因为在这些方程组中,为了缩短运行时间而牺牲精度是合理的。迭代方法仅通过矩阵向量积或抽象线性算子间接使用系数矩阵。迭代方法可以用于任何矩阵,但它们通常适用于直接求解速度较慢的大型稀疏矩阵。与直接法相比,用间接法求解线性系统的速度对系数矩阵填充模式的依赖性不强。但是,使用迭代方法通常需要针对每个特定问题调整参数。
例如,这段代码解决了一个具有对称正定系数矩阵的大型稀疏线性系统。
一个= delsq (numgrid (“我,400)); b=一(尺寸(A,1),1);x=pcg(A,b,[],1000);范数(b-A*x)
pcg在迭代796收敛到一个相对残差为9.9e-07的解。ans=3.4285e-04
MATLAB实现了各种迭代方法,根据系数矩阵的属性,这些方法具有不同的优缺点一个
.
如果有足够的存储空间,直接方法通常比间接方法更快,更普遍适用。一般情况下,您应该尝试使用x = A \ b
第一。如果直接求解太慢,那么你可以尝试使用迭代方法。
求解线性方程组的大多数迭代算法遵循类似的过程:
从解向量的初始猜测开始x0
.(这通常是一个零向量,除非你指定一个更好的猜测。)
计算剩余范数res =规范(b * x0)
.
将残余量与规定的公差进行比较。如果res<=tol
,结束计算并返回的计算答案x0
.
应用A * x0
并更新向量的大小和方向x0
根据剩余的价值和其他计算的数量。这是完成大多数计算的步骤。
重复步骤2到步骤4,直到x0
是足以满足宽容的。
迭代方法的不同之处在于它们如何更新数据的大小和方向x0
在第4步中,有些人在第2步和第3步中的收敛标准略有不同,但这抓住了所有迭代解算器遵循的基本过程。
MATLAB有几个函数可以实现线性方程组的迭代方法。这些方法旨在解决一个x=b或者最小化范数||b- - - - - -一个x||.其中一些方法有相似之处,并且基于相同的底层算法,但每种算法在某些情况下都有好处[1],[2].
描述 |
笔记 |
---|---|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
MATLAB中迭代求解器的流程图给出了每个求解器有用的大致情况。您通常可以使用巨磁电阻
对于几乎所有的平方非对称问题。在某些情况下,双共轭梯度算法(bicg
,bicgstab
,研究生院理事会
等等)比巨磁电阻
,但它们不可预测的收敛行为往往会导致巨磁电阻
一个更好的初始选择。
迭代法的收敛速度依赖于系数矩阵的谱(特征值)。因此,通过将线性系统转换为更有利的谱(聚类特征值或1附近的条件数),可以提高大多数迭代方法的收敛性和稳定性。这种转换是通过应用第二个矩阵a来实现的预调节器,给系统。这个过程转换了线性系统
转化为等效系统
理想预调节器对系数矩阵进行变换一个转化为单位矩阵,因为任何迭代方法在有这样一个前置条件的一次迭代中都会收敛。在实践中,找到一个好的预调节器需要权衡。转换以三种方式之一进行:左预处理、右预处理或分裂预处理。
第一种情况叫做离开了预处理由于预条件矩阵米出现在。的左边一个:
这些迭代解算器使用左预处理:
在正确的预处理,米出现在…的右边一个:
这些迭代求解器使用了正确的预处理:
最后,对于对称系数矩阵一个
,分割预处理确保转换后的系统仍然是对称的
获取拆分,并且因子显示在一个:
分裂预条件系统的求解算法是基于上述方程,但在实际应用中不需要计算H。解算器算法乘以米
直接。
这些迭代解算器使用分割预处理:
在所有情况下,都是前置条件米,以加速迭代法的收敛。当迭代解决方案的残余误差在迭代之间停滞不前或进展甚微时,这通常意味着您需要生成一个预处理矩阵来合并到问题中。
MATLAB中的迭代求解器允许您指定单个预条件器矩阵米,或两个预条件矩阵因子米=米1米2.这使得以因式分解的形式指定前置条件变得很容易,例如米=lU. 注意,在拆分预处理的情况下,其中米=HHT也认为,两者之间没有关系M1
和平方米
输入和H因素。
在某些情况下,预置条件在给定问题的数学模型中很自然地出现。在没有自然预因子的情况下,您可以使用该表中的一个不完全分解来生成预因子矩阵。不完全因子分解本质上是不完全的直接解,计算起来很快。
看到不完全分解有关伊鲁
和ichol
.
考虑拉普拉斯方程在二维方域上的五点有限差分近似。下面的命令使用带有预处理因子的预处理共轭梯度法M = L * L '
,在那里l
的零填充不完全Cholesky因子一个
.对于这个系统,pcg
不能在不指定预处理矩阵的情况下找到解。
一个= delsq (numgrid (“年代”, 250));b = 1(大小(A, 1), 1);托尔= 1 e - 3;麦克斯特= 100;L = ichol(一个);x = pcg (A, b,托尔,麦克斯特,L, L ');
pcg在迭代92收敛到相对残差为0.00076的解。
pcg
需要92次迭代才能达到指定的公差。但是,使用不同的预处理程序可以产生更好的结果。例如,使用ichol
构造修改的不完全Cholesky允许pcg
仅在39次迭代后满足规定公差。
L=ichol(A,结构(“类型”,“nofill”,“michol”,“上”));x=pcg(A,b,tol,maxit,L,L');
PCG在迭代39时收敛到一个相对残差为0.00098的解。
对于计算困难的问题,您可能需要一个比由伊鲁
或ichol
直接。例如,您可能希望生成一个质量更好的预处理程序或最小化正在进行的计算量。在这些情况下,您可以使用平衡使系数矩阵对角占优(可以得到更好的预处理器)重新排序最小化矩阵因子中的非零数(这可以减少内存需求,并可能提高后续计算的效率)。
如果你同时使用平衡和重新排序来生成一个预处理程序,过程是:
下面是一个使用平衡和重新排序来生成稀疏系数矩阵的预调节器的例子。
创建系数矩阵一个
和一个1的向量b
对于线性方程的右边。计算条件数的估计一个
.
负载west0479;一个= west0479;b = 1(大小(A, 1), 1);康德(A)
ans=1.4244e+12
使用平衡
改进了系数矩阵的条件数。
[P,R,C]=平衡态(A);重新=R*P*A*C;b新=R*P*b;冷凝(重新)
ans=5.1042e+04
重新排列平衡矩阵使用解剖
.
q=解剖(重新);重新=重新(q,q);B新=B新(q);
使用不完全逻辑单元分解生成预处理程序。
[L U] = ilu(重新);
求解线性方程组巨磁电阻
使用预条件矩阵,公差为1平台以及
,最大外部迭代次数为50次,内部迭代次数为30次。
托尔= 1平台以及;麦克斯特= 50;重启= 30;[xnew, flag, relres] = gmres(new,bnew,restart,tol,maxit,L,U);x (q) = xnew;x = C * x (:);
现在,比较一下relres
返回的相对残差巨磁电阻
(包括预处理子)到没有预处理子的相对残差新的
和不平衡的相对剩余量res
结果表明,即使线性系统都是等价的,但不同的方法对每个元素应用不同的权重,这会显著影响残差的值。
reres renew = norm(new*xnew - bnew) / norm(bnew) res = norm(A*x - b) / norm(b)
Relres = 8.7537e-11 resnew = 3.6805e-08 res = 5.1415e-04
MATLAB中的迭代求解器不支持需要你提供了一个数值矩阵一个
.因为求解器执行的计算使用矩阵向量乘法的结果A*x
或‘* x
,您可以提供一个函数来计算这些线性操作的结果。计算这些量的函数通常称为A线性算子.
除了用线性算子代替系数矩阵一个
,也可以使用线性运算符而不是矩阵作为预条件器米
.在这种情况下,函数需要计算M \ x
或M ' \ x
,如解算器参考页上所示。
使用线性运算符可以利用一个
或米
与解算器显式使用矩阵执行完整矩阵向量乘法相比,更高效地计算线性运算的值。这也意味着您不需要内存来存储系数或预处理矩阵,因为线性运算符通常计算矩阵向量乘法的结果n而根本不形成矩阵。
例如,考虑系数矩阵
A = [2 -1 0 0 0 0;1 2 1 0 0 0;0 -1 2 -1 0 0;0 0 -1 2 -1 0;0 0 0 -1 2 -1;0 0 0 1 2];
当一个
乘以一个向量,得到的向量中的大多数元素都是零。结果中的非零元素与向量中的非零三对角元素相对应一个
. 对于给定的向量x
,线性算子函数只需将三个向量相加即可计算出A*x
:
函数y=linearPropertora(x)y=-1*[0;x(1:end-1)]…+2*x…+-1*[x(2:end);0];结束
大多数迭代求解器都要求线性算子函数一个
返回的值A*x
.同样地,对于预处理矩阵米
,函数通常必须计算M \ x
.为解决最小二乘QR分解算法
,qmr
,及bicg
,线性运算符函数还必须返回‘* x
或M ' \ x
应要求。有关线性运算符函数的示例和说明,请参见迭代解算器参考页。
[1] Barrett, R., M. Berry, t.f. Chan, et al.,线性系统解的模板:迭代方法的构建块,暹罗,费城,1994年。
[2]萨阿德,尤瑟夫稀疏线性方程的迭代方法.出版公司,1996。