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

合集 - 《计算方法》-陈丽娟(8)

1.2023-06-16 《计算方法》- 陈丽娟 - 绪论.md2023-06-162.2023-06-16《计算方法》- 陈丽娟 - 插值法（一）2023-06-163.2023-06-27《计算方法》- 陈丽娟 - 线性方程组的直接解法.md2023-06-274.2023-06-19《计算方法》- 陈丽娟 - 方程的近似解法（注解）2023-06-195.2023-06-18《计算方法》- 陈丽娟 - 方程的近似解法.md2023-06-186.2023-06-17《计算方法》- 陈丽娟 - 插值法（二）.md2023-06-17

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

8.2023-06-28《计算方法》- 陈丽娟 - 向量和矩阵基础.md2023-06-28

收起

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次（渐进）收敛的. 若或且，则称点列是超线性收敛的。

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

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

渐进收敛速度定义为

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

1. 
2. 
3. 

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"]) 
    

    ---

结果如下所示

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

定义
设, 则有：

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

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

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

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

7. 待研究情况

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