% 3DOF_shear_building_modal_analysis.
m
% Fully-ready MATLAB script for modal analysis of the 3-DOF shear building
% Computes natural frequencies, mode shapes, steady-state harmonic
response,
% impulse response, and includes damping (modal damping ratio).
% Plots modal coordinates and physical DOF responses.
clear; close all; clc;
%% Parameters
m1 = 5000; m2 = 5000; m3 = 5000; % kg
k1 = 2000; k2 = 2000; k3 = 2000; % N/m
M = diag([m1 m2 m3]);
K = [k1+k2, -k2, 0;
-k2, k2+k3, -k3;
0, -k3, k3];
% External force vector (applied at DOF1)
F0 = [300; 0; 0]; % amplitude for harmonic forcing (N)
%% Eigenvalue problem (undamped)
[Phi, D] = eig(K, M); % generalized eigenproblem K*phi =
lambda*M*phi
omega_n2 = diag(D);
omega_n = sqrt(omega_n2);
[omega_n, idx] = sort(omega_n);
Phi = Phi(:, idx); % reorder modes by ascending frequency
omega_n = omega_n(:);
% Mass-normalize modes: Phi^T * M * Phi = I
for j = 1:size(Phi,2)
scale = sqrt(Phi(:,j)' * M * Phi(:,j));
Phi(:,j) = Phi(:,j) / scale;
end
% Verify normalization
MTMP = Phi' * M * Phi; % should be identity
KTKP = Phi' * K * Phi; % should be diagonal with omega_n.^2
%% Print natural frequencies
fprintf('Natural frequencies (rad/s):\n');
for j=1:3
fprintf(' mode %d: %8.4f rad/s (f = %6.4f Hz)\n', j, sqrt(KTKP(j,j)),
sqrt(KTKP(j,j))/(2*pi));
end
%% Modal forcing (mass-normalized)
f_modal = Phi' * F0; % modal force amplitudes (for harmonic forcing)
%% (1) Steady-state harmonic response (undamped) for F(t) = 300*cos(20 t)
omega = 20; % forcing frequency (rad/s)
% Modal steady-state amplitude (undamped): q_j = f_j / (omega_n^2 -
omega^2) * cos(omega t)
Q_modal_amp = f_modal ./ (diag(KTKP) - omega^2);
% Reconstruct physical steady-state amplitude (cosine-phase)
x_amp = Phi * Q_modal_amp; % amplitudes for each DOF (cosine term)
% Time vector for plotting
T = 4*pi/omega; % show two full forcing periods-ish
t = linspace(0, T, 2001);
x_t_undamped = zeros(3, numel(t));
for k=1:numel(t)
x_t_undamped(:,k) = x_amp * cos(omega * t(k));
end
%% (2) Impulse response: F(t) = 300 * delta(t)
% For impulse of magnitude F_imp (N*s) applied at t=0, jump in velocity
x_dot(0+) = M^{-1} * F_imp
F_impulse = 300; % treat as impulse magnitude (N*s)
x_dot_0 = M \ (F_impulse * [1;0;0]); % initial velocity vector
q_dot_0 = Phi' * (M * x_dot_0); % modal initial velocities (since
q_dot = Phi' * M * x_dot)
% Modal free response for each mode: q_j(t) = (q_dot_j(0+)/omega_n_j) *
sin(omega_n_j * t)
x_t_impulse = zeros(3, numel(t));
for j=1:3
qj = (q_dot_0(j) / sqrt(KTKP(j,j))) * sin(sqrt(KTKP(j,j)) * t);
x_t_impulse = x_t_impulse + Phi(:,j) * qj;
end
%% (3) Include small damping (modal damping ratio zeta = 0.02 for all
modes)
zeta = 0.02; % modal damping ratio (same for all modes)
% Steady-state harmonic response with damping (complex amplitudes)
omega_n_diag = sqrt(KTKP);
H_modal = 1 ./ (omega_n_diag - omega^2 + 2i * zeta * sqrt(omega_n_diag) *
omega); % Note: denom uses omega_n^2 form
% Correction: omega_n_diag currently stores omega_n^2; let's recompute
properly
omega_n_sq = diag(KTKP);
omega_n_vals = sqrt(omega_n_sq);
H_modal = 1 ./ (omega_n_sq - omega^2 + 2i * zeta * omega_n_vals * omega);
Q_modal_complex = f_modal .* H_modal; % complex modal amplitudes
% Physical complex response: x(t) = Re{ Phi * (Q_modal_complex .* exp(i
omega t)) }
x_t_damped = zeros(3, numel(t));
for k=1:numel(t)
x_complex = Phi * (Q_modal_complex .* exp(1i * omega * t(k)));
x_t_damped(:,k) = real(x_complex);
end
%% (4) Damped impulse response (modal damping applied)
% Modal initial velocities are the same as undamped (impulse instant).
Solve damped modal eqns:
% q_j'' + 2*zeta*omega_n_j*q_j' + omega_n_j^2 * q_j = 0 with q_dot(0+) and
q(0)=0
x_t_impulse_damped = zeros(3, numel(t));
for j=1:3
wn = omega_n_vals(j);
wd = wn * sqrt(1 - zeta^2);
qj = (q_dot_0(j) / wd) * exp(-zeta * wn * t) .* sin(wd * t); % damped
modal response
x_t_impulse_damped = x_t_impulse_damped + Phi(:,j) * qj;
end
%% Plots
figure('Name','Steady-state harmonic response (undamped vs damped)');
subplot(3,1,1);
plot(t, x_t_undamped(1,:), 'LineWidth', 1.2); hold on;
plot(t, x_t_damped(1,:), '--', 'LineWidth', 1.2); grid on;
title('DOF 1 response: Undamped (solid) vs Damped (dashed)'); xlabel('t
(s)'); ylabel('x_1 (m)');
legend('Undamped','Damped');
subplot(3,1,2);
plot(t, x_t_undamped(2,:), 'LineWidth', 1.2); hold on;
plot(t, x_t_damped(2,:), '--', 'LineWidth', 1.2); grid on;
title('DOF 2 response'); xlabel('t (s)'); ylabel('x_2 (m)');
subplot(3,1,3);
plot(t, x_t_undamped(3,:), 'LineWidth', 1.2); hold on;
plot(t, x_t_damped(3,:), '--', 'LineWidth', 1.2); grid on;
title('DOF 3 response'); xlabel('t (s)'); ylabel('x_3 (m)');
% Impulse response plots
figure('Name','Impulse response (undamped vs damped)');
for i=1:3
subplot(3,1,i);
plot(t, x_t_impulse(i,:), 'LineWidth', 1.2); hold on;
plot(t, x_t_impulse_damped(i,:), '--', 'LineWidth', 1.2); grid on;
title(sprintf('DOF %d impulse response', i)); xlabel('t (s)');
ylabel('x (m)');
legend('Undamped','Damped');
end
%% Modal coordinates plotting
q_modal_undamped = zeros(3, numel(t));
for j=1:3
q_modal_undamped(j,:) = (f_modal(j) ./ (omega_n_sq(j) - omega^2)) *
cos(omega * t);
end
figure('Name','Modal coordinates (undamped steady-state)');
for j=1:3
subplot(3,1,j);
plot(t, q_modal_undamped(j,:), 'LineWidth', 1.2); grid on;
title(sprintf('Modal coord q_%d (undamped steady-state)', j));
xlabel('t (s)'); ylabel(sprintf('q_%d', j));
end
%% End of script
fprintf('\nScript complete.\n');.
Output code is
Natural frequencies (rad/s):
mode 1: 0.2815 rad/s (f = 0.0448 Hz)
mode 2: 0.7887 rad/s (f = 0.1255 Hz)
mode 3: 1.1396 rad/s (f = 0.1814 Hz)