Adjust the given code based on Newmark Beta Method based on input base acceleration file.

clear all;

clc;

close all;

%% Solution 1 & 2:

% Define parameters

m = 1; % Mass of each story (kg)

k = 1; % Stiffness of each story (N/m)

% Mass matrix

M = m * eye(3);

% Stiffness matrix

K = [2*k, -k, 0; -k, 2*k, -k; 0, -k, k];

% Solve the eigenvalue problem

[eigenvectors, eigenvalues] = eig(K, M);

% Extract eigenvalues

omega_squared = diag(eigenvalues);

% Convert to natural frequencies (rad/s)

natural_frequencies = sqrt(omega_squared);

% Convert to natural frequencies (Hz)

natural_frequencies_Hz = natural_frequencies / (2 * pi);

% Display results

disp('Natural Frequencies (Hz):');

disp(natural_frequencies_Hz);

% Adjust parameters to get the first natural frequency to 1 Hz

% This requires iterating and adjusting m and k until the first natural frequency is 1 Hz

desired_frequency = 1; % Hz

desired_omega = desired_frequency * 2 * pi; % rad/s

desired_omega_squared = desired_omega^2;

tolerance = 1e-3; % Tolerance for convergence

max_iterations = 1000; % Maximum number of iterations

iteration = 0;

while iteration < max_iterations

iteration = iteration + 1;

% Update stiffness matrix with current stiffness

K = [2*k, -k, 0; -k, 2*k, -k; 0, -k, k];

% Solve the eigenvalue problem again

[eigenvectors, eigenvalues] = eig(K, M);

% Extract the first eigenvalue

first_omega_squared = eigenvalues(1, 1);

% Adjust mass and stiffness

k = k * (desired_omega_squared / first_omega_squared);

M = m * eye(3);

% Check convergence

if abs(first_omega_squared - desired_omega_squared) < tolerance

break;

end

end

% Display final parameters

disp('Final Parameters:');

disp(['Mass per story (kg): ', num2str(m)]);

disp(['Stiffness per story (N/m): ', num2str(k)]);

% Display final natural frequencies

natural_frequencies = sqrt(diag(eigenvalues)) / (2 * pi);

disp('Final Natural Frequencies (Hz):');

disp(natural_frequencies);

% Normalize mode shapes for plotting

mode_shapes = eigenvectors ./ max(abs(eigenvectors), [], 1);

% Plot mode shapes

figure;

title('Mode Shapes of the 3-Story Shear Building');

%% Solution-3: Mode Shapes Plotting:

stories = [0; 1; 2; 3]; % Include ground level

for i = 1:3

subplot(1, 3, i); % Change to a 1x3 grid

plot([0; mode_shapes(:,i)], stories, '-o');

title(['Mode Shape ', num2str(i)]);

xlabel('Amplitude');

ylabel('Story');

grid on;

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