11/1/24 2:54 PM MATLAB Command Window 1 of 2
% Given values
k = 1000; % Spring constant in N/m
m = 10; % Mass in kg
delta = 1e-3; % Initial displacement for m/2 mass in meters
% Define the Mass Matrix (M)
M = [m 0 0;
0 m 0;
0 0 m/2];
% Define the Stiffness Matrix (K)
K = [3*k -2*k 0;
-2*k 3*k -k;
0 -k 3*k];
% (c) Calculate Natural Frequencies and Mode Shapes
% Solve the generalized eigenvalue problem K*Phi = lambda*M*Phi
[Phi, Lambda] = eig(K, M);
omega = sqrt(diag(Lambda)); % Natural frequencies
% Display Natural Frequencies
disp('Natural Frequencies (rad/s):');
disp(omega);
% (d) Plot mode shapes
figure;
plot([1 2 3], Phi(:,1), '-o', 'DisplayName', 'Mode 1'); hold on;
plot([1 2 3], Phi(:,2), '-o', 'DisplayName', 'Mode 2');
plot([1 2 3], Phi(:,3), '-o', 'DisplayName', 'Mode 3');
legend;
xlabel('Mass Index');
ylabel('Amplitude');
title('Mode Shapes');
% (e) Define initial conditions and solve for the time response
x0 = [0; 0; delta]; % Initial displacement (only third mass has initial
displacement)
v0 = [0; 0; 0]; % Initial velocity (all masses start at rest)
initial_conditions = [x0; v0]; % Combine initial displacements and velocities
% Define time span
tspan = [0 1]; % Time from 0 to 1 second
% Solve the system using ode45
[t, X] = ode45(@(t, X) system_dynamics(t, X, M, K), tspan, initial_conditions);
% (f) Plot the response of each mass over time
figure;
plot(t, X(:,1), 'r', 'DisplayName', 'Mass 1'); hold on;
plot(t, X(:,2), 'g', 'DisplayName', 'Mass 2');
plot(t, X(:,3), 'b', 'DisplayName', 'Mass 3');
legend;
xlabel('Time (s)');
11/1/24 2:54 PM MATLAB Command Window 2 of 2
ylabel('Displacement (m)');
title('Response of Each Mass Over Time');
% Function to define the system dynamics for ode45
function dXdt = system_dynamics(~, X, M, K)
n = length(X)/2; % Number of degrees of freedom
dXdt = zeros(2*n, 1); % Initialize derivative vector
% Split the state vector into displacements and velocities
displacements = X(1:n);
velocities = X(n+1:end);
% First half of dXdt is velocities
dXdt(1:n) = velocities;
% Second half of dXdt is accelerations
dXdt(n+1:end) = -M \ (K * displacements);
end
