fem.m

% Ilość elementów
N = 1000;


h = 3 / N;
G = 6.6743e-11;

function result = integrate_half(a, b, f1, i1, f2, i2)
    % Liczba węzłów: 2
    % Funkcje bazowe są liniowe na całkowanych przedziałach (połowa przedziału dla elementu), więc kwadratura jest dokładna
    t1 = (a + b) / 2 + (b - a) / 2 * (-1/sqrt(3));
    t2 = (a + b) / 2 + (b - a) / 2 * (1/sqrt(3));

    if nargin < 5
        % Dla jednej funkcji
        result = ((b - a) / 2) * (f1(i1, t1) + f1(i1, t2));
    else
        % Dla iloczynu dwóch funkcji
        result = ((b - a) / 2) * (function_product(f1, i1, f2, i2, t1) + ...
                                  function_product(f1, i1, f2, i2, t2));
    end
end

% Funkcja do całkowania metodą Gaussa-Legendre'a
function result = integrate(a, b, f1, i1, f2, i2)

    if nargin < 5
        % Dla jednej funkcji
        result = integrate_part(a, (a+b)/2, f1, i1) + integrate_half((a+b)/2, b, f1, i1);
    else
        % Dla iloczynu dwóch funkcji
        result = integrate_part(a, (a+b)/2, f1, i1, f2, i2) + integrate_half((a+b)/2, b, f1, i1, f2, i2);
    end
end

% Iloczyn dwóch funkcji
function result = function_product(f1, i1, f2, i2, t)
    result = f1(i1, t) * f2(i2, t);
end

% Funkcja e_i(x)
function result = e(i, x, h)
    if x > i * h - h && x < i * h
        result = x / h - i + 1;
    elseif x > i * h && x < i * h + h
        result = -x / h + i + 1;
    else
        result = 0;
    end
end

% Pochodna e_i(x)
function result = de(i, x, h)
    if x > i * h - h && x < i * h
        result = 1 / h;
    elseif x > i * h && x < i * h + h
        result = -1 / h;
    else
        result = 0;
    end
end

% Funkcja B(e_i, e_j)
function result = B(i, j, h)
    if abs(j - i) > 1
        result = 0;
        return;
    end

    if abs(i - j) == 1
        a = min(i, j) * h;
        b = max(i, j) * h;
    else
        a = i * h - h;
        b = i * h + h;
    end

    result = integrate(a, b, @(k, x) de(k, x, h), i, @(k, x) de(k, x, h), j);
end

% Funkcja L^(e_j)
function result = L_tilde(j, h, G)
    if j * h < 1 || j * h > 2
        result = -(1 / 3) * integrate(j * h - h, j * h + h, @(k, x) de(k, x, h), j);
    else
        a = max(1, (j - 1) * h);
        b = min(2, (j + 1) * h);

        result = 4 * pi * G * integrate(a, b, @(k, x) e(k, x, h), j) - ...
                 (1 / 3) * integrate(j * h - h, j * h + h, @(k, x) de(k, x, h), j);
    end
end

% Tworzenie macierzy A
function A = create_A(N, h)
    A = zeros(N - 1, N - 1);

    for i = 1:(N - 1)
        A(i, i) = B(i, i, h); % Diagonalne
        if i < N - 1
            A(i, i + 1) = B(i, i + 1, h); % Górny element
            A(i + 1, i) = A(i, i + 1);   % Dolny element (symetryczne)
        end
    end
end

% Tworzenie macierzy B
function result = create_B(N, h, G)
    result = zeros(N - 1, 1);
    for j = 1:(N - 1)
        result(j) = L_tilde(j, h, G);
    end
end


% Tworzenie i wyświetlanie macierzy
A = create_A(N, h);
C = create_B(N, h, G);

disp('Równanie: A * X = B')
disp('Macierz A:');
disp(A);

disp('Macierz B:');
disp(C);

disp('Macierz X:')
X = A \ C;
disp(X)

% Obliczanie wartości phi(x)
function result = phi(x, X, h)

    result = zeros(size(x));
    
    for k = 1:length(x)
        temp = 0;
        for i = 1:length(X)
            temp = temp + X(i) * e(i, x(k), h);
        end
        result(k) = temp - (1 / 3) * x(k) + 5;
    end
end


x = linspace(0, 3, 100); % Przedział od 0 do 5 z 100 punktami
y = phi(x, X, h); % Obliczenie wartości funkcji phi(x) dla tych punktów

plot(x, y);
grid on;
xlabel('x');
ylabel('phi(x)');
title('Wykres funkcji phi(x)');  
disp(y)