How to use lagrange equations for pendulum

7 次查看(过去 30 天)
Below is the code for symbolically simulating a pendulum, the plot produce doesn't seem to be the response of a pendulum swinging back and forth.
% Use Lagrange equations
% d / dt (d L / d (d qi / dt)) - d L / d qi = Q
%
% Simple pendulum
%%Without constraings
% Use angle as dof
% v = L d theta / dt
% K = 1/2 m v^2
% V = -mg L cos(theta)
% L = K - V = 1/2 m v^2 + m g L cos(theta)
syms Len d theta(t) m g
arc = Len * theta;
v = diff(arc,t);
K = 1/2 * m * v^2;
V = m*g*Len*(1-cos(theta));
L = K - V;
syms dtheta_dt
L1 = subs(L,diff(theta(t), t), dtheta_dt);
L2 = subs(diff(L1,dtheta_dt), dtheta_dt, diff(theta,t));
L3 = diff(L2,t);
syms thta
L4 = subs(L, theta, thta);
L5 = diff(L4, thta);
L6 = subs(L5, thta, theta);
eqn_pend = L3 + L6 == 0
[eqs_pend,vars_pend] = reduceDifferentialOrder(eqn_pend,theta(t))
[Mpend,Fpend] = massMatrixForm(eqs_pend,vars_pend)
syms Dtheta_Vart(t) dthta_dt;
MM = matlabFunction(Mpend, 'vars', {t, [thta; dthta_dt], Len,g,m})
Fpend1 = subs(Fpend, theta(t), thta);
Fpend2 = subs(Fpend1, Dtheta_Vart(t), dthta_dt);
FF = matlabFunction(Fpend2,'vars',{t,[dthta_dt;thta],Len,g,m})
MM_Fixed = @(t, in2)MM(t, in2, 5, 9.8, 100);
FF_Fixed = @(t, in2)FF(t, in2, 5, 9.8, 100);
opt = odeset('Mass', MM_Fixed);
[ts, ys]=ode15s(FF_Fixed, [0,10], [.0001; 0], opt);
figure;
plot(ts, ys);
legend('angle','rate');
  2 个评论
John D'Errico
John D'Errico 2017-12-8
编辑:John D'Errico 2017-12-8
No. This is arbitrary code that you obtained from some source. That it truly simulates a pendulum is not proven at all. Contact the source that provided the code, and ask them.
John
John 2017-12-8
John, I wrote the code trying to follow 3 refs. I recognize it's a simple problem, but I need to start somewhere.

请先登录,再进行评论。

回答(1 个)

John
John 2017-12-8
Thanks for the clues by John D Errico. There were two errors, the eqn_pend should have a minus sign, and the FF line reversed the order of the inputs. Here is the corrected code.
% Use Lagrange equations
% d / dt (d L / d (d qi / dt)) - d L / d qi = Q
%
% Simple pendulum
%%Without constraings
% Use angle as dof
% v = L d theta / dt
% K = 1/2 m v^2
% V = -mg L cos(theta)
% L = K - V = 1/2 m v^2 + m g L cos(theta)
syms Len d theta(t) m g
arc = Len * theta;
v = diff(arc,t);
K = 1/2 * m * v^2;
V = m*g*Len*(1-cos(theta));
L = K - V;
syms dtheta_dt
L1 = subs(L,diff(theta(t), t), dtheta_dt);
L2 = subs(diff(L1,dtheta_dt), dtheta_dt, diff(theta,t));
dL_d_theta_dt = diff(L2,t);
syms thta
L4 = subs(L, theta, thta);
L5 = diff(L4, thta);
dL_d_theta = subs(L5, thta, theta);
eqn_pend = dL_d_theta_dt - dL_d_theta == 0
[eqs_pend,vars_pend] = reduceDifferentialOrder(eqn_pend,theta(t))
[Mpend,Fpend] = massMatrixForm(eqs_pend,vars_pend)
syms dthta_dt;
MM = matlabFunction(Mpend, 'vars', {t, [thta; dthta_dt], Len,g,m})
Fpend1 = subs(Fpend, theta(t), thta);
Fpend2 = subs(Fpend1, vars_pend(2), dthta_dt);
FF = matlabFunction(Fpend2,'vars',{t,[thta;dthta_dt],Len,g,m})
MM_Fixed = @(t, in2)MM(t, in2, 5, 9.8, 100);
FF_Fixed = @(t, in2)FF(t, in2, 5, 9.8, 100);
opt = odeset('Mass', MM_Fixed);
[ts, ys]=ode15s(FF_Fixed, [0,100], [0; .1], opt);
figure;
plot(ts, ys);
legend('angle','rate');
end

类别

Help CenterFile Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by