MATLAB: Plotting Displacement, Velocity and Acceleration vs. Time

#earthquake analysis

Hi, I am a new MATLAB user, so please bear with me. I'm trying to plot the values of displacement(u), velocity(v), and acceleration(a) versus time(tn) using the Newmark's Beta Method for a SDOF system. But all I get is an empty plot. I would really appreciate if someone could help me. Thank you in advance. 
%%
%Newmark-Beta Method
%% Given Values
clear all; 
clc; 
close all;
load elcentro.dat           % load the Seismic data
t = elcentro(:,1);          %time values
Ag = elcentro(:,2);         %acceleration values
g = 9.81;                   %acceleration due to gravity
Dt = 0.02;                  %time step
m = 1;                      %mass
xi = 0.02;                  %damping ratio
endp = 31.18;               %terminating value of t
tn = 0:Dt:31.18;
T = 1;                      %Period
%% Solver
    omega = 2*pi/T;
    k = (omega^2)*m;     %stiffness
    c = 2*m*omega*xi;    %damping coefficient
    
    %parameters
    gamma = 1/2;
    beta = 1/6;
    
    u(1) = 0;               %displacement initial condition
    v(1) = 0;               %velocity initial condition
    a(1) = 0;               %acceleration initial condition
    
    keff = k +gamma*c/(beta*Dt) + m/(beta*Dt*Dt);
    A = m/(beta*Dt)+gamma*c/beta;
    B = m/(2*beta)+(Dt*c*((0.5*gamma/beta)-1));
    
    u0 = u;
    v0 = v;
    a0 = a;
    
    %loop for every time step
    for i = 1:(length(t)-Dt)
        DF = -m*g.*(Ag(i+1)-Ag(i));
            [t,u,v,a] = Newmark(t,A,B,DF,Dt,keff,u0,v0,a0,gamma,beta);
        u_t1(:,i) = u(:,1)';
        v_t1(:,i) = v(:,1)';
        a_t1(:,i) = a(:,1)';
            
        u0 = u;
        v0 = v;
        a0 = a;
        t0 = t;
    end
    
    umax = max(abs(u_t1));
    vmax = max(abs(v_t1));
    amax = max(abs(a_t1));
    fprintf('\n-------NEWMARK BETA METHOD-------\n');
    fprintf('\nMaximum Response of a SDOF System\n\n');
    fprintf('Maximum Displacement = %f\n',umax);
    fprintf('Maximum Velocity     = %f\n',vmax);
    fprintf('Maximum Acceleration = %f\n',amax);
    
%% Plotting the Figures for Response
figure(1)
plot(tn,u);
title('Displacement Response')
xlabel('Period(sec)');
ylabel('Displacement (m)');
% Save picture


saveas(gca,'Displacement2.tif');
figure(2)
plot(tn,v);
title('Velocity Response')
xlabel('Period(sec)');
ylabel('Velocity (m/s)');
% Save picture
saveas(gca,'Velocity2.tif');
figure(3)
plot(tn,a);
title('Acceleration Response')
xlabel('Period(sec)');
ylabel('Acceleration (m/s^2)');
% Save picture
saveas(gca,'Acceleration2.tif');
%% Function for Newmark-Beta
function [t,u,v,a] = Newmark (t,A,B,DF,Dt,keff,u0,v0,a0,gamma,beta)
    DFbar = DF + A*v0+B*a0;
    Du = DFbar/keff;
    Dudot = gamma*Du/(beta*Dt)-((gamma*v0)/beta)+ Dt*a0*(1-0.5*(gamma/beta));
    Dudotdot = Du/(beta*Dt*Dt)-v0/(beta*Dt)-a0/(2*beta);
    
    u = u0+Du;
    v = v0+Dudot;
    a = a0+Dudotdot;
    t = t+Dt; 
end

Best Answer

function Myfunction()
%% @ Jeanette M. Sarra, CE 165
%Newmark-Beta Method
%% Given Values
clear all;
clc;
close all;
elcentro = load('elcentro.dat') ;           % load the Seismic data
t = elcentro(:,1);          %time values
Ag = elcentro(:,2);         %acceleration values
g = 9.81;                   %acceleration due to gravity
Dt = 0.02;                  %time step
m = 1;                      %mass
xi = 0.02;                  %damping ratio
endp = 31.18;               %terminating value of t
tn = 0:Dt:31.18;
T = 1;                      %Period
%% Solver
omega = 2*pi/T;
k = (omega^2)*m;     %stiffness
c = 2*m*omega*xi;    %damping coefficient
%parameters
gamma = 1/2;
beta = 1/6;
u(1) = 0;               %displacement initial condition
v(1) = 0;               %velocity initial condition
a(1) = 0;               %acceleration initial condition
keff = k +gamma*c/(beta*Dt) + m/(beta*Dt*Dt);
A = m/(beta*Dt)+gamma*c/beta;
B = m/(2*beta)+(Dt*c*((0.5*gamma/beta)-1));
u0 = u;
v0 = v;
a0 = a;
%loop for every time step
u_t1 = zeros(1,length(t)) ; u_t1(1) = u0 ; 
v_t1 = zeros(1,length(t)) ; v_t1(1) = v0 ; 
a_t1 = zeros(1,length(t)) ; a_t1(1) = a0 ; 
for i = 2:(length(t)-Dt)
    DF = -m*g.*(Ag(i+1)-Ag(i));
    [t,u,v,a] = Newmark(t,A,B,DF,Dt,keff,u0,v0,a0,gamma,beta);
    u_t1(:,i) = u(:,1)';
    v_t1(:,i) = v(:,1)';
    a_t1(:,i) = a(:,1)';
    
    u0 = u;
    v0 = v;
    a0 = a;
    t0 = t;
end
umax = max(abs(u_t1));
vmax = max(abs(v_t1));
amax = max(abs(a_t1));
fprintf('\n-------NEWMARK BETA METHOD-------\n');
fprintf('\nMaximum Response of a SDOF System\n\n');
fprintf('Maximum Displacement = %f\n',umax);
fprintf('Maximum Velocity     = %f\n',vmax);
fprintf('Maximum Acceleration = %f\n',amax);
%% Plotting the Figures for Response
figure(1)
plot(tn,u_t1);
title('Displacement Response')
xlabel('Period(sec)');
ylabel('Displacement (m)');
% Save picture


saveas(gca,'Displacement2.tif');
figure(2)
plot(tn,v_t1);
title('Velocity Response')
xlabel('Period(sec)');
ylabel('Velocity (m/s)');
% Save picture
saveas(gca,'Velocity2.tif');
figure(3)
plot(tn,a_t1);
title('Acceleration Response')
xlabel('Period(sec)');
ylabel('Acceleration (m/s^2)');
% Save picture
saveas(gca,'Acceleration2.tif');
end
%% Function for Newmark-Beta
function [t,u,v,a] = Newmark (t,A,B,DF,Dt,keff,u0,v0,a0,gamma,beta)
DFbar = DF + A*v0+B*a0;
Du = DFbar/keff;
Dudot = gamma*Du/(beta*Dt)-((gamma*v0)/beta)+ Dt*a0*(1-0.5*(gamma/beta));
Dudotdot = Du/(beta*Dt*Dt)-v0/(beta*Dt)-a0/(2*beta);
u = u0+Du;
v = v0+Dudot;
a = a0+Dudotdot;
t = t+Dt;
end
Related Question