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