exp3_4.m

使用matlab软件编制的计算程序

% exp3_4.m --- SOR迭代法收敛速度受松驰因子的影响试验

function try_sor_and_relaxation
% 参见 P64例8 和 P65表3-4
A = [ 2  -1   0   0 
     -1   2  -1   0
      0  -1   2  -1
      0   0  -1   2];
b = [1 0 1 0]';
tol = 1e-6;
maxiter = 1000;
P = [1 0 1 0]';

clc
fprintf('** SOR迭代:松驰因子与迭代次数的关系 **\n'),
fprintf('\n    松驰因子         迭代次数')
fprintf('\n---------------------------------')

for W = 1:0.1:1.9
    [X,iternum] = sor(A,b,tol,maxiter,P,W);
    fprintf('\n      %3.2f            %4.0f',W,iternum)
end
% 通过以上结果,你认为最佳松驰因子大致为多少?

% ----------  说明  ----------
% 对于对称正定三对角矩阵(上面矩阵就是),其最佳松驰因子是可以求得的
%                  Wopt = 2 / ( 1+sqrt(1-rho^2) )
%                  其中rho是Jacobi迭代矩阵 Mj 的谱半径
% 下面求之,应该与我们实验的结果差不多吧,看看
D = diag(diag(A));
Mj = eye(4) - inv(D)*A;
rho = max(abs(eig(Mj)));
Wopt = 2 / ( 1+sqrt(1-rho^2) );
fprintf('\n---------------------------------')
fprintf('\n理论上最佳松驰因子是 Wopt = %3.2f',Wopt)



% ---------- SOR 迭代法 -----------
function [X,iternum] = sor(A,B,tol,maxiter,P,W)
% [X,iternum] = sor(A,B,tol,maxiter,P,W) SOR迭代法解线性方程组 AX=B
% 输入 ---- tol:     相对残向量的容差,当norm(B-AX)/norm(B) <= tol 终止迭代
%           maxiter: 最大迭代次数;
%           P:       初值
%           W:       松驰因子,必须 0 < W < 2(当 W = 1 时为G-S迭代法)
% 输出 ---- X:       解向量
%           iternum: 收敛迭代次数

N = length(B);  X = P;  tol = tol*norm(B);
for k = 1:maxiter
    for i=1:N
        if i == 1
            X(1) = (B(1)-A(1,2:N)*P(2:N))/A(1,1);
        elseif i == N
            X(N) = (B(N)-A(N,1:N-1)*X(1:N-1))/A(N,N);
        else
            X(i) = (B(i)-A(i,1:i-1)*X(1:i-1)-A(i,i+1:N)*P(i+1:N))/A(i,i);
        end
        X(i) = (1-W)*P(i) + W*X(i);
    end
    if norm(B-A*X) <= tol, break,end;
    P=X;
end
iternum=k;
return
% --------------------------------------------