2023

您所在的位置：网站首页 › 非线性计算不收敛 › 2023

2023

2023-07-01 13:01| 来源: 网络整理| 查看: 265

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

结果如下所示 enter description here