公告

View Post

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

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迭代中与做一个凸组合作为下一步迭代的值，即是

上式整理之后可得

同样可以得到对应的

迭代格式为

其中不同取值时：

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

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

4-2. SSOR 迭代法

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

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

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

5. 收敛性分析

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

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

矩阵序列的收敛

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

则称

收敛到

, 即

是

的极限，记为

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

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

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

则称点列

是p次（渐进）收敛的. 若

或

且

，则称点列是超线性收敛的。

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

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

渐进收敛速度定义为

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

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

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

function [x0, xlog] = JacobiLE(x0, A, b, epsilon, maxit) 
    if(nargin < 5) 
        maxit = 1e3; 
    end 
    iter = 1; 
    xlog = zeros(maxit,1); 
    errors =  norm(A*x0 - b); 
    D = diag(diag(A)); 
    LU = D - A; 
    D = inv(D); 
    G = D * LU; 
    F = D * b; 
    while iter <= maxit && errors >= epsilon 
        %x1 = G * x0 + F; 
        x1 = x0 + D * (b - A*x0); 
        errors =  norm(A * x1 - b); 
        xlog(iter) = errors; 
        x0 = x1; 
        iter = iter + 1; 
    end 
    xlog = xlog(1:(iter-1)); 
end 

G-S迭代法

function [x0, xlog] = GSLE(x0, A, b, epsilon, maxit) 
    if(nargin < 5) 
        maxit = 1e3; 
    end 
    iter = 1; 
    xlog = zeros(maxit,1); 
    errors =  norm(A*x0 - b); 
    D = diag(diag(A)); 
    L = D - tril(A); 
    U = D - triu(A); 
    DL = inv(D-L); 
    G = DL * U; 
    F = DL * b; 
    while iter <= maxit && errors >= epsilon 
        x1 = G * x0 + F; 
        %x1 = x0 + D * (b - A*x0); 
        errors =  norm(A * x1 - b); 
        xlog(iter) = errors; 
        x0 = x1; 
        iter = iter + 1; 
    end 
    xlog = xlog(1:(iter-1)); 
end 

SOR迭代法

function [x0, xlog] = SOR(x0, A, b, omega, epsilon, maxit) 
    if(nargin < 6) 
        maxit = 1e3; 
    end 
    iter = 1; 
    xlog = zeros(maxit,1); 
    errors =  norm(A*x0 - b); 
    D = diag(diag(A)); 
    L = D - tril(A); 
    U = D - triu(A); 
    DL = inv(D-omega *L); 
    G = DL * ((1-omega)*D+omega*U); 
    F = omega * DL * b; 
    while iter <= maxit && errors >= epsilon 
        x1 = G * x0 + F; 
        %x1 = x0 + D * (b - A*x0); 
        errors =  norm(A * x1 - b); 
        xlog(iter) = errors; 
        x0 = x1; 
        iter = iter + 1; 
    end 
    xlog = xlog(1:(iter-1)); 
end 

SSOR迭代法

function [x0, xlog] = SSOR(x0, A, b, omega, epsilon, maxit) 
    if(nargin < 6) 
        maxit = 1e3; 
    end 
    iter = 1; 
    xlog = zeros(maxit,1); 
    errors =  norm(A*x0 - b); 
    D = diag(diag(A)); 
    L = D - tril(A); 
    U = D - triu(A); 
    DL = inv(D-omega*L); 
    DU = inv(D-omega*U); 
    GL = DL * ((1-omega)*D+omega*U); 
    GU = DU * ((1-omega)*D+omega*L); 
    FL = omega * DL * b; 
    FU = omega * DU * b; 
    while iter <= maxit && errors >= epsilon 
        x1 = GL * x0 + FL; 
        x1 = GU * x1 + FU; 
        %x1 = x0 + D * (b - A*x0); 
        errors =  norm(A * x1 - b); 
        xlog(iter) = errors; 
        x0 = x1; 
        iter = iter + 1; 
    end 
    xlog = xlog(1:(iter-1)); 
end 

下面给出例子

clear 
n = 1000; 
x0 = rand(n,1); 
A = rand(n,n)*10; 
A = A + diag(sum(A,2)); 
max(abs(eig(A))); %谱半径 
b = rand(n,1); 
epsilon = 1e-6; 
maxit = 1e3; 
[JacpbiX, JacobiXlog] = JacobiLE(x0, A, b, epsilon, maxit); 
plot(log10(JacobiXlog)) 
xlim([0 100]) 
hold on 
[GSS, GSSlog] = GSLE(x0, A, b, epsilon, maxit); 
plot(log10(GSSlog)) 
hold on 
omega = 1.5; 
[SORX, SORXlog] = SOR(x0, A, b, omega, epsilon, maxit); 
plot(log10(SORXlog)) 
hold on 
omega = 1.5; 
[SSORX, SSORXlog] = SSOR(x0, A, b, omega, epsilon, maxit); 
plot(log10(SSORXlog)) 
legend(["Jacobi" "G-S" "SOR" "SSOR"]) 