Part4_Bonus.m

% ME 3255 Computational Mechanics
% Peter Joseph Damian, Group 9
% Part 4 Rewritten for bonus problem
clear; clc;

%% Please see issue "Derivation" on github for the derivations of these equations
%% Equation 16: Setting up eigenvalue problem to solve for lowest natural frequency

%% Constants and Properties
% Material Properties
E = 70*(10^6);          % Young's modulus, Pa
p = 2700;               % Density, kg/m^3

% Geometry
b = 0.1;                % base width, m
h = 0.01;               % height, m
L = 1;                  % bar length, m
A = b*h;                % area, m^2
I = (b*(h^3))/12;       % second moment of inertia, m^4

% Frequency squared coefficient
EIpa = (E*I)/(p*A);     % Inverse coefficient of omega^2

%% 20 Segment System of Equations

% Number of Segments
n20 = 20;

% Write Values of P
P = linspace(-1000,1000,2000);
omega = zeros(length(P),1);
for i = 1:length(P)
    Q = P(i);
    % Write w coefficient matrix (see above for details)
    [ w20 ] = Coefficients( n20,EIpa,L,Q,E,I );
    % Calculate eigenvalues of system
    Lambda20 = eig(w20);
    % Calculate natural frequencies of system
    omega20 = Lambda20.^(1/2);
    % Store lowest natural frequency for use outside of loop
    omega(i) = omega20(1);
end

%% Plotting
setdefaults
figure(1)
plot(P,omega)
title('Lowest Nat Freq vs Load P')
xlabel('Load P,[N]')
ylabel('Omega, [1/s]')