%%% Finite Elements: Linear (1-D) Diffusion with value-value B.C.
%%% BENG 221 example, 11/17/2011

L = 1;    % length
n = 100;  % number of finite elements
ne = 10;  % number of displayed eigenmodes
don = diag(ones(n,1),0); % on-diagonal, all ones
doff = diag(ones(n-1,1),1); % off-diagonal, all ones

dx = 1 / (n+1); % grid spacing, uniform
x = (0:n+1)' * L / (n+1);  % vertex points, including 0 and L boundary points
K = 1/dx * (2*don - doff - doff'); % stiffness matrix
M = dx/6 * (4*don + doff + doff'); % mass matrix


% finite differences: eigenvalue-eigenmode problem K*U = Lambda*dx*U with Lambda diagonal
help eig
% [U, Lambda] = eig(K/dx);
[U, Lambda] = eig(K, dx*don);
lambda_fd = diag(Lambda); % eigenvalues
phi_fd = [zeros(1,n); U; zeros(1,n)]; % eigenmodes; including 0 and L boundary points
phi_fd = phi_fd * diag(sign(phi_fd(2,:))); % flip eigenmodes with negative slope at zero

% finite elements: generalized eigenvalue-eigenmode problem K*U = Lambda*M*U with Lambda diagonal
[U, Lambda] = eig(K, M);
lambda_fe = diag(Lambda); % eigenvalues
phi_fe = [zeros(1,n); U; zeros(1,n)]; % eigenmodes; including 0 and L boundary points
phi_fe = phi_fe * diag(sign(phi_fe(2,:))); % flip eigenmodes with negative slope at zero

% analytical solution; ground truth
lambda_anal = (pi * (1:n) / L).^2; % analytical eigenvalues (separation of variables)
phi_anal = sqrt(2) * sin(pi * x * (1:n) / L); % analytical eigenmodes (separation of variables)


% plot eigenvalues
figure(1)
semilogy(lambda_fd, 'k:')
hold on
semilogy(lambda_fe, 'k-')
semilogy(lambda_anal, 'k--')
hold off
xlabel('Mode{\it k}')
ylabel('Eigenvalue {\it \lambda_k}')
legend('Finite Differences', 'Finite Elements', 'Analytical', 'Location', 'NorthWest')
print -dpdf 'linfinelem_eigenvalues.pdf'

% plot eigenmodes
figure(2)
plot(x, phi_fd(:,1:ne), ':')
hold on
plot(x, phi_fe(:,1:ne), '-')
plot(x, phi_anal(:,1:ne), '--')
hold off
xlabel('{\it x/L}')
ylabel('Eigenmodes {\it \phi_k}')
print -dpdf 'linfinelem_eigenmodes.pdf'
