Problem based on Modal period & Newmark-Beta Method

I’ve generated an MATLAB code based on these pointed problems, such that:
Q: Consider a 2D 3-story shear building model.
1. Form the mass and stiffness matrices for this building model. No damping in needed for this problem.
2. Using trial and errors, select meaningful story mass and stiffness parameters so that the first natural frequency of the building is 1 sec. Consider that the story mass and stiffness values are similar among the three stories; therefore, you need to set only one stiffness and one mass parameter. Use your engineering judgement and explain your logic. Find and plot the mode shapes of the building model.
3. Implement a Newmark-beta algorithm for time integration of the dynamic equilibrium equation. Use a unidirectional earthquake base excitation as input and solve for the acceleration response at the three DOFs. For the earthquake base excitation, use this input file with a sampling frequency of 100 Hz. Plot the time history of absolute acceleration response at each floor.
4. Resample the input signal at 20, 10, and 2 Hz and redo steps 4-5 with integration time step size of 0.05, 0.1, and 0.5 sec, respectively. Compare the results and explain your observations.

MATLAB CODE
—————————————————————————————————————–
clear all;
clc;
close all;
%% INPUT Parameters:
m = 1; % Initial arbitrary mass of each story
k = 1; % Initial arbitrary stiffness of each story
%% Solution-1: Form Mass and Stiffness Matrices
M = m * eye(3); % Mass matrix
K = [ 2*k, -k, 0; -k, 2*k, -k; 0, -k, k]; % Stiffness matrix
%% Solution-2: Adjust Mass and Stiffness based on target frequency
target_frequency = 1; % Natural Frequency (Hz)
target_omega = 2 * pi * target_frequency; % Angular Frequency (rad/s)
% Adjust stiffness based on the target frequency
k = target_omega^2 * m / 2;
K = [ 2*k, -k, 0; -k, 2*k, -k; 0, -k, k]; % Updated Stiffness matrix
%% Solution-3: Find and Plot Mode Shapes
[eigVectors, eigValues] = eig(K, M); % Solve eigenvalue problem
omega = sqrt(diag(eigValues)); % Natural frequencies (rad/s)
figure;
for i = 1:3
subplot(3, 1, i);
plot([0 eigVectors(:,i)’], ‘o-‘); % Plot mode shapes
title([‘Mode Shape ‘, num2str(i)]);
xlabel(‘Story’);
ylabel(‘Displacement’);
end
%% Solution-4: Implement Newmark-beta algorithm for dynamic analysis
% Parameters for Newmark-beta
beta = 1/4;
gamma = 1/2;
% Load earthquake data
data = load(‘CNP100.txt’); % Load earthquake data
dt = 0.01; % Sampling interval for 100 Hz
time = (0:length(data)-1)*dt;
% Initial conditions
u = zeros(3, 1); % Initial displacement
v = zeros(3, 1); % Initial velocity
a = M ((data(1) * ones(3, 1)) – K*u); % Initial acceleration
% Time integration using Newmark-beta method
u_history = zeros(3, length(time));
a_history = zeros(3, length(time));
a_history(:, 1) = a;
for i = 1:length(time)-1
% Predictors
u_pred = u + (dt * v) + ((dt^2) * (0.5-beta) * a);
v_pred = v + (dt * (1-gamma) * a);

% Solve for acceleration
a_new = M ((data(i+1) * ones(3, 1)) – (K * u_pred));

% Correctors
u = u_pred + (beta * (dt^2) * a_new);
v = v_pred + (gamma * dt * a_new);

% Store response
u_history(:, i+1) = u;
a_history(:, i+1) = a_new;
a = a_new;
end
%% Solution-5: Plot time history of absolute acceleration response
figure;
for i = 1:3
subplot(3, 1, i);
plot(time, a_history(i, :));
title([‘Absolute Acceleration Response Sampled at 100 Hz at DOF ‘, num2str(i)]);
xlabel(‘Time (s)’);
ylabel(‘Acceleration (m/s^2)’);
end
%% Solution-6: Resample input signal and repeat the analysis
frequencies = [20, 10, 2];
time_steps = [0.05, 0.1, 0.5];
for j = 1:length(frequencies)
new_fs = frequencies(j);
new_dt = 1/new_fs;
new_data = resample(data, new_fs, 100);
new_time = (0:length(new_data)-1)*new_dt;

% Repeat the Newmark-beta integration with new data
u = zeros(3, 1); % Initial displacement
v = zeros(3, 1); % Initial velocity
a = M ((new_data(1) * ones(3, 1)) – K*u); % Initial acceleration

u_history = zeros(3, length(new_time));
a_history = zeros(3, length(new_time));
a_history(:, 1) = a;

for i = 1:length(new_time)-1
% Predictors
u_pred = u + (new_dt * v) + ((new_dt^2) * (0.5-beta) * a);
v_pred = v + (new_dt * (1-gamma) * a);

% Solve for acceleration
a_new = M ((new_data(i+1) * ones(3, 1)) – K * u_pred);

% Correctors
u = u_pred + (beta * (new_dt^2) * a_new);
v = v_pred + (gamma * new_dt * a_new);

% Store response
u_history(:, i+1) = u;
a_history(:, i+1) = a_new;
a = a_new;
end

% Plot time history of absolute acceleration response
figure;
for i = 1:3
subplot(3, 1, i);
plot(new_time, a_history(i, :));
title([‘Absolute Acceleration Response at DOF ‘, num2str(i), ‘ (Resampled at ‘, num2str(new_fs), ‘ Hz)’]);
xlabel(‘Time (s)’);
ylabel(‘Acceleration (m/s^2)’);
end
end
——————————————————————————————————————–
I’m getting error in obtaining:
The first modal period as 1 second.
The correct Newmark Beta plots.I’ve generated an MATLAB code based on these pointed problems, such that:
Q: Consider a 2D 3-story shear building model.
1. Form the mass and stiffness matrices for this building model. No damping in needed for this problem.
2. Using trial and errors, select meaningful story mass and stiffness parameters so that the first natural frequency of the building is 1 sec. Consider that the story mass and stiffness values are similar among the three stories; therefore, you need to set only one stiffness and one mass parameter. Use your engineering judgement and explain your logic. Find and plot the mode shapes of the building model.
3. Implement a Newmark-beta algorithm for time integration of the dynamic equilibrium equation. Use a unidirectional earthquake base excitation as input and solve for the acceleration response at the three DOFs. For the earthquake base excitation, use this input file with a sampling frequency of 100 Hz. Plot the time history of absolute acceleration response at each floor.
4. Resample the input signal at 20, 10, and 2 Hz and redo steps 4-5 with integration time step size of 0.05, 0.1, and 0.5 sec, respectively. Compare the results and explain your observations.