% EB_Beam_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:
% I*E/rho/A * d^4 u(x)/dx^4 = omega^2 u(x)
% a <= x <= b
% u(a) = u(b) = 0
%
% IMPORTANT REMARK
% In the case of the flow stability problems, you might have to assemble the operator from multiple matrices
% or also, might have to assemble multiple operators (each on which the
% boundary consitions must be enforced), and find the eigenvalues using the
% command
%%
clear all
close all
clc
%% parameters
% Chebyshev polynomial degree - the number of points is actually N+1
N = 41;
N = 81;
N = 151;
% domain boundaries
a = 0;
b = 1;
% material properties
ly = 1e-2; % heigth of the beam
lz = 1e-2; % transversal size of the beam
A = ly*lz; % area
I = lz*ly^3/12; % second moment of area
rho = 7800; % density
E = 200*1e9; % Elasticity modulus
matpar = I*E/rho/A; % combining parameters into a single constant
%% 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;
D3 = D1^3;
D4 = D1^4;
%% boundary conditions and eigenvalue calculation
% our differential operator is the 4th derivative, multiplied with the material constans
A = D4*matpar;
% Boundary conditions enforced with the equation C*w=f
% 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 2 boundary conditions at each of the boundaries, we need to
% remove 4 points from the calculation
% the best way to do this is to remove 2-2 points at the 2 ends of the domain
rind = [1,2,N,N+1]; % BC enforced at the two endpoints
kind = 3:N-1; % 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 = [ ...
D2(1,:); ... % u''(b) = 0; enforced at the 1st point
D3(1,:); ... % u'''(b) = 0; enforced at the 2nd point
D1(end,:); ... % u'(a) = 0; enforced at the (N)th point
zeros(1,N),1 ... % u(a) = 0; enforced at the (N+1)th point
];
% since the boundary conditions are homogeneous, the function is zero at the ends
% this is not used to simplify the code
f = zeros(4,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;
% creating the matrix, which contains the boundary conditions
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');
% usind index vector for sorting the eigenfunctions
V = V(:,ind);
% the boundary points can be calculated by ur = Guk + Hf
% note that f is zero
Vr = zeros(4,N-3);
for i = 1:N-3
Vr(:,i) = G*V(:,i) + H*f;
end
V = [Vr(1:2,:);V;Vr(end-1:end,:)];
%% plotting the eigenvalues
figure;
subplot(1,2,1);
semilogx(real(eigvals),imag(eigvals),'Marker','o','Linestyle','none');
hold on;
xlabel('$\mathrm{log}(\mathcal{R}(\lambda))$','interpreter','latex');
ylabel('$\mathcal{I}(\lambda)$','interpreter','latex');
title('eigenvalue spectrum');
axis([0 1e10 -1 1]);
set(gcf,'Position',[680 647 931 331]);
% with left click, you can continue to plot the eigenvalues, right clicking stops this
buttonvar = 1;
while buttonvar==1
subplot(1,2,1);
% getting an eigenvalue
[xx,yy,buttonvar] = ginput(1);
% calculating complex eigenvalue from values
tempeigv = xx + 1i*yy;
% finding eigenvalue closest to the click
[~, tempind] = min(abs(tempeigv - eigvals));
% delete previous plot
if exist('plth') delete(plth); end
% plotting selected eigenvalue
plth = plot(real(eigvals(tempind)),imag(eigvals(tempind)),'x','LineWidth',1.5,'MarkerSize',11,'Color','k');
% plotting selected eigenvector
subplot(1,2,2);
plot(x,V(:,tempind),'LineWidth',1.3, 'Marker','x');
temp_omega = sqrt(real(eigvals(tempind)));
title(['\omega = ', num2str(temp_omega),', f = ', num2str(temp_omega/2/pi)]);
end
%% 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