Mass spring damper system not working

I'm solving an ode45 for a mass spring damper system as shown below. But the results i'm getting don't oscillate as much as they should. I have seen other questions about this (link) but can not get it to work. All help is apprecriated, thanks.

%% Initialisation
clc;
close all;
clear; 

load("RealWorldData.mat")

m_total =[ 363.5  463.7   363.5 ; 563.9  663.1   563.9 ; 764.3  864.5  764.3]; % Different masses
m = m_total(2,:); % Masses used

t_end = 20; % Time in seconds
w = 5; % Staitionary frequency
k = 40; % Spring constant
b = 0.2; % Friction constant

%% Calculations
t = [0 t_end]; % Time period
x0 = [0, 0, 0, 0, 0, 0, 0]; % Initial conditions

Change = ischange(RealX); 
index = find(Change,1,'first'); % Starting point realdata
index2 = find(Change,1,'last'); % Ending point realdata

RealX = RealX(index:index2,:)*0.001; 
RealT = (1:length(RealX)).*0.001;

[Time,X] = ode45(@(t,x)MassSpringFunction(t,x,m,w,k,b),t,x0);
displacement = X(:,1:3);
velocity = X(:,4:6);

%% Plotting
figure
plot(Time,displacement)
hold on 
plot(RealT,RealX)
legend("RealWorldData","U1" , "U2" , "U3")

%% System
function dxdt = MassSpringFunction(t, x, m, w, k, b)

% Acceleration of 10 hz/s
Freq = 2*pi*10*t;

if w > Freq % Acceleration of frequency
    xin = sin(Freq.*t);
else % Stationairy input
    xin = sin(2*pi*w.*t);
end

dxdt = zeros(6,1);

% Velocities
dxdt(1) = x(4);
dxdt(2) = x(5);
dxdt(3) = x(6);

% Positions
dxdt(4) = -k/m(1) .* (x(1) - xin)  - k/m(1).*(x(1)-x(2)) - b/m(1).*x(4);
dxdt(5) = -k/m(2) .* (x(2) - x(1)) - k/m(2).*(x(2)-x(3)) - b/m(2).*x(5);
dxdt(6) = -k/m(3) .* (x(3) - x(2)) - k/m(3).*(x(3)     ) - b/m(3).*x(6);
end

[Sam Chak]

Hi @Dirk te Brake

I didn't review your model as you didn't provide the derivation steps. However, I made a change by switching the signs from 'plus' to 'minus'. This adjustment ensures that the system's energy is dissipated through negative kinetic energy flow, aligning with the principles of physics.

% Mass
m_total = [363.5  463.7 363.5; 
           563.9  663.1 563.9; 
           764.3  864.5 764.3];
m       = m_total(2,:);

% Staitionary frequency
w       = 5;

% Constants for springs
k       = 80;
b       = 0.2;

t       = [0 40];
x0      = [0, 0, 0, 0, 0, 0];

[Time,x]= ode45(@(t, x) MassSpringFunction(t, x, m, w, k, b), t, x0);

figure(1)
plot(Time,x(:,1))
hold on 
plot(Time,x(:,3))
plot(Time,x(:,5))
legend("U1" , "U2" , "U3")
hold off
grid on
xlabel('t')
ylabel('\bf{U}')

%% System
function dxdt = MassSpringFunction(t, x, m, w, k, b)
    % Acceleration of 10 hz/s
    Freq = 2*pi*10*t;
    
    if w > Freq % Acceleration of frequency
        xin = sin(Freq.*t);
    else % Stationairy input
        xin = sin(2*pi*w.*t);
    end
    
    dxdt    = zeros(6,1);
   
    % Velocities
    dxdt(1) = x(4);
    dxdt(2) = x(5);
    dxdt(3) = x(6);
    
    % Positions
    dxdt(4) = - k/m(1) .* (x(1) - xin ) - k/m(1).*(x(1) - x(2)) - b/m(1).*x(4);
    dxdt(5) = - k/m(2) .* (x(2) - x(1)) - k/m(2).*(x(2) - x(3)) - b/m(2).*x(5);
    dxdt(6) = - k/m(3) .* (x(3) - x(2)) - k/m(3).*(x(3)       ) - b/m(3).*x(6);
end

[Sam Chak]

@Dirk te Brake, Please note that the External Input Force is not provided.

[Dirk te Brake]

Sorry maybe my question was a bit rushed the minus signs are indeed correct. I've gotten the same plot as you showed before. I want to compare my model with real world data. But the model doesn't match the data. I forgot to change the amplitude of the input to 0.05 m, but this doesnt impact the model (only scales it linearly). The derivation steps are shown below:

Code to read in the data:

direction = 'Path';
fileStructure = dir(fullfile(direction,'*.xlsx'));

name = fileStructure.name;
file = fullfile(direction, name);
data = readtable(file,'FileType','spreadsheet');

%Select either 'Links', 'Midden', 'Rechts'
xResults = fileStructure.data.Rechts;

plot(1:length(xResults),xResults)

[Sam Chak]

Hi @Dirk te Brake

Could you please extract the data and save it in the mat-file? Additionally, ensure that the code is updated where amplitude = 0.05 should be. If the simulated data doesn't match the real-world data, it's highly likely that your proposed mathematical model doesn't accurately capture the real-world dynamics. While the block masses can be measured fairly accurately with the right instruments, the dashpot coefficient (b) and the spring stiffness (k) need to be determined experimentally or through system identification tools.

@Star Strider has extensive experience in this topic and has provided numerous solutions. It is advised that you revise or update the question with additional context based on their expertise.

The solutions provided in the following two links may be helpful for your questions:

• Finding Unknown Parameters from Monod Equations: link
• Parameters Identification Providing Derivative: link

Since you have real-world data, the codes provided in these solutions could be beneficial.

[Dirk te Brake]

Hello @Sam Chak

Thanks for your answer, i will look into the links you have provided. I updated the orginal code given and added the .mat file. As you can see the data doesnt overlap that well, unless you drasticly change the spring stiffness. But the spring stiffness was given by my teacher so i don't think he would be of by a factor of 10^3 or something.

[Sam Chak]

Hi @Dirk te Brake

Thank you for the update. I have a question about the 7th state equation (dxdt(7) = xin). It doesn't appear to represent actuator dynamics. Could you please provide some clarification?

[Dirk te Brake]

As you mentioned before i dont provided a F_in but instead i give an postion as input.

So if i would like to compare the output with the input i can easly get the input. It is just an easy way of getting the input out of the function, atleast that was my thought process.

Thank you for the quick replies!

[Sam Chak]

I'm not sure I understand this part. An additional ODE, dx7 = xin, is included in the three sets of second-order motion systems. The coupled motion system consists of three displacements (x1 - x3) and three velocities (x4 - x6). The extra state, x7, represents the integral of xin. Could you please explain its significance?

[Dirk te Brake]

I forgot about it representing the integral, forget about that part.

[Sam Chak]

Hi @Dirk te Brake

Based on the laws of physics, the derived model of the mass-spring-damper system with the external sinusoidal force, cannot accurately reproduce the response that matches the provided real-world measured unscaled data.

Firstly, there appears to be a delay in the real-world data, but this delay is not accounted for in the mathematical model. Secondly, when coupled mass systems are subjected to an external undamped sinusoidal input signal, they typically exhibit repetitive oscillatory patterns at specific periodic intervals. However, no such repetitive patterns are observed in the real-world data.

Hence, there might be some missing information. I recommend checking with your teacher for further clarification.

load("RealWorldData.mat")
plot(RealX), grid on, xlabel('Number of Samples'), ylabel('Amplitude'), title('Real-world data')

[Dirk te Brake]

Hello @Sam Chak,

Thanks for the help! I will take a closer look at the mathematical model.