Solution of 2nd order differential equation

I am trying to validate the plots of a research paper. But whenever I try to solve I don't get the appropriate plot. There is a difference in the plot. I am providing the equation and the values of different constants. I am also providing the plot which is given in the research paper. Kindly guide me for this.

The Highlighted one is the governing equation.

Plot given in the research paper given below.

I m providing you the code I have tried and the result obtained below.

% Constants

m = 0.89694;

m0 = 0.012992;

c = 0.0689097;

g = 9.81;

F0 = 124.728;

alpha = 54.5841;

e = 0.050;

f = 5;

omega = 2 * pi * f;

t = linspace(0, 100, 100000);

Fm = @(x) F0 ./ (1 + alpha * x).^4;

y0 = [0; 0];

[t, sol] = ode45(@(t, y) model(t, y, m, m0, c, g, Fm, omega, e), t, y0);

x = sol(:, 1);

x_dot = sol(:, 2);

% Plot time-displacement

figure(1);

plot(t, x,'r');

xlabel('Time (s)');

ylabel('Displacement (m)');

title('Time vs Displacement');

grid on;

% Define the differential equation

function dydt = model(t, y, m, m0, c, g, Fm, omega, e)

x = y(1); % Displacement

x_dot = y(2); % Velocity

total_mass = m + m0;

% Acceleration term

x_ddot = (Fm(x) + m0 * e * (omega^2) * sin(omega * t) - c * x_dot - total_mass * g) / total_mass;

dydt = [x_dot; x_ddot]; % Return [velocity; acceleration]

end

[Shashi Kiran]

Hi @Parthajit,

Upon reviewing the code you provided, I recommend two improvements:

1. Assign a small initial displacement to match the oscillation amplitude, such as setting y0 = [0.01; 0]
2. Adjust the frequency f by reducing it or modifying it as needed to better fit the data, for example, using f = 4.5.

[Torsten]

You don't want to include the "big middle term" of equation (32) ?

[Parthajit]

The plot is for equation no 31.

[Sam Chak]

Have you checked if the initial values are correct?

tspan   = [0 100];
x0      = [0.025; 0.016];   % <-- Educated guess only
option  = odeset('RelTol', 1e-2, 'AbsTol', 1e-5);
[t, x]  = ode45(@model, tspan, x0, option);

% Plot time-displacement
figure(1);
plot(t, x(:,1));
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Time vs Displacement');
grid on;

% Define the differential equation
function dydt = model(t, y)
    % Constants
    m       = 0.89694;
    m0      = 0.012992;
    c       = 0.0689097;
    g       = 9.81;
    F0      = 124.728;
    alpha   = 54.5841;
    e       = 0.050;
    f       = 5;
    omega   = 2*pi*f;

    x       = y(1);         % Displacement
    x_dot   = y(2);         % Velocity
    total_mass = m + m0;
    
    Fm      = F0/(1 + alpha*x)^4;
    
    % Acceleration term
    x_ddot  = (Fm + m0*e*(omega^2)*sin(omega*t) - c*x_dot - total_mass*g)/total_mass;
    
    dydt    = [x_dot; x_ddot];  % Return [velocity; acceleration]
end

[Parthajit]

Thank you for your suggestion sir @Shashi Kiran. I have tried with with what you have suggested. I get quite similar result as in provided in the research paper.

But with initial condition 0.01 unit and frequency 5 I am not getting similar profile. But in the paper they clearly mentioned that frequency taken as 5. Please shareif you have any comment regarding this.

Thank you once again, Sir.

[Parthajit]

@Sam Chak All the constant values are given in the paper and I have checked it again and again. These are correct and talking about the initial conditions they did not mention , only thing is that its a forced vibration system, where external harmonic excitation are provided. Initial conditions may be provided for displacement as 0.01744m some positive value and velocity 0.

[VBBV]

@Parthajit, You can keep the frequency value at 5 as given in paper, but check the imbalance term value, e in the paper

% Constants
m = 0.89694;
m0 = 0.012992;
c = 0.0689097;
g = 9.81;
F0 = 124.728;
alpha = 54.5841;
e = 0.004;
f = 5;
omega = 2 * pi * f;
t = linspace(0, 100, 100000);
Fm = @(x) F0 ./ (1 + alpha * x).^4;
y0 = [0.01; 0];
[t, sol] = ode45(@(t, y) model(t, y, m, m0, c, g, Fm, omega, e), t, y0);
x = sol(:, 1);
x_dot = sol(:, 2);
% Plot time-displacement
figure(1);subplot(411)
plot(t, x,'r');xlim([35 40])
subplot(412)
plot(t, x,'r');xlim([15 30])
subplot(413)
plot(t(1:200:end), x(1:200:end),'r');
 
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Time vs Displacement');
grid on;

% Define the differential equation
function dydt = model(t, y, m, m0, c, g, Fm, omega, e)
    x = y(1);         % Displacement
    x_dot = y(2);     % Velocity
    total_mass = m + m0;
    % Acceleration term
    x_ddot = (Fm(x) + m0 * e * (omega^2) * sin(omega * t) - c * x_dot - total_mass * g) / total_mass;
    dydt = [x_dot; x_ddot];  % Return [velocity; acceleration]
end

[Parthajit]

@VBBV No sir, in the paper they provided "e" as 5 mm.

[VBBV]

@Parthajit Ok, that means it is 0.005 ~ 5mm, and not 0.05 ~ 5cm which you had used earlier.

[Parthajit]

@VBBV Sorry Sir, I mistakenly wrote it 5 mm, e is given as 50 mm.

[VBBV]

@Parthajit the values given show that the graphs are computed using Eq. 32 and not Eq. 31

[Parthajit]

@VBBV Okay sir I am trying with that, actually no Xst value is provided in the paper, but this Xst value is equal to the initial value of x , which is also not provided directly (can be seen in the plot only).