Trouble animating a double pendulum

9 次查看（过去 30 天）
Ronak
Ronak 2025-4-25
编辑： Torsten 2025-4-25
Hello! For a project in my physics class, I am animating a double pendulum, but I keep having trouble with the animation. It keeps accelerating and it doesn't seem accurate. I showed the code to my physics teacher, and he said that it all seemed fine and couldn't understand why it isn't working. Should I stick with the method I am using right now or switch to using ode45?
syms a_1 a_2
theta_1 = input('What do you want Theta 1 to equal? ')*pi/180;
theta_2 = input('What do you want Theta 2 to equal? ')*pi/180;
l_1 = 1;
l_2 = 1;
g = 9.8;
m_1 = 1;
m_2 = 1;
dt = .01;
w_1 = 0;
w_2 = 0;
for n=1:1000
clf
a_1=-(g*(2*m_1+m_2)*sin(theta_1)+m_2*g*sin(theta_1-2*theta_2)+2*sin(theta_1-theta_2)*m_2*(((w_2)^2)*l_2+((w_1)^2)*l_1*cos(theta_1-theta_2)))/(l_1*2*m_1+m_2-m_2*cos(2*(theta_1-theta_2)));
a_2 = (2*sin(theta_1-theta_2)*((w_1)^2)*l_1*(m_1+m_2)+g*(m_1+m_2)*cos(theta_1)+((w_2)^2)*l_2*m_2*cos(theta_1-theta_2))/(l_2*(2*m_1+m_2-m_2*cos(2*theta_1-2*theta_2)));
w_1 = w_1+a_1*dt;
w_2 = w_2+a_2*dt;
theta_1 = theta_1+w_1*dt;
theta_2 = theta_2+w_2*dt;
x_1 = l_1*sin(theta_1);
x_2 = l_1*sin(theta_1)+l_2*sin(theta_2);
y_1 = -l_1*cos(theta_1);
y_2 = -l_1*cos(theta_1)-l_2*cos(theta_2);
plot(0, 0, 'O')
hold on
plot(x_1, y_1, '*')
hold on
line([0 x_1], [0 y_1])
hold on
plot(x_2, y_2, 'o')
hold on
line([x_1 x_2], [y_1 y_2])
axis([-10 10 -10 10])
pause(.01)
end
  1 个评论
Torsten
Torsten 2025-4-25
编辑：Torsten 2025-4-25
Should I stick with the method I am using right now or switch to using ode45?
If you solve a problem with a self-written numerical code, you should always compare with a professional solver if possible.
Concerning your code:
Don't define a1 and a2 as symbolic.
And don't plot within the loop. Save x_1, x_2, y_1 and y_2 in an array of length 1001 and plot after the loop.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2025-4-25
编辑：Torsten 2025-4-25
Ran in:
theta_1 = 19*pi/180;
theta_2 = 30*pi/180;
l_1 = 1;
l_2 = 1;
g = 9.8;
m_1 = 1;
m_2 = 1;
M = m_1 + m_2;
dt = .001;
w_1 = 0;
w_2 = 0;
x_1 = zeros(10000,1);
y_1 = zeros(10000,1);
x_2 = zeros(10000,1);
y_2 = zeros(10000,1);
x_1(1) = l_1*sin(theta_1);
x_2(1) = x_1(1) + l_2*sin(theta_2);
y_1(1) = -l_1*cos(theta_1);
y_2(1) = y_1(1) - l_2*cos(theta_2);
for n=2:10000
% Source: https://physics.umd.edu/hep/drew/pendulum2.html
dtheta = theta_1 - theta_2;
alpha = m_1 + m_2*sin(dtheta)^2;
a_1 = (-sin(dtheta)*(m_2*l_1*w_1^2*cos(dtheta)+m_2*l_2*w_2^2)-...
g*(M*sin(theta_1)-m_2*sin(theta_2)*cos(dtheta)))/(l_1*alpha);
a_2 = (sin(dtheta)*(M*l_1*w_1^2+m_2*l_2*w_2^2*cos(dtheta))+...
g*M*(sin(theta_1)*cos(dtheta)-sin(theta_2)))/(l_2*alpha);
%dtheta = theta_1 - theta_2;
%Mat = [M*l_1,m_2*l_2*cos(dtheta);l_1*cos(dtheta),l_2];
%rhs = [-m_2*l_2*w_2^2*sin(dtheta)-M*g*sin(theta_1);l_1*w_1^2*sin(dtheta)-g*sin(theta_2)];
%sol = Mat\rhs;
%a_1 = sol(1);
%a_2 = sol(2);
w_1_new = w_1+a_1*dt;
w_2_new = w_2+a_2*dt;
theta_1_new = theta_1+w_1*dt;
theta_2_new = theta_2+w_2*dt;
w_1 = w_1_new;
w_2 = w_2_new;
theta_1 = theta_1_new;
theta_2 = theta_2_new;
x_1(n) = l_1*sin(theta_1);
x_2(n) = x_1(n) + l_2*sin(theta_2);
y_1(n) = -l_1*cos(theta_1);
y_2(n) = y_1(n) - l_2*cos(theta_2);
end
hold on
plot(x_1,y_1)
plot(x_2,y_2)
hold off
grid on
  1 个评论
Sam Chak
Sam Chak 2025-4-25
A double pendulum is a chaotic system, and the second pendulum will exhibit unpredictable patterns of motion. If the simulation responses appear predictable, either the equations are incorrect or the ratios of lengths and masses are carefully tuned, resulting in more predictable behavior.
Chaotic behavior:
Predictable behavior:

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Physics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by