% LaplaceEigenvalues_Diffop_Example.m
%
% example code for flow stability projects using the method of discretizing
% differential operators
%
% in this code the 1D Laplace eigenvalue problem is solved:
% u''(x) = lambda u(x)
% a <= x <= b
% u(a) = u(b) = 0

% the eigenvalues are: lambda = -k^2 pi^2, k = 1,2,...
% the eigenvectors: u_k = sin(k*pi*x), k = 1,2,...
%%
clear all
close all
clc

%% parameters
% Chebyshev polynomial degree - the number of points is actually N+1
N = 81;
% domain boundaries
a = 0;
b = 1;

%% chebyshev differentiation matrices
% chebyshev differentiation matrix and grid, on [-1,1]
[Dcheb,xcheb] = cheb(N);
% mapping chebyshev grid to the computational domain [a,b]
% then mapping differentiation matrix
% algebraic mapping
x=(xcheb+1)*(b-a)/2 + a;
D1 = Dcheb*(2/(b-a));

% rational mapping, in the case the point need to be refined toward one of
% the boundaries
% half the points lie on the interval [a, yhalf], and the other half on [yhalf,b];
% yhalf = 0.2;
% c1 = b*yhalf/(b-2*yhalf);
% c2 = 1 + 2*c1/(b-a);
% y = c1*(1+xcheb)./(c2-xcheb) + a;
% dxi = (c2-xcheb).^2 / (c1*(c2+1));
% D1 = diag(dxi) * Dcheb;

% differentiation matrix for the higher derivatives
D2 = D1^2;

%% boundary conditions and eigenvalue calculation
% our differential operator is simply D2, the second derivative
A = D2;

% First, we need vectors for indexing
%       the points r, at which the boundary conditions are enforced
%       These points are removed from the calculation.
%       At the remaining points with index k, we will do the actual eigenvalue calculation
% 
% since we have 1-1 boundary condition at each of the boundaries, the boundary conditions
%       are enforced near the boundaries, 1-1 point at the end of the computational domain are chosen
rind = [1,N+1]; % BC enforced at the two endpoints
kind = 2:N; % internal points



% Second, we need to assembly the matrix C
% In the matrix C, the boundary conditions are prescribed in the form: C*w = f
% note the y goes from 'b' to 'a' (in descending order), so we should
% prescribe the boundary conditions in this order
% 
% Dirichlet boundary conditions can be prescribed at a given row of the matrix C 
%       by zeroing out that row, and setting the corresponding coloumn of that row to 1
% Boundary conditions for a derivative of the function at a given point can be prescribed
%        by making equal a row of C with the corresponding row of the differentiation matrix
% Inhomogeneous boundary conditions can be handled by prescribing nonzero values in the vector f
C = [ ...
    1, zeros(1,N); ...           % u(b) = 0;
    zeros(1,N), 1  ...           % u(a) = 0;
    ];

% since the boundary conditions are homogeneous, the function is zero at the
% ends; for simplicity, this is not used
f = zeros(2,1);
% reordering C, simply by using the index vectors
C_r = C(:,rind);
C_k = C(:,kind);
% calculating the matrix H and G
H = inv(C_r);
G = -H*C_k;

A_hat = A(kind,kind) + A(kind,rind)*G;


% calculating eigenvalues and eigenfunctions
[V,D] = eig(A_hat);

% sorting eigenvalues, outputs are the sorted eigenvalues and the indices
[eigvals,ind] = sort(diag(D),'descend');

% calculating k integers
k = sqrt(-eigvals)/pi;

% usind index vector for sorting the eigenfunctions
V = V(:,ind);

% the boundary points can be calculated by ur = Guk + Hf
% this is just demonstation, in our case the values are zero
Vr = zeros(2,N-1);
for i = 1:N-1
    Vr(:,i) = G*V(:,i) + H*f;
end
V = [Vr(1,:);V;Vr(2,:)];


%% plotting eigenfunctions
nplot = 5;
figure;
for i = 1:nplot
    subplot(nplot,1,i);
    plot(x,V(:,i));
    title(['k = ',num2str(k(i)),'; The eigenvalue: ',num2str(eigvals(i))]);
    ylabel('u(x)');
end
xlabel('x');
% success: we succesfully recivered the integers k, and also the eigendunctions

%% Chebyshev differentiation function
function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)';
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';
D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries
D = D - diag(sum(D')); % diagonal entries
end