How to use lagrange equations for pendulum

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

[John D'Errico]

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, I wrote the code trying to follow 3 refs. I recognize it's a simple problem, but I need to start somewhere.

Answer 1

[John]

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