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=(C1+K1) *eps*cos(Omeg*t)+K2* 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’)% dynamic analysis using direct integration method
% % the porbleme is: Mx”+Cx’+Kx=(C1+K1) *eps*cos(Omeg*t)+K2* 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];

tt=0;

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’) % dynamic analysis using direct integration method
% % the porbleme is: Mx”+Cx’+Kx=(C1+K1) *eps*cos(Omeg*t)+K2* 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];

tt=0;

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

Cart

Cart