Hello everyone i am solving multi degree freedom of system by using Newmarks beta method in which loading function F(x,v). When i am ruing the code, the result i

% dynamic analysis using direct integration method

% % the porbleme is: Mx''+Cx'+Kx=(C1x'+K1x) *eps*cos(Omeg*t)+K2x* eps^2cos(Omeg*t)^2

clc;

clear all

format short e

%close all

m=[18438.6 0 0;0 13761 0;0 0 9174];

disp(' mass matrix')

m;

[ns,ms]=size(m);

% if forces are acting at degrees of freedom

m=[40000 0 0;0 20000 0;0 0 20000];

c0=[0 0 0;0 -8000 8000;0 -8000 8000];

c1=[0 0 0;0 -4000 8000;0 -8000 4000];

k0=[30000 -10000 0;-10000 20000 -10000;0 -10000 10000];

k1=[30000 -10000 0;-10000 50000 -10000;0 -10000 50000];

k2=[30000 -10000 0;-10000 20000 -40000;0 -40000 10000];

% disp(' force at various degrees of freedom');

% f;

%if base ground acceleration is given

% dis='disp.dat'

% di=load(dis)

% % convert to equivalent nodal loads

% for i=1:ns

% f(:,i)=-di*m(i,i)

% end

disp(' damping matrix')

c0;

disp(' stiffness matrix')

k0;

format long;

kim=inv(k0)*m;

[evec,ev]=eig(kim);

for i=1:ns

omega(i)=1/sqrt(ev(i,i));

end

disp(' natural frequencies')

omega;

% give gamma=0.5 and beta=0.25 for Newmark average accln method

%gama=0.5;

%beta=0.25;

%give gamma=0.5 and beta=0.1667 for Newmark linear accln method

gama=0.5;

beta=0.167;

%give initial conditions for displacements

u0=[0 0.01 0.05];

disp(' initial displacements')

u0;

%give initial condition for velocities

v0=[0. 0. 0.];

%y0=[0;0.01;0.05;0;0;0];

disp(' initial velocities')

v0;

om=5; eps=0.01;

IM=inv(m);

X=u0'; Z=v0';

tt=0;

% S1=-IM*k0+IM*k1*eps.*cos(om*tt)+IM*k2*eps^2.*cos(om*tt)^2;

% S2=-IM*c0+IM*c1*eps.*cos(om*tt);

dt=0.02;

S=k2*dt^2*beta*eps^2*cos(om*tt)^2+k1*dt^2*beta*eps*cos(om*tt)+...

c1*gama*dt*eps*cos(om*tt)+k0*dt^2*beta+c0*gama*dt+m;

S1=k1*eps.*cos(om*tt)+k2*eps^2.*cos(om*tt)^2;

S2=c1*eps.*cos(om*tt);

f(1,:)=(S1*X+S2*Z);

%for i=1:ns

a0=-inv(m)*(f(1,:)'+c0*v0'+k0*u0');

%end

kba=k0+(gama/(beta*dt))*c0+(1/(beta*dt*dt))*m;

kin=inv(kba);

aa=(1/(beta*dt))*m+(gama/beta)*c0;

bb=(1/(2.0*beta))*m+dt*(gama/(2.0*beta)-1)*c0;

u(1,:)=u0;

v(1,:)=v0;

a(1,:)=a0;

t=linspace(0,5,251);

for i=2:10%251

X=u(i-1,:)'; Z=v(i-1,:)';

S1=k1*eps.*cos(om*tt)+k2*eps^2.*cos(om*tt)^2;

S2=c1*eps.*cos(om*tt);

f(i,:)=(S1*X+S2*Z); %%%%% ??????

df=f(i,:)-f(i-1,:);

dfb(i,:)=df+v(i-1,:)*aa'+a(i-1,:)*bb';

du(i,:)=dfb(i,:)*kin;

dv(i,:)=(gama/(beta*dt))*du(i,:)-(gama/beta)*v(i-1,:)+dt*...

(1-gama/(2.0*beta))*a(i-1,:);

da(i,:)=(1/(beta*dt^2))*du(i,:)-(1/(beta*dt))*v(i-1,:)-(1/(2.0*beta))*a(i-1,:);

u(i,:)=u(i-1,:)+du(i,:);

v(i,:)=v(i-1,:)+dv(i,:);

a(i,:)=a(i-1,:)+da(i,:);

end

%figure(1);

hold on

plot(tt,u(:,1),'k');

xlabel(' time in secs');

ylabel(' roof displacement');

title(' displacement response of the roof');

%figure(2);

hold on

plot(tt,u(:,2),'k');

xlabel(' time in secs');

ylabel(' roof velocity');

title('velocity response of the roof');

%figure(3);

hold on

plot(tt,u(:,3),'k');

xlabel(' time in secs');

ylabel(' roof acceleration');

title(' acceleration response of the roof')

[Ashok]

Hi abdellhadi,

Your question seems to be incomplete. I observed that the figure is blank. Is that the error you are referring to?

Answer 1

[Zinea]

Hi abdellhadi,

The code needs the following changes to give accurate results for the Newmark-beta method for solving the dynamic response of a multi-degree-of-freedom system with a specific loading function:

1. The iteration loop goes from ‘i’ = 2 to 10, but the time vector ‘t’ is defined from 0 to 5 seconds with 251 points. The loop needs to be adjusted to start from ‘i’ = 2 to 251 to cover the entire time range.
2. ‘tt’ must be updated within the loop to reflect the current time step: tt = t(i)
3. The plotting commands are currently inside the loop and may not properly display the results. They should be moved outside the loop along with the usage of separate figures for displacement, velocity and acceleration.
4. It would be good to use a vectorized approach to compute results for multiple degrees of freedom simultaneously.
5. Pre-allocating arrays (u,v,a) for storing results and initializing them with zeros or initial conditions in a single step.

Here is the MATLAB code for your reference:

% Define mass, damping, and stiffness matrices 
m = [40000 0 0; 0 20000 0; 0 0 20000]; 
c0 = [0 0 0; 0 -8000 8000; 0 -8000 8000]; 
c1 = [0 0 0; 0 -4000 8000; 0 -8000 4000]; 
k0 = [30000 -10000 0; -10000 20000 -10000; 0 -10000 10000]; 
k1 = [30000 -10000 0; -10000 50000 -10000; 0 -10000 50000]; 
k2 = [30000 -10000 0; -10000 20000 -40000; 0 -40000 10000]; 

% Calculate natural frequencies 
kim = inv(k0) * m; 
[evec, ev] = eig(kim); 
omega = 1 ./ sqrt(diag(ev))'; % Vectorized calculation of natural frequencies 
disp('Natural frequencies:'); 
disp(omega); 

% Newmark-beta method parameters 
gama = 0.5; 
beta = 0.167; 

% Initial conditions 
u0 = [0 0.01 0.05]; 
v0 = [0 0 0]; 

% System parameters 
om = 5;  
eps = 0.01; 
IM = inv(m); 

% Time settings 
dt = 0.02; 
t = linspace(0, 5, 251); 

% Pre-allocate response arrays 
u = zeros(length(t), length(u0)); 
v = zeros(length(t), length(v0)); 
a = zeros(length(t), length(u0)); 

% Initial conditions 
u(1, :) = u0; 
v(1, :) = v0; 

% Initial acceleration 
a(1, :) = (-IM * (c0 * v0' + k0 * u0'))'; 

% Calculate effective stiffness matrix 
kba = k0 + (gama / (beta * dt)) * c0 + (1 / (beta * dt^2)) * m; 
kin = inv(kba); 
aa = (1 / (beta * dt)) * m + (gama / beta) * c0; 
bb = (1 / (2.0 * beta)) * m + dt * (gama / (2.0 * beta) - 1) * c0; 

% Precompute constant matrices 
S1_base = eps * k1; 
S2_base = eps * c1; 
S1_quad = eps^2 * k2; 

% Time integration loop 
for i = 2:length(t) 
    tt = t(i); 
    
    % Update force vector using vectorized operations 
    cos_om_tt = cos(om * tt); 
    S1 = S1_base * cos_om_tt + S1_quad * cos_om_tt^2; 
    S2 = S2_base * cos_om_tt; 
    
    % Calculate effective force increment 
    f = S1 * u(i-1, :)' + S2 * v(i-1, :)'; 
    dfb = f + aa * v(i-1, :)' + bb * a(i-1, :)'; 
    
    % Newmark-beta updates, vectorized 
    du = kin * dfb; 
    dv = (gama / (beta * dt)) * du - (gama / beta) * v(i-1, :)' + dt * (1 - gama / (2.0 * beta)) * a(i-1, :)'; 
    da = (1 / (beta * dt^2)) * du - (1 / (beta * dt)) * v(i-1, :)' - (1 / (2.0 * beta)) * a(i-1, :)'; 
    
    % Update responses 
    u(i, :) = u(i-1, :) + du';  % Transpose du to match 1-by-3 
    v(i, :) = v(i-1, :) + dv'; 
    a(i, :) = a(i-1, :) + da'; 
end 

% Plotting results for displacement (moved outside loop) 
figure; 
plot(t, u(:, 1), 'k'); 
xlabel('Time (s)'); 
ylabel('Displacement'); 
title('Displacement Response of the First Degree of Freedom'); 
grid on 
% Simlarly plot for velocity and acceleration 

Here is the output from the above code and the plotted figure:

You may refer to the following documentation links for more information on vectorization and pre-allocation in MATLAB :

1. https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html
2. https://www.mathworks.com/help/matlab/matlab_prog/preallocating-arrays.html

Hope this helps!

[abdellhadi]

Thank you so much