How to use lagrange equations for pendulum

11 views (last 30 days)
John
John on 8 Dec 2017
Commented: John on 8 Dec 2017
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 Comments
John D'Errico
John D'Errico on 8 Dec 2017
Edited: John D'Errico on 8 Dec 2017
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 on 8 Dec 2017
John, I wrote the code trying to follow 3 refs. I recognize it's a simple problem, but I need to start somewhere.

Sign in to comment.

Answers (1)

John
John on 8 Dec 2017
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

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!