ODE 45 gives unexpected result on pendulum tuned mass damper

Question

Baker 2023-8-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2009187-ode-45-gives-unexpected-result-on-pendulum-tuned-mass-damper

回答： Sam Chak 2023-8-19

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

Torsten 2023-8-16

在 MATLAB Online 中打开

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 2023-8-16

在 MATLAB Online 中打开

@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 2023-8-17

编辑：Baker 2023-8-17

在 MATLAB Online 中打开

@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

Torsten 2023-8-18

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 2023-8-18

编辑：Baker 2023-8-18

@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 2023-8-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2009187-ode-45-gives-unexpected-result-on-pendulum-tuned-mass-damper#answer_1289732

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

ODE 45 gives unexpected result on pendulum tuned mass damper

9 个评论
显示 7更早的评论隐藏 7更早的评论

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

ODE 45 gives unexpected result on pendulum tuned mass damper

9 个评论 显示 7更早的评论隐藏 7更早的评论

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

9 个评论
显示 7更早的评论隐藏 7更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论