ODE 45 gives unexpected result on pendulum tuned mass damper

I'm trying to obtain a numerical solution for the motion of a pendulum mass damper using nonlinear differential equations obtained from Lagrange's equation.

I am changing the mass ratio to observe the resulting performance. It is expected that as the ratio gets closer to 1, the displacement of the main structure should be smaller.

However, the ode45 solver is providing completely different results. What could be wrong with my code?

I have used the state variable from an essay, so it should be correct. I have replaced "u" with "M*ag".

I also replaced the symbols like this: "x1" remains "x1", "x2" remains as "x2", "theta1" with "x3", and "theta2" with "x4".

clear;

ratio = [1,0.5,0.25];

for r = 1:numel(ratio)

    M = 0.400;
    m = M*ratio(r);
    
    [t,y] = ode45(@(t,y) f(t,y,m), [0 20], [0,0,0,0]);


    %%Finding Acceleration 
    accelerationDifference = diff(y(:, 2));
    timeDifference = diff(t);
    a = accelerationDifference ./ timeDifference;
    a(end+1) = a(end-1);

    subplot(4,1,1); 
    hold on;
    plot(t, y(:,1), 'DisplayName', [num2str(ratio(r)), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('Displacement vs Time graph of primary structure');
    legend('show');
    
    
    subplot(4,1,2); 
    hold on;
    plot(t, y(:,2), 'DisplayName', [num2str(ratio(r)), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Velocity (m/s)');
    title('Velocity vs Time graph of primary structure');
    legend('show');
    
    subplot(4, 1, 3);
    hold on;
    plot(t, a, 'DisplayName', [num2str(ratio(r)), ' Ratio']); 
    xlabel('Time (s)');
    ylabel('Acceleration (m/s^2)');
    title('Acceleration vs Time graph of primary structure');
    legend('show'); 
        
    subplot(4,1,4);
    hold on;
    plot(t, y(:,3), 'DisplayName', [num2str(m), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('PTMD Displacement vs Time ');
    legend('show'); 

end

function dydt = f(t,y,m)

    L = 0.058;
    M = 0.400;
    b = 0.1353; 
    d = 0.15;
    c = 100*0.000098*60/(20*2*pi); 
    g = 9.807;
    k = 67; 
    ag = 0.025*(2*pi/0.9)^2*sin(2*pi/0.9*t);

    x1 = y(1); % x1 = displacement of main structure
    x2 = y(2); % x2 = velocity of main structure
    theta1 = y(3); %theta1 = angular displacement of pendulum
    theta2 = y(4); %theta1 = angular velocity of pendulum
    
    x3 = ( M*ag + m*g*cos(theta1)*sin(theta1) + cos(theta1)*(c/L+d*L*cos(theta1)^2)*theta2 - k*x1 - b*x2 + m*L*theta2^2*sin(theta1) ) / ( M + m*(1 - cos(theta1)^2) );
    theta3 = ( -(M+m)/m*(c/L+d*L*cos(theta1)^2)*theta2 - (M+m)*g*sin(theta1) - cos(theta1)*(m*L*theta2^2*sin(theta1) + M*ag - k*x1 + b*x2) ) / ( L*(M + m*(1 - cos(theta1)^2)) );
   
   
    dydt = [y(2);x3;y(4);theta3];
end

[Torsten]

I interpreted F(t) = M*ag.

ratio = [1,0.5,0.25];

for r = 1:numel(ratio)

    M = 0.400;
    m = M*ratio(r);
    
    [t,y] = ode45(@(t,y) f(t,y,m), [0 20], [0,0,0,0]);

    a = zeros(size(t));
    for i = 1:numel(t)
      dydt = f(t(i),y(i,:),m);
      a(i) = dydt(2); 
    end
    
    %%Finding Acceleration 
    %accelerationDifference = diff(y(:, 2));
    %timeDifference = diff(t);
    %a = accelerationDifference ./ timeDifference;
    %a(end+1) = a(end-1);

    subplot(4,1,1); 
    hold on;
    plot(t, y(:,1), 'DisplayName', [num2str(ratio(r)), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('Displacement vs Time graph of primary structure');
    legend('show');
    
    
    subplot(4,1,2); 
    hold on;
    plot(t, y(:,2), 'DisplayName', [num2str(ratio(r)), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Velocity (m/s)');
    title('Velocity vs Time graph of primary structure');
    legend('show');
    
    subplot(4, 1, 3);
    hold on;
    plot(t, a, 'DisplayName', [num2str(ratio(r)), ' Ratio']); 
    xlabel('Time (s)');
    ylabel('Acceleration (m/s^2)');
    title('Acceleration vs Time graph of primary structure');
    legend('show'); 
        
    subplot(4,1,4);
    hold on;
    plot(t, y(:,3), 'DisplayName', [num2str(m), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('PTMD Displacement vs Time ');
    legend('show'); 

end

function dydt = f(t,y,m)

    L = 0.058;
    M = 0.400;
    b = 0.1353; 
    d = 0.15;
    c = 100*0.000098*60/(20*2*pi); 
    g = 9.807;
    k = 67; 
    ag = 0.025*(2*pi/0.9)^2*sin(2*pi/0.9*t);

    x1 = y(1); % x1 = displacement of main structure
    x2 = y(2); % x2 = velocity of main structure
    theta1 = y(3); %theta1 = angular displacement of pendulum
    theta2 = y(4); %theta1 = angular velocity of pendulum
    
    amat(1,1) = M+m;
    amat(1,2) = m*L*cos(theta1);
    amat(2,1) = cos(theta1);
    amat(2,2) = L;
    bvec(1,1) = M*ag -k*x1 - b*x2 + m*L*theta2^2*sin(theta1);
    bvec(2,1) = -g*sin(theta1) - 1/m*(c/L+d*L*cos(theta1)^2)*theta2;
    
    result = amat\bvec;
    
    dydt = [y(2);result(1);y(4);result(2)];
    
    %x3 = ( M*ag + m*g*cos(theta1)*sin(theta1) + cos(theta1)*(c/L+d*L*cos(theta1)^2)*theta2 - k*x1 - b*x2 + m*L*theta2^2*sin(theta1) ) / ( M + m*(1 - cos(theta1)^2) );
    %theta3 = ( -(M+m)/m*(c/L+d*L*cos(theta1)^2)*theta2 - (M+m)*g*sin(theta1) - cos(theta1)*(m*L*theta2^2*sin(theta1) + M*ag - k*x1 + b*x2) ) / ( L*(M + m*(1 - cos(theta1)^2)) );
   
   
    %dydt = [y(2);x3;y(4);theta3];
end

[Sam Chak]

@Baker, The Jacobian matrix shows that the pendulum is stable under small perturbations around the equilibrium point , as all the eigenvalues of the matrix have negative real parts. Therefore, the displacement trajectories converge to 0 even when the force .

Why do you think that the displacement should get smaller when ?

ratio = [1, 0.5, 0.25];

for r = 1:numel(ratio)

    M = 0.400;
    m = M*ratio(r);
    
    [t,y] = ode45(@(t,y) f(t,y,m), [0 20], [0.1, 0, 0, 0]);

    a = zeros(size(t));
    for i = 1:numel(t)
      dydt = f(t(i),y(i,:),m);
      a(i) = dydt(2); 
    end
    
    %%Finding Acceleration 
    %accelerationDifference = diff(y(:, 2));
    %timeDifference = diff(t);
    %a = accelerationDifference ./ timeDifference;
    %a(end+1) = a(end-1);

    % subplot(4,1,1); 
    hold on;
    plot(t, y(:,1), 'DisplayName', [num2str(ratio(r)), 'Ratio']); 
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('Displacement vs Time graph of primary structure');
    legend('show');
    
    
    % subplot(4,1,2); 
    % hold on;
    % plot(t, y(:,2), 'DisplayName', [num2str(ratio(r)), 'Ratio']); 
    % xlabel('Time(s)');
    % ylabel('Velocity (m/s)');
    % title('Velocity vs Time graph of primary structure');
    % legend('show');
    % 
    % subplot(4, 1, 3);
    % hold on;
    % plot(t, a, 'DisplayName', [num2str(ratio(r)), ' Ratio']); 
    % xlabel('Time (s)');
    % ylabel('Acceleration (m/s^2)');
    % title('Acceleration vs Time graph of primary structure');
    % legend('show'); 
    % 
    % subplot(4,1,4);
    % hold on;
    % plot(t, y(:,3), 'DisplayName', [num2str(m), 'Ratio']); 
    % xlabel('Time(s)');
    % ylabel('Displacement(m)');
    % title('PTMD Displacement vs Time ');
    % legend('show'); 

end

function dydt = f(t,y,m)

    L  = 0.058;
    M  = 0.400;
    b  = 0.1353; 
    d  = 0.15;
    c  = 100*0.000098*60/(20*2*pi); 
    g  = 9.807;
    k  = 67; 
    ag = 0.025*(2*pi/0.9)^2*sin(2*pi/0.9*t);

    x1 = y(1); % x1 = displacement of main structure
    x2 = y(2); % x2 = velocity of main structure
    theta1 = y(3); %theta1 = angular displacement of pendulum
    theta2 = y(4); %theta1 = angular velocity of pendulum
    
    amat(1,1) = M+m;
    amat(1,2) = m*L*cos(theta1);
    amat(2,1) = cos(theta1);
    amat(2,2) = L;
    F         = 0;
    bvec(1,1) = F - k*x1 - b*x2 + m*L*theta2^2*sin(theta1);
    bvec(2,1) = - g*sin(theta1) - 1/m*(c/L + d*L*cos(theta1)^2)*theta2;
    
    result = amat\bvec;
    
    dydt = [y(2);
            result(1);
            y(4);
            result(2)];
    
    %x3 = ( M*ag + m*g*cos(theta1)*sin(theta1) + cos(theta1)*(c/L+d*L*cos(theta1)^2)*theta2 - k*x1 - b*x2 + m*L*theta2^2*sin(theta1) ) / ( M + m*(1 - cos(theta1)^2) );
    %theta3 = ( -(M+m)/m*(c/L+d*L*cos(theta1)^2)*theta2 - (M+m)*g*sin(theta1) - cos(theta1)*(m*L*theta2^2*sin(theta1) + M*ag - k*x1 + b*x2) ) / ( L*(M + m*(1 - cos(theta1)^2)) );
   
   
    %dydt = [y(2);x3;y(4);theta3];
end

[Baker]

@Sam Chak

I think will decrease when because a pendulum tuned mass damper should reduce the amplitude of the primary structure. The author of the essay used linearization to process the differential equations, and this approach resulted in a solution that decreases as m increases. Do you think the process of linearization affects the outcome?

[Torsten]

Why don't you solve the linearized equations and see what happens ?

[Baker]

@Torsten I tried and it gives really big values that data points appear to be a stright line.

clear;close;

ratio = [1,0.5,0.25];
averageMaxAcceleration = zeros(1, numel(ratio));

for r = 1:numel(ratio)
    asum=0;
    avga=0;
    M = 0.400;
    m = M*ratio(r);
    
    % Solve ODE for current value of m
    [t,y] = ode45(@(t,y) f(t,y,m), [0 10], [0,0,0,0]);

    accelerationDifference = diff(y(:, 2));
    timeDifference = diff(t);
    a = accelerationDifference ./ timeDifference;
    a(end+1) = a(end-1);


    subplot(4,1,1);  
    
    hold on;
    plot(t, y(:,1), 'DisplayName', [num2str(m), 'Ratio']); % Add label with ratio
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('Displacement vs Time graph of primary structure');
    legend('show'); % Show legend with set labels
    
    
    subplot(4,1,2); % Create second subplot for velocity
    hold on;
    plot(t, y(:,2), 'DisplayName', [num2str(m), 'Ratio']); % Add label with ratio
    xlabel('Time(s)');
    ylabel('Velocity (m/s)');
    title('Velocity vs Time graph of primary structure');
    legend('show'); % Show legend with set labels
    
    subplot(4, 1, 3); % Create third subplot for acceleration
    hold on;
    plot(t, a, 'DisplayName', [num2str(ratio(r)), ' Ratio']); % Add label with ratio
    xlabel('Time (s)');
    ylabel('Acceleration (m/s^2)');
    title('Acceleration vs Time graph of primary structure');
    legend('show'); % Show legend with set labels

    subplot(4,1,4);
    hold on;
    plot(t, y(:,3), 'DisplayName', [num2str(m), 'Ratio']); % Add label with ratio
    xlabel('Time(s)');
    ylabel('Displacement(m)');
    title('PTMD Displacement vs Time ');
    legend('show'); % Show legend with set labels
        

end

function dydt = f(t,y,m)

    L = 0.058; 
    M = 0.400;
    b = 0.1353; 
    d = 0.15;
    c = 100*0.000098*60/(20*2*pi); 
    g = 9.807;
    k = 67; 
    ag = 0.025*(2*pi/0.9)^2*sin(2*pi/0.9*t);


    F = M*ag;
    m_r = m/M;
    w = sqrt(k/M);
    r = b/ (2*w*M);
    w_d = sqrt(g/L);
    r_d = (c/L+d*L) / (2*w_d*m*L);



    % Unpack state variables
    x1 = y(1); % x1 = displacement of main structure
    x2 = y(2); % x2 = velocity of main structure
    y1 = y(3); % y1 = linearised displacement of main structure
    y2 = y(4); % y1 = linearised velocity of main structure
    
    %Essay Direct Copy - I made x3 and y3 as subject 
    x3 = 2*m_r*r_d*w_d*y2 + m_r*w_d^2*y1 -1/M*F -2*r*w*x2 -w^2*x1;
    y3 = 2*r*w*x2 + w^2*x1 +1/M*F -2*(1+m_r)*r_d*w_d*y2 - (1+m_r)*w_d^2*x1;
   
    dydt = [y(2);x3;y(4);y3];

 
    
    
   

end

[Baker]

The author obtains a sum of homogenous and non-homogenous solutions tho, maybe I should work on that.

[Sam Chak]

Hi @Baker,

You are correct in stating that a larger mass ratio of Tuned Mass Damper (TMD) offers significant advantages in controlling the displacement response of the main structure, as mentioned by Kang & Peng in 2019. You can also refer to Figure 15 in the following link: https://www.hindawi.com/journals/sv/2019/8376781/fig15/.

However, I recommend verifying the equations. It is crucial to accurately derive the Euler–Lagrange equations. Could you please provide the expressions for kinetic energy (T) and potential energy (V) of the system based on the Lagrangian ?

Furthermore, you can utilize MATLAB to confirm the accuracy of the derivation of Eq. (3-15).

[Torsten]

Can you please explain what are you refering in expression (3-15) and (3-18) with the arrow?

Should be -b/(L*M) instead of b/(L*M) in the matrix.

[Baker]

@Torsten Oic. Thank you for explaining.

@Sam Chak oops, I accidently removed my comment.

The following picture is a part of the essay. And I derived the equations again and i got the same expressions.

Answer 1

[Sam Chak]

Hi @Baker

You can use the odeToVectorField() function to derive the dynamical model from the Euler–Lagrange equations. The accuracy of the derivation depends on the defined Lagrangian.

L  = 0.058;
M  = 0.400;
mu = 1;
m  = M*mu;
g  = 9.807;
k  = 67; 
b  = 0.1353; 
c  = 100*0.000098*60/(20*2*pi); 
d  = 0.15;
sympref('AbbreviateOutput', false);

syms x(t) theta(t)

% definitions
Dx     = diff(x, 1);
Dtheta = diff(theta, 1);

% Langrangian
T      = 1/2*M*Dx^2 + 1/2*m*((Dx + L*Dtheta*cos(theta))^2 + (L*Dtheta*sin(theta))^2);
V      = 1/2*k*x^2 + m*g*L*(1 - cos(theta));
L      = T - V;
ag     = 0.025*(2*pi/0.9)^2*sin(2*pi/0.9*t);
F      = M*ag;

% Euler–Lagrange equations
eqn1   = diff(diff(L, Dx), 1)     - diff(L, x)     == F - b*Dx;
eqn2   = diff(diff(L, Dtheta), 1) - diff(L, theta) == - c*Dtheta - d*Dtheta*(L*cos(theta))^2;
eqns   = [eqn1 eqn2];
[V, S] = odeToVectorField(eqns)

V = 

S = 

dydt   = matlabFunction(V, 'vars', {'t', 'Y'});
tspan  = [0 20];                % time span of simulation
y0     = [0 0 0.1 0];         % initial values
[t, Y] = ode45(@(t, Y) dydt(t, Y), tspan, y0);

plot(t, Y(:,3)), grid on
xlabel('Time(s)');
ylabel('Displacement(m)');
title('Displacement vs Time graph of primary structure');