SOR法求解方程组Ax=b

原理

矩阵形式的迭代公式：

x^(k+1) = (D - ωL)^-1[(1 - ω)D + ωU] x^(k)+ ω(D - ωL)^-1b

MATLAB实现

function [x,cnt,errory] = sor(A,b,w,x0)
% sor 
% 输入参数：矩阵A b 初值x0
% 输出参数：计算结果x，迭代次数cnt，误差下降图
sizeA = size(A);
if sizeA(1) ~= sizeA(2)
    error('A不是方阵');
end
%% 初始化参数
e = 1e-4;
D = diag(diag(A));
L = D - tril(A);
U = D - triu(A);
M = (D - w * L) \ ((1 - w) * D + w * U);
f = w * ((D - w * L) \ b);
%% 迭代计算
cnt = 1;
x1 = x0;
x2 = M * x1 + f;
errory(cnt) = norm(x2 - x1,1);
while errory(cnt) > e
    x1 = x2;
    x2 = M * x1 + f;
    cnt = cnt + 1;
    errory(cnt) = norm(x2 - x1,1);
end
x = x2;
errory(cnt + 1) = norm(x2 - x1,1);
%% 收敛速度
errorx = 0:1:cnt;
errory = exp(errory);
figure
semilogy(errorx,errory);
end