%metoda Jacobiego
clc
clear all
%funkcja
F = @(x,y) -cos(x+y)-cos(x-y);
%rozw analityczne
z=@(x,y) cos(x).*cos(y);

%warunki brzegowe
u1 = @(x) cos(x);
u2 = @(y) -cos(y);
u3 = @(x) 0;
u4 = @(y) cos(y);

%przedział omega
a=0;
b=pi;
c=0;
d=pi/2;

n=5;
m=5;

h=(b-a)/(n+1);
k=(d-c)/(b-a)*(n+1)-1;
x=[a:h:b];
y=[c:h:d];

tol=1e-6; %tolerancja
error = 1; %błąd
licznik=0; %licznik iteracji

%wartość lambda
dlgx=length(x)
dlgy=length(y)
lambda = (0.25)*((cos(pi/(dlgx-1))+cos(pi/(dlgy-1)))^2);
%wartość omega
w = 1+(lambda/((1+sqrt(1-lambda))^2));

%tworzenie macierzy
U1(1:n+2) = u1(x);
U2(1:k+2) = u2(y(1:(k+2)));
U3(1:n+2) = u3(x);
U4(1:k+2) = u4(y(1:(k+2)));

U(1,:) = U1(1:n+2);
U(k+2,:) = U3(1:n+2);
U(:,1) = U4(1:k+2);
U(:,n+2) = U2(1:k+2);

Uk=U;



while error>tol
    licznik = licznik+1;
    for i=2:k+1
        for j=2:n+1
             Uk(i,j) =(0.25*w*(U(i,j+1)+U(i+1,j)+Uk(i,j-1)+Uk(i-1,j)))+((1-w)*U(i,j))-0.25*h^2*F(x(j),y(i));
        end
    end

   % error = sqrt(sum(sum((Uk-U).^2)))
    error = max(max(abs(Uk-U)));

    U=Uk;
end
U;
[X,Y] = meshgrid(x,y);
subplot(1,2,1)
surf(X,Y,U)
title('Obliczone')
subplot(1,2,2)
surf(X,Y,(z(X,Y)))
title('Analitycznie')

bb = abs(U-z(X,Y));
max(max(bb))
licznik