%
% @file gs_fd.m
%
% @brief
%
% @copyright
% Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
% SPC (Swiss Plasma Center)
%
% spclibs is free software: you can redistribute it and/or modify it under
% the terms of the GNU Lesser General Public License as published by the Free
% Software Foundation, either version 3 of the License, or (at your option)
% any later version.
%
% spclibs is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <https://www.gnu.org/licenses/>.
%
% @authors
% (in alphabetical order)
% @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
%
clear all

N=1024; 
L=100;
n0=1;
Narr=[128 256 512 1024];

for i=1:length(Narr);
    N=Narr(i);
    h=L/N; h2=h^2;
    title_str=sprintf(['N=%d, L=%.1f, n0=%.1f'], N, L, n0);
    v0=-2.+n0*h2;
    u0=1.;
    dom(i)= 2*abs(u0)/abs(v0);
    v=v0.*ones(N,1);
    u=u0*ones(N-1,1);
    mata=diag(u,1) + diag(u,-1) + diag(v);
    matl= tril(mata,0); % D+L
    lambda = eig(-triu(mata,1),matl);
    rho(i)=max(abs(lambda));
    
% $$$     figure
% $$$     plot(lambda,'o')
% $$$     xlabel('Real eigenvalue'), ylabel('Imag eigenvalue')
% $$$     grid on
% $$$     axis equal
% $$$     title_str=sprintf(['N=%d, L=%.1f, n0=%.1f, Spec. Radius=%.4f, dom=%.4f'], N, L, ...
% $$$                       n0, rho(i),dom(i));
% $$$     title(title_str);
    fprintf(1, 'Spectral Radius of GS relaxation matrix = %.4f, dom=%.4f\n', ...
            rho(i),dom(i))
end

figure
plot(Narr, rho,'o-')
xlabel('N'); ylabel('GS spectral radius')
title_str=sprintf(['FD scheme, L=%.1f, n0=%.1f'], L, n0);
title(title_str);
grid on