Prepare an MATLAB code based on Newmark Beta Algorithm problem

The problem is attached herewith along with needed text file.

[Steven Lord]

This sounds like a homework assignment. If it is, show us the code you've written to try to solve the problem and ask a specific question about where you're having difficulty and we may be able to provide some guidance.

If you aren't sure where to start because you're not familiar with how to write MATLAB code, I suggest you start with the free MATLAB Onramp tutorial to quickly learn the essentials of MATLAB.

If you aren't sure where to start because you're not familiar with the mathematics you'll need to solve the problem, I recommend asking your professor and/or teaching assistant for help.

[Parvesh Deepan]

I've written this code.

%%%%%%%%%%%%%%%%%%%%%%%%% Assignment Problem %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;

clc;

close all;

%% INPUT Parameters:

m = 1; % mass of each story (arbitrary value, to be adjusted)

k = 1; % stiffness of each story (arbitrary value, to be adjusted)

%% Solution-1:

M = m * eye(3); % Mass matrix

K = [ 2*k, -k, 0;-k, 2*k, -k;0, -k, k]; % Stiffness matrix

%% Solution-2:

target_frequency = 1; % Natural Frequency (Hz)

target_omega = 2 * pi * target_frequency; % Angular Frequency (rad/s)

k = target_omega^2 * m / 2; % Mass and stiffness Adjustment Based on the target frequency

K = [ 2*k, -k, 0;-k, 2*k, -k;0, -k, k]; % Stiffness matrix Updation with the new k value

%% Solution-3:

[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:

% Parameters for Newmark-beta

beta = 1/4;

gamma = 1/2;

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);

v = zeros(3, 1);

a = M \ (data(1) * ones(3, 1)); % initial acceleration

% Time integration

u_history = zeros(3, length(time));

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;

end

%% Solution-5: Plot time history of absolute acceleration response --

figure;

for i = 1:3

subplot(3, 1, i);

plot(time, u_history(i, :));

title(['Absolute Acceleration Response at DOF ', num2str(i)]);

xlabel('Time (s)');

ylabel('Acceleration (m/s^2)');

end

%% Solution-6:

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);

v = zeros(3, 1);

a = M \ (new_data(1) * ones(3, 1)); % initial acceleration

u_history = zeros(3, length(new_time));

for i = 1:length(new_time)-1

u_pred = u + new_dt * v + new_dt^2 * (0.5-beta) * a;

v_pred = v + new_dt * (1-gamma) * a;

a_new = M \ (new_data(i+1) * ones(3, 1) - K * u_pred);

u = u_pred + beta * new_dt^2 * a_new;

v = v_pred + gamma * new_dt * a_new;

u_history(:, i+1) = u;

end

% Plot time history of absolute acceleration response

figure;

for i = 1:3

subplot(3, 1, i);

plot(new_time, u_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

[Parvesh Deepan]

But the time history responses are incorrect.

[Umar]

Hi YQ, The line a = M \ (data(1) * ones(3, 1)); calculates the initial acceleration using the earthquake data at the first time step. However, this calculation should be based on the initial displacement and velocity of the system, not the earthquake data.

By initializing the initial acceleration as zero and calculating it based on the initial displacement and velocity, you can ensure that the time history responses are calculated correctly.

Hope this will help resolve your issue.