The graphs of position, velocity and acceleration are wrong.

Hello Everyone. Following is the code I'm using to calculate the response of coupled MDOF system. But the graphs I am getting is undesirable as it should damp but instead it's overshooting. Please Help me. Thanks.

clear all; close all; clc
tspan=[0 2]; %time in seconds
ic=[0.1 0.2 0.1 -0.1 0 0 0 0];
[t,x]=ode45(@diff,tspan,ic);
xf1=x(:,1); xf2=x(:,2); xf3=x(:,3); xf4=x(:,4);
fx1=x(:,5); fx2=x(:,6); fx3=x(:,7); fx4=x(:,8);
m1=95; m2=945; m3=10; m4=80; %Mass Values in kgs
k1=730000; k2=50000; k3=1000; k4=10000; %Stiffness Values in kg/sec^2
c1=10; c2=8000; c3=1000; c4=750; %Damping Values in kg/m.sec^2
M=[m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4]; %Mass Matrix
K=[-(k1+k2) k2 0 0; k2 -(k2+k3) k3 0; 0 k3 -(k3+k4) k4; 0 0 k4 -k4]; %Stiffness Matrix
C=[-(c1+c2) c2 0 0; c2 -(c2+c3) c3 0; 0 c3 -(c3+c4) c4; 0 0 c4 -c4]; %Damping Matrix
A=[-(k1+k2)/m1 k2/m1 0 0; k2/m2 -(k2+k3)/m2 k3/m2 0; 0 k3/m3 -(k3+k4)/m3 k4/m3; 0 0 k4/m4 -k4/m4]; %Characteristic Matrix
[e,d]=eig(A); %Eigen Values and Vector
w=sqrt(abs(d));  %Natural Frequency Matrix
F=1320; %Input force in Newton
a1=(K(1,1)*x(:,1)+K(1,2)*x(:,2)+C(1,1)*x(:,5)+C(1,2)*x(:,6))/M(1,1); %Acc M1
a2=(K(2,1)*x(:,1)+K(2,2)*x(:,2)+K(2,3)*x(:,3)+C(2,1)*x(:,5)+C(2,2)*x(:,6)+C(2,3)*x(:,7))/M(2,2); %Acc M2
a3=(K(3,2)*x(:,2)+K(3,3)*x(:,3)+K(3,4)*x(:,4)+C(3,2)*x(:,6)+C(3,3)*x(:,7)+C(3,4)*x(:,8))/M(3,3); %Acc M3
a4=(K(4,3)*x(:,3)+K(4,4)*x(:,4)+C(4,3)*x(:,7)+C(4,4)*x(:,8))/M(4,4); %Acc M4
subplot(3,1,1)
hold on;
plot(t,xf1,'r')
plot(t,xf2,'g')
plot(t,xf3,'b')
plot(t,xf4,'black')
ylabel('Displacement (m)')
subplot(3,1,2)
hold on;
plot(t,fx1,'r')
plot(t,fx2,'g')
plot(t,fx3,'b')
plot(t,fx4,'black')
ylabel('Velocity (m/s)')
title('Disp, Vel & Acc wrt T')
subplot(3,1,3)
hold on;
plot(t,a1,'r')
plot(t,a2,'g')
plot(t,a3,'b')
plot(t,a4,'black')
legend('m1','m2','m3','m4')
ylabel('Acceleration (m/s^2)')
xlabel('Time (s)')
clc
disp('Natural frequencies (rad/sec)')
disp(w)
disp('Velocity (m/s)')
disp(max(x(:,5)))
disp(max(x(:,6)))
disp(max(x(:,7)))
disp(max(x(:,8)))
disp('Acceleration (m/s^2)')
disp(max(a1))
disp(max(a2))
disp(max(a3))
disp(max(a4))
function yt = diff(t,x)
    m1=95; m2=945; m3=10; m4=80; %Mass Values in kgs
    k1=730000; k2=50000; k3=1000; k4=10000; %Stiffness Values in kg/sec^2
    c1=10; c2=8000; c3=1000; c4=750; %Damping Values in kg/m.sec^2
    M=[m1 0 0 0; 0 m2 0 0; 0 0 m3 0; 0 0 0 m4]; %Mass Matrix
    K=[-(k1+k2) k2 0 0; k2 +(k2+k3) k3 0; 0 k3 -(k3+k4) k4; 0 0 k4 -k4]; %Stiffness Matrix
    C=[-(c1+c2) c2 0 0; c2 -(c2+c3) c3 0; 0 c3 -(c3+c4) c4; 0 0 c4 -c4]; %Damping Matrix
    A=[-(k1+k2)/m1 k2/m1 0 0; k2/m2 -(k2+k3)/m2 k3/m2 0; 0 k3/m3 -(k3+k4)/m3 k4/m3; 0 0 k4/m4 -k4/m4]; %Characteristic Matrix
    [e,d]=eig(A) %Eigen Values and Vector
    w=sqrt(abs(d))  %Natural Frequency Matrix
    F=1320; %Input force in Newton
    yt=zeros(8,1);
    yt(1)=x(5); %Velocity M1
    yt(2)=x(6); %Velocity M2
    yt(3)=x(7); %Velocity M3
    yt(4)=x(8); %Velocity M4
    yt(5)=(K(1)*x(1)+K(2)*x(2)+C(1)*x(5)+C(2)*x(6))/M(1); %Acc M1
    yt(6)=(K(5)*x(1)+K(6)*x(2)+K(7)*x(3)+C(5)*x(5)+C(6)*x(6)+C(7)*x(7))/M(6); %Acc M2
    yt(7)=(K(10)*x(2)+K(11)*x(3)+K(12)*x(4)+C(10)*x(6)+C(11)*x(7)+C(12)*x(8))/M(11); %Acc M3
    yt(8)=(K(15)*x(3)+K(16)*x(4)+C(15)*x(7)+C(16)*x(8))/M(16); %Acc M4
end

[darova]

please show original equations

[Ankit Joshi]

@darova I've attached the image with original question. Please find attached.

Answer 1

[Mathieu NOE]

hello

I didn't find any mistake in your original equations , but I am lazy and I found your code quite long and complicated... why not make simply a state space model and use all available tools to run a time domain simulation (using either step , impulse or lsim ( lsim Simulate time response of dynamic systems to arbitrary inputs.))

this is my code that I double checked and seems to work as expected. You can easily add the eigenvalues computation of the A matrix if you want

also here I plotted only the acceleration responses of the 4 masses for each step input

good luck

clearvars; close all; clc

M1=95; M2=945; M3=10; M4=80; %Mass Values in kgs
K1=730000; K2=50000; K3=1000; K4=10000; %Stiffness Values in kg/sec^2
C1=10; C2=8000; C3=1000; C4=750; %Damping Values in kg/m.sec^2

[A,B,C,D] = model_4masses(M1,M2,M3,M4,C1,C2,C3,C4,K1,K2,K3,K4);

% 4 plots showing a step response on each mass (separetely)
% plotted data : acceleration 
for ci = 1:4
    figure(ci)
    step(A,B,C,D,ci)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,B,C,D] = model_4masses(M1,M2,M3,M4,C1,C2,C3,C4,K1,K2,K3,K4)

% state space model of the stacked 4 masses - springs - dash pots

A = [0  1  0  0  0  0  0  0; 
 -(K1+K2)/M1  -(C1+C2)/M1  K2/M1  C2/M1  0  0  0  0 ;
 0  0  0  1  0  0  0  0;
 K2/M2  C2/M2  (-K3-K2)/M2  (-C3-C2)/M2  K3/M2  C3/M2  0  0 ;
 0  0  0  0  0  1  0  0 ;
 0  0  K3/M3  C3/M3  (-K4-K3)/M3  (-C4-C3)/M3  K4/M3  C4/M3 ;
 0  0  0  0  0  0  0  1;
 0  0  0  0  K4/M4  C4/M4  -K4/M4  -C4/M4];


B = [0  0  0  0  0  0 ;
 1/M1  0  0  0  0  0 ;
 0  0  0  0  0  0 ;
 0  1/M2  0  0  0  0 ;
 0  0  0  0  0  0 ;
 0  0  1/M3  0  0  0 ;
 0  0  0  0  0  0 ;
 0  0  0  1/M4  0  0];

C = [-(K1+K2)/M1  -(C1+C2)/M1  K2/M1  C2/M1  0  0  0  0 
 K2/M2  C2/M2  (-K3-K2)/M2  (-C3-C2)/M2  K3/M2  C3/M2  0  0 
 0  0  K3/M3  C3/M3  (-K4-K3)/M3  (-C4-C3)/M3  K4/M3  C4/M3 
 0  0  0  0  K4/M4  C4/M4  -K4/M4  -C4/M4];

D = [1/M1  0  0  0  0  0 ;
 0  1/M2  0  0  0  0 ;
 0  0  1/M3  0  0  0 ;
 0  0  0  1/M4  0  0];

end