% 
% 测试 SOR 对参数 omega 的敏感性
%

clear all;
close all;
n = 16*2;
h = 1/(n+1); % step size

% initial guess
u0 = zeros(n+2,n+2);

% boundary condition
x = 0:h:1;  y = 0:h:1;
u0(1,:) = 0.25 * (y .* y);
u0(n+2,:) = 0.25 * (1 + y .* y);
u0(:,1) = 0.25 * (x .* x);
u0(:,n+2) = 0.25 * (x .* x + 1);

f = -1; % right-hand side f(x,y)

% true solution
[X,Y] = meshgrid(x);
ut = 0.25*(X.*X + Y.*Y);
norm_ut = norm(ut(:));

tol = 1e-6; % stopping criteria
iter_max = 1000; % maximun number of iterations

hh = h*h;  % h^2

% values of omega
Omega = 1.3:0.01:1.9;

Iter = zeros(size(Omega));  % iteration number for different omega

% start SOR
for k_omega = 1 : length(Omega)
    omega = Omega(k_omega);
    u = u0;
    for k = 1 : iter_max
        for i = 2:n+1
            for j = 2:n+1
                vij_old = u(i,j);
                u(i,j) = 0.25* (hh*f + u(i+1,j) + u(i-1,j) ...
                    + u(i,j+1) + u(i,j-1));
                u(i,j) = (1-omega)*vij_old + omega*u(i,j);
            end
        end
        relerr = norm(u(:) - ut(:)) / norm_ut;
        if relerr < tol
            break;
        end
    end
    Iter(k_omega) = k;
end

% plot the result
plot(Omega,Iter,'o-');
axis([1.3,1.9,0,1000])
xlabel('\omega'); ylabel('iteration number');
legend([int2str(n),' x ', int2str(n)]);
title('SOR with different \omega');


