eigenvalue.m

%This function solves for the natural frequencies of the beam when P=0 and
%q=0 for 6,10, and 20 segments 
function [f1,f2,f3] = eigenvalue(N,P)
E = 70E9; %Young's Modulus (Pa)
Rho=2700; %Density (kg/m^3) 
b=0.1; %width of beam(m)
h=0.01; %Beam height (m)
L=1; %beam length (m)
I =(b*(h^3))/12; %Second MOI (m^4)  
Area = b*h; %Area (m^2) 

seg_length = 1/N; %segment length for L=1 

%Diagonal components 
Diagonal = ((P*2)/((N^2)*E*I))+6; %Diagonal 
Off_Diag = (-P/((N^2)*E*I)) - 4; %Off Diagonal 

A = zeros(N-1,N-1); %Matrix A of zeros 
%For loop to set up terms 
for row = 1:(N-1)
    for column = 1:(N-1)
        %solve for diagonal terms
        if row == column
            A(row,column) = Diagonal;
        end
        %solve for off diagonal terms
        if column == row - 1
            A(row,column) = Off_Diag;
        end
        if column == row + 1
            A(row,column) = Off_Diag;
        end
        %Solve for off-off diagonal terms
        if column == row - 2
            A(row,column) = 1;
        end
        if column == row + 2
            A(row,column) = 1;
        end
    end
end

A(1,1) = ((2*P)/((N^2)*(E*I)))+5;
A(N-1,N-1) = ((2*P)/((N^2)*(E*I)))+5;
eigenvalue = eig(A);

dx = L/(N+1);
x = linspace(dx,L-dx,N+1);

%Solving for natural frequency with corresponding values 
f1 = sqrt((eigenvalue(1)*E*I)/Rho/Area/(seg_length^4));
f2 = sqrt((eigenvalue(2)*E*I)/Rho/Area/(seg_length^4));
f3 = sqrt((eigenvalue(3)*E*I)/Rho/Area/(seg_length^4));
setdefaults

%Plot result based on input segment and load 
plot(x,sin(pi*x),x,sin(2*pi*x),x,sin(3*pi*x))
xlabel('x distance (m)')
ylabel('Deflection')

end