Add Hayden HW4 files and program flow block diagrams

This commit is contained in:
deamonkai
2026-04-04 17:45:30 -05:00
parent c4d1a468f1
commit 19e5383a43
3 changed files with 297 additions and 0 deletions

207
AERO3220_HW4.m Normal file
View File

@@ -0,0 +1,207 @@
% Hayden Molloy
% AERO3220 HW4
% Due 18 Feb 2026
clear; close all; clc;
%% Initial Parameters
k1 = 100; % N/m
k2 = 50; % N/m
m1 = 500; % kg
m2 = 250; % kg
% Time Parameters
t0 = 0;
tf = 200;
dt = 0.01;
tspan = t0:dt:tf;
% Initial conditions [z1 z1dot z2 z2dot]
x0 = [0; 0; 0; 0];
% Initialize Force Functions
F_1 = @(t) 100; % N
F_2 = @(t) 100*sin(0.25*t); % N
F_3 = @(t) 100*sin(0.50*t); % N
%% Define Cases and ODEs
% CASE 1:
b = 100;
F = F_1;
[t1,x1] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);
% Case 2:
b = 100;
F = F_2;
[t2,x2] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);
% CASE 3:
b = 0;
F = F_2;
[t3,x3] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);
% Case 4:
b = 100;
F = F_3;
[t4,x4] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);
% Case 5:
b = 0;
F = F_3;
[t5,x5] = ode45(@(t,x) eom_two_mass(t,x,m1,m2,k1,k2,b,F), tspan, x0);
%% Report Results
% Create Legend
lgnd = { ...
'Case 1: F=100, b=100', ...
'Case 2: F=100sin(0.25t), b=100', ...
'Case 3: F=100sin(0.25t), b=0', ...
'Case 4: F=100sin(0.50t), b=100', ...
'Case 5: F=100sin(0.50t), b=0'};
% Plot 1: Top mass displacement vs. Time
figure(1);
plot(t1,x1(:,1),'LineWidth',1.3);
hold on;
plot(t2,x2(:,1),'LineWidth',1.3);
plot(t3,x3(:,1),'LineWidth',1.3);
plot(t4,x4(:,1),'LineWidth',1.3);
plot(t5,x5(:,1),'LineWidth',1.3);
grid on;
xlabel('Time (s)');
ylabel('z_1 (m)');
title('Top Mass Displacement z_1(t)');
legend(lgnd,'Location','southoutside');
% Plot 2: Bottom mass displacement vs. Time
figure(2);
plot(t1,x1(:,3),'LineWidth',1.3);
hold on;
plot(t2,x2(:,3),'LineWidth',1.3);
plot(t3,x3(:,3),'LineWidth',1.3);
plot(t4,x4(:,3),'LineWidth',1.3);
plot(t5,x5(:,3),'LineWidth',1.3);
grid on;
xlabel('Time (s)');
ylabel('z_2 (m)');
title('Bottom Mass Displacement z_2(t)');
legend(lgnd,'Location','southoutside');
% Plot 3: Top mass velocity vs. Time
figure(3);
plot(t1,x1(:,2),'LineWidth',1.3);
hold on;
plot(t2,x2(:,2),'LineWidth',1.3);
plot(t3,x3(:,2),'LineWidth',1.3);
plot(t4,x4(:,2),'LineWidth',1.3);
plot(t5,x5(:,2),'LineWidth',1.3);
grid on;
xlabel('Time (s)');
ylabel('dz_1/dt (m/s)');
title('Top Mass Velocity z_1''');
legend(lgnd,'Location','southoutside');
% PLOT 4: Bottom mass velocity vs. Time
figure(4);
plot(t1,x1(:,4),'LineWidth',1.3);
hold on;
plot(t2,x2(:,4),'LineWidth',1.3);
plot(t3,x3(:,4),'LineWidth',1.3);
plot(t4,x4(:,4),'LineWidth',1.3);
plot(t5,x5(:,4),'LineWidth',1.3);
grid on;
xlabel('Time (s)');
ylabel('dz_2/dt (m/s)');
title('Bottom Mass Velocity z_2''');
legend(lgnd,'Location','southoutside');
% PLOT 5: z1(t) and z2(t)
% (solid = z1, dashed = z2)
c1 = [0 114 178]/255;
c2 = [213 94 0]/255;
c3 = [230 159 0]/255;
c4 = [117 112 179]/255;
c5 = [0 158 115]/255;
figure(5);
plot(t1,x1(:,1),'-','LineWidth',1.2, 'Color', c1);
hold on;
plot(t1,x1(:,3),'--','LineWidth',1.2, 'Color', c1);
plot(t2,x2(:,1),'-','LineWidth',1.2, 'Color', c2);
plot(t2,x2(:,3),'--','LineWidth',1.2, 'Color', c2);
plot(t3,x3(:,1),'-','LineWidth',1.2, 'Color', c3);
plot(t3,x3(:,3),'--','LineWidth',1.2, 'Color', c3);
plot(t4,x4(:,1),'-','LineWidth',1.2, 'Color', c4);
plot(t4,x4(:,3),'--','LineWidth',1.2, 'Color', c4);
plot(t5,x5(:,1),'-','LineWidth',1.2, 'Color', c5);
plot(t5,x5(:,3),'--','LineWidth',1.2, 'Color', c5);
grid on;
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Top and Bottom Mass Displacements');
legend({ ...
'Case 1 z_1', ...
'Case 1 z_2', ...
'Case 2 z_1', ...
'Case 2 z_2', ...
'Case 3 z_1', ...
'Case 3 z_2', ...
'Case 4 z_1', ...
'Case 4 z_2', ...
'Case 5 z_1', ...
'Case 5 z_2'}, ...
'Location','eastoutside');
%% Local ODE function
function dx = eom_two_mass(t, x, m1, m2, k1, k2, b, Ffun)
z1 = x(1);
z1dot = x(2);
z2 = x(3);
z2dot = x(4);
u = Ffun(t);
% m1*z1ddot + k1*z1 - k1*z2 = 0
z1ddot = (-k1/m1)*z1 + (k1/m1)*z2;
% m2*z2ddot + b*z2dot + (k1+k2)*z2 - k1*z1 = u
z2ddot = (k1/m2)*z1 - ((k1+k2)/m2)*z2 - (b/m2)*z2dot + (1/m2)*u;
dx = [z1dot; z1ddot; z2dot; z2ddot];
end