526互联

2023-06-30《计算方法》- 陈丽娟 - 线性方程组的迭代解法.md

发布时间 2023-06-30 18:55:26作者: XNEF

2023-06-30《计算方法》- 陈丽娟 - 线性方程组的迭代解法

Matlab计算方法JacobiGauss-SeidelSORSSOR定常迭代法

所谓迭代法实际上是求解一个关于映射的不动点问题：


然后利用构造一个迭代格式

这里表示T的一个复合函数, 其可能随迭代次数而改变，最终目标即是得到. 下面我们给出收敛的定义：
定义1
设, , 若, , 则称点列收敛于, 记作.

定理1
序列依坐标收敛(依范数收敛)于的充分条件是

1. 雅可比迭代法

对于非奇异线性方程组


令

其中, 为的下三角矩阵，为的上三角矩阵。
则有

给定, 下述迭代格式称为雅可比迭代：

注意上式很容易可以写出个分量的迭代格式：

2. 高斯-赛德尔迭代

在雅可比迭代中，如果先被计算出来，则我们可以将分量替换带入到的计算中，即是


其矩阵形式为

即是

3. Jacobi和G-S迭代的统一框架

矩阵分裂
设非奇异，称


为的一个矩阵分裂，其中非奇异。

对于


经过矩阵分裂可得

移项可得

由上可知雅可比迭代就是取, , 而G-S迭代是取, , 因此都可称为矩阵分裂迭代法。

4. 其他矩阵分裂迭代法

4-1. SOR 迭代法（Successive Overrelaxation）

为了加速迭代的收敛，SOR方法将G-S迭代中与做一个凸组合作为下一步迭代的值，即是


上式整理之后可得

同样可以得到对应的


迭代格式为

其中不同取值时：

  1. 当时，即是G-S迭代法；
  2. 当时，称为低松弛方法；
  3. 当时，称为超松弛方法（大部分情况下超松弛更好）。

注意到一开始SOR方法实际上是一个MANN迭代（）,但是MANN迭代需要等条件，且MANN迭代的提出并非是为了加快收敛速度，而是为了解决Picard迭代在非扩张映射下不一定收敛而提出的。

4-2. SSOR 迭代法

将 SOR 方法中的 和 相互交换位置, 则可得迭代格式


结合上式和SOR则得到SSOR（对称超松弛迭代法）

注
对于某些特殊问题, SOR 方法不收敛, 但仍然可能构造出收敛的 SSOR 方法.
一般来说, SOR 方法的渐进收敛速度对参数比较敏感, 但 SSOR 对参数不是很敏感.

5. 收敛性分析

从不动点理论的角度来看，上述所有的迭代格式都可以看作是, 其中, 即所有的迭代都是Picard迭代，那么根据Banach不动点定理，我们知道当是压缩映射时(即)，序列收敛到的不动点，也就是线性方程组的解。

下面根据华东师范大学潘建瑜老师《矩阵计算讲义》列出矩阵下的收敛性分析，我们可以看到虽然描述的方式不一样，但是和不动点理论的结果却是相似的。

矩阵序列的收敛

设是中的一个矩阵序列。如果存在矩阵使得


则称收敛到, 即是的极限，记为
.

定理
设矩阵，则当且仅当.
证明：
设，则存在矩阵范数使得。因此
必要性：由即可得。

定理
设, 则对任意矩阵范数， 有

收敛速度
设点列收敛，且. 若存在一个有界常数, 使得


则称点列是p次（渐进）收敛的. 若或且，则称点列是超线性收敛的。

收敛性定理
对任意迭代初始向量, 矩阵分裂迭代法收敛充要条件是.

定义
设是迭代矩阵，则矩阵分裂迭代法的平均收敛速度定义为


渐进收敛速度定义为

定理
考虑矩阵分裂迭代法，如果存在某个算子范数使得,则

6. 各迭代法之间的数值比较

下面给出的算法包含求逆运算，用行迭代是否可以加速计算有待深入。
Jacobi迭代法

  1. function [x0, xlog] = JacobiLE(x0, A, b, epsilon, maxit) 
  2. if(nargin < 5) 
  3. maxit = 1e3; 
  4. end 
  5. iter = 1; 
  6. xlog = zeros(maxit,1); 
  7. errors = norm(A*x0 - b); 
  8. D = diag(diag(A)); 
  9. LU = D - A; 
  10. D = inv(D); 
  11. G = D * LU; 
  12. F = D * b; 
  13. while iter <= maxit && errors >= epsilon 
  14. %x1 = G * x0 + F; 
  15. x1 = x0 + D * (b - A*x0); 
  16. errors = norm(A * x1 - b); 
  17. xlog(iter) = errors; 
  18. x0 = x1; 
  19. iter = iter + 1; 
  20. end 
  21. xlog = xlog(1:(iter-1)); 
  22. end 

G-S迭代法

  1. function [x0, xlog] = GSLE(x0, A, b, epsilon, maxit) 
  2. if(nargin < 5) 
  3. maxit = 1e3; 
  4. end 
  5. iter = 1; 
  6. xlog = zeros(maxit,1); 
  7. errors = norm(A*x0 - b); 
  8. D = diag(diag(A)); 
  9. L = D - tril(A); 
  10. U = D - triu(A); 
  11. DL = inv(D-L); 
  12. G = DL * U; 
  13. F = DL * b; 
  14. while iter <= maxit && errors >= epsilon 
  15. x1 = G * x0 + F; 
  16. %x1 = x0 + D * (b - A*x0); 
  17. errors = norm(A * x1 - b); 
  18. xlog(iter) = errors; 
  19. x0 = x1; 
  20. iter = iter + 1; 
  21. end 
  22. xlog = xlog(1:(iter-1)); 
  23. end 

SOR迭代法

  1. function [x0, xlog] = SOR(x0, A, b, omega, epsilon, maxit) 
  2. if(nargin < 6) 
  3. maxit = 1e3; 
  4. end 
  5. iter = 1; 
  6. xlog = zeros(maxit,1); 
  7. errors = norm(A*x0 - b); 
  8. D = diag(diag(A)); 
  9. L = D - tril(A); 
  10. U = D - triu(A); 
  11. DL = inv(D-omega *L); 
  12. G = DL * ((1-omega)*D+omega*U); 
  13. F = omega * DL * b; 
  14. while iter <= maxit && errors >= epsilon 
  15. x1 = G * x0 + F; 
  16. %x1 = x0 + D * (b - A*x0); 
  17. errors = norm(A * x1 - b); 
  18. xlog(iter) = errors; 
  19. x0 = x1; 
  20. iter = iter + 1; 
  21. end 
  22. xlog = xlog(1:(iter-1)); 
  23. end 

SSOR迭代法

  1. function [x0, xlog] = SSOR(x0, A, b, omega, epsilon, maxit) 
  2. if(nargin < 6) 
  3. maxit = 1e3; 
  4. end 
  5. iter = 1; 
  6. xlog = zeros(maxit,1); 
  7. errors = norm(A*x0 - b); 
  8. D = diag(diag(A)); 
  9. L = D - tril(A); 
  10. U = D - triu(A); 
  11. DL = inv(D-omega*L); 
  12. DU = inv(D-omega*U); 
  13. GL = DL * ((1-omega)*D+omega*U); 
  14. GU = DU * ((1-omega)*D+omega*L); 
  15. FL = omega * DL * b; 
  16. FU = omega * DU * b; 
  17. while iter <= maxit && errors >= epsilon 
  18. x1 = GL * x0 + FL; 
  19. x1 = GU * x1 + FU; 
  20. %x1 = x0 + D * (b - A*x0); 
  21. errors = norm(A * x1 - b); 
  22. xlog(iter) = errors; 
  23. x0 = x1; 
  24. iter = iter + 1; 
  25. end 
  26. xlog = xlog(1:(iter-1)); 
  27. end 

下面给出例子

  1. clear 
  2. n = 1000; 
  3. x0 = rand(n,1); 
  4. A = rand(n,n)*10; 
  5. A = A + diag(sum(A,2)); 
  6. max(abs(eig(A))); %谱半径 
  7. b = rand(n,1); 
  8. epsilon = 1e-6; 
  9. maxit = 1e3; 
  10. [JacpbiX, JacobiXlog] = JacobiLE(x0, A, b, epsilon, maxit); 
  11. plot(log10(JacobiXlog)) 
  12. xlim([0 100]) 
  13. hold on 
  14. [GSS, GSSlog] = GSLE(x0, A, b, epsilon, maxit); 
  15. plot(log10(GSSlog)) 
  16. hold on 
  17. omega = 1.5; 
  18. [SORX, SORXlog] = SOR(x0, A, b, omega, epsilon, maxit); 
  19. plot(log10(SORXlog)) 
  20. hold on 
  21. omega = 1.5; 
  22. [SSORX, SSORXlog] = SSOR(x0, A, b, omega, epsilon, maxit); 
  23. plot(log10(SSORXlog)) 
  24. legend(["Jacobi" "G-S" "SOR" "SSOR"]) 

结果如下所示
enter description here

随机产生的矩阵很有可能在迭代过程中导致迭代发散，因此这里将随机矩阵改成了严格行对角占优矩阵，相关补充如下。

定义
设, 则有：

  1. 如果的元素满足, 则称为严格(行)对角占优矩阵。（相应的可以定义列对角占优矩阵）
  2. 如果的元素满足, 且至少对一个不等式严格成立，则称为弱对角占优矩阵

定义
设，如果存在置换矩阵， 使得


其中为r阶方阵，为阶方阵，则称为可约矩阵，否则则称A为不可约矩阵。

定理
设，如果为严格对角占优矩阵，或为弱对角占优矩阵且是不可约矩阵，则Jacobi和G-S迭代均收敛。

7. 待研究情况

  1. 除了定理所示情况，还有哪些矩阵能使迭代收敛。
  2. 不可约矩阵的判定方法。
  3. 最佳的松弛因子。
  4. 加速的、或更广泛的迭代算法。