三种迭代法 Three iteration

clear
n=10;
for i = 1:n
    for j = 1:n
        A(i,j)=vpa(1/(i+j-1)); % 创建希尔伯特矩阵
    end
end
b=ones(n,1); % 创建全1的右端向量
I=1000;w=0.5; % 设置最大迭代次数和松弛因子

% 调用三种迭代方法
[x,t,it] = Iteration(A,b,I,eps)
norm(x)
[x,t,it] = GaussSeidel(A,b,I,eps)
norm(x)
[x,t,it,w] = SOR(A,b,I,eps,w)
norm(x)

function [x,t,it] = Iteration(A,b,I,eps)
% Jacobi迭代法
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% eps: 收敛精度
% 输出：
% x: 解矩阵
% t: 计算时间
% it: 实际迭代次数

if nargin < 4
    eps = 1e-6;
end

tic
[n,~] = size(A);
B = zeros(size(A));
x = zeros(n,1);
f = zeros(n,1);

% 构造迭代矩阵B和向量f
for r = 1:n
    for l = 1:n
        if l == r
            B(r,r) = 0;
        else
            B(r,l) = -A(r,l)/A(r,r);
        end
    end
    f(r) = b(r)/A(r,r);
end

x_exact = A\b; % 精确解（用于比较）
it = 1;
for k = 1:I-1
    x = B*x+f;
    if norm(x-x_exact)>eps
        it = it+1;
    end
end
t = toc;
end

function [x,t,it] = GaussSeidel(A,b,I,eps)
% 高斯-赛德尔(Gauss-Seidel)迭代法
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% eps: 收敛精度
% 输出：
% x: 解矩阵
% t: 计算时间
% it: 实际迭代次数

if nargin < 4
    eps = 1e-6;
end

tic
[n,~] = size(A);
x = zeros(n,1);

% 矩阵分解
D = diag(diag(A)); % 对角矩阵
L = -tril(A,-1);   % 严格下三角部分（取负）
U = -triu(A,1);    % 严格上三角部分（取负）
G = (D-L)\U;       % 迭代矩阵
f = (D-L)\b;       % 迭代向量

x_exact = A\b; % 精确解（用于比较）
it = 1;
for k = 1:I-1
    x = G*x+f;
    if norm(x_exact-x)>eps
        it = it+1;
    end
end

t = toc;
end

function [x,t,it,w] = SOR(A,b,I,eps,w)
% 逐次超松驰迭代法(successive over relaxation method)
% 输入：
% A: 系数矩阵
% b: 载荷矩阵
% I: 最大迭代次数
% eps: 收敛精度
% w: 松弛因子(w=1 时即为 Gauss-Seidel 迭代法)
% 输出：
% x: 解矩阵
% t: 计算时间
% it: 实际迭代次数
% w: 实际使用的松弛因子

tic
[n,~] = size(A);
x = zeros(n,1);

% 矩阵分解
D = diag(diag(A)); % 对角矩阵
L = -tril(A,-1);   % 严格下三角部分（取负）
U = -triu(A,1);    % 严格上三角部分（取负）

% 计算最佳松弛因子
w_opt = 2/(1+sqrt(1-(max(abs(eig(D\(L+U))))))^2);

if nargin < 4
    eps = 1e-6;
    w = w_opt;
end
if nargin < 5
    w = w_opt;
end

% 构造迭代矩阵和向量
Lw = (D-w*L)\((1-w)*D+w*U);
f = w*((D-w*L)\b);

x_exact = A\b; % 精确解（用于比较）
it = 1;
for k = 1:I-1
    x = Lw*x+f;
    if norm(x-x_exact)>eps
        it = it+1;
    end
end

t = toc;
end