%%%       ******************************************
%%% Hilbert Matrix Analysis: Hilbert.m
%%%       ******************************************
%%% MATLAB m-file: Hilbert.m
%%% Date: 26 December 97
%%% Author: Ronald Yannone
%%% For: NOESIS - The Journal of the Mega Society
%%%      ******************************************

flag = 1;

if flag == 1
    grand_sum = zeros(1,20);

    for g = 1:20
        grand_sum(g) = sum(sum(sym(hilb(g)))');
    end

    figure(1)
    clf
    plot(grand_sum)
    grid
    title('Sum of all the elements in the nth Hilbert Matrix')
    xlabel('Nth Hilbert Matrix')
    ylabel('SUM of all Elements in Hilbert Matrix')
    polyfit(g,grand_sum,2);
end

if flag == 4
    gamma = 0.577215664901532860606512; % The Euler-Mascheroni number
    true_harm = zeros(1,50); % Place-holder for first 50
    appr_harm = zeros(1,50);
    true_harm(l) = 1;
    appr_harm(l) = log(l) + gamma + 1/(2*1);

    for k=2:50
        true_harm(k) = true_harm(k-1) + 1/k;
        appr_harm(k) = log(k) + gamma + 1/(2*k);
    end

    figure(1)
    clf
    plot(true_harm)
    hold on
    plot(appr_harm,lr+')
    grid
    title('True [solid line] and "approximation to" [+ line] Harmonic Numbers')
    xlabel('Kth Harmonic Number')
    ylabel('True and Approximate Harmonic NUMBER')
    hold off

    error = true_harm - appr_harm;

    figure(2)
    clf
    plot(error)
    grid title('ERROR: True - Approximation Harmonic Number')
    xlabel('Kth Harmonic Number')
    ylabel('ERROR')

    min_error = min(error)
    max_error = max(error)
end

if flag==3
    suml 1 + 1/(7/3) + 1/(101/30);

    for k = 4:15
        hh sym(hilb(k));
        hh(2:k-1,2:k-l)=zeros(k-2);
        suml = suml + (1/sum(sum(hh)'));
    end
end

if flag == 2
    d = zeros(1,10);

    d(l) = det(sym(hilb(l)));

    for k=2:11
        d(k) = d(k-1) + det(sym(hilb(k)));
    end

    d(l) det(sym(hilb(l)))
    figure(1)
    clf
    plot(d)
    grid
    title('Cumsum of Determinant of Hilbert Matrix')

%     w = zeros(1,20);

%     for k=1:20
%         k;
%         hilb_matrix = sym(hilb(k));
%         w(k) = sum(sum(sym(hilb(k))'));
%         pretty(w);
%     end

%     figure(1)
%     clf
%     plot(w)
%     grid
%     title('The sum of all elements in the Hilbert Matrix')
%     xlabel('Hilbert Matrix of Size N')
%     ylabel('Sum of all elements')

end