I have been working on a model for a pendulum(using a Lagrangian) and I am using NDSolve. The solver gives me outputs for alpha an l yet I can´t us these to recalculate x and y.
In[2]:= ClearAll["Global`*"]
Definition of parameters
In[3]:= R = 0
Subscript[l, 0] = 0.5
m = 1
M = 2.394
g = 9.81
Subscript[\[Alpha], 0] = 1/2*Pi
Expression of position(x,y) as a function of [Alpha]
In[9]:= x[\[Alpha]_, l_, t_] := -Cos[\[Alpha][t]]*R +
Sin[\[Alpha][t]]*(l[t] - \[Alpha][t]*R)
y[\[Alpha]_, l_, t_] :=
Sin[\[Alpha][t]]*R + Cos[\[Alpha][t]]*(l[t] - \[Alpha][t]*R)
Calculation of the components and the norm of the velocity
Computation of kinetic energy, potential energy and Lagrangian
In[11]:= dx[\[Alpha]_, l_, t_] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((x[\[Alpha], l, t])\)\)
dy[\[Alpha]_, l_, t_] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((y[\[Alpha], l, t])\)\)
v[\[Alpha]_, l_, t_] :=
Sqrt[dx[\[Alpha], l, t]^2 + dy[\[Alpha], l, t]^2] // FullSimplify;
Ekin[\[Alpha]_, l_, t_] := 1/2*m*v[\[Alpha], l, t]^2 + 1/2*M*(l'[t])^2
Epot[\[Alpha]_, l_, t_] := m*g*y[\[Alpha], l, t] + M*g (l[t] - Subscript[l, 0])
L1[\[Alpha]_, l_, t_] := Ekin[\[Alpha], l, t] - Epot[\[Alpha], l, t]
Euler-Lagrange equation
In[17]:= eq1[\[Alpha]_, l_, t_] =
D[D[L1[\[Alpha], l, t], \[Alpha]'[t]], t] -
D[L1[\[Alpha], l, t], \[Alpha][t]] // FullSimplify;
eq2[\[Alpha]_, l_, t_] =
D[D[L1[\[Alpha], l, t], l'[t]], t] - D[L1[\[Alpha], l, t], l[t]] //
FullSimplify;
eq = {eq1[\[Alpha], l, t] == 0, eq2[\[Alpha], l, t] == 0}
Out[19]= {l[t] (-9.81 Sin[\[Alpha][t]] +
2 Derivative[1][l][t] Derivative[1][\[Alpha]][t] +
l[t] (\[Alpha]^\[Prime]\[Prime])[t]) == 0,
23.4851 + 9.81 Cos[\[Alpha][t]] - l[t] Derivative[1][\[Alpha]][t]^2 +
3.394 (l^\[Prime]\[Prime])[t] == 0}
In[20]:=
Solution of equation of motion [AliasDelimiter]with given initial conditions
In[22]:= sol = First@
NDSolve[{eq, \[Alpha][0] == Subscript[\[Alpha], 0], \[Alpha]'[0] == 0,
l[0] == Subscript[l, 0], l'[0] == 0}, {\[Alpha][t], l[t]}, {t, 0, 20}]
Plot
In[23]:= Plot[\[Alpha][t] /. sol, {t, 0, 20}]
In[24]:= Plot[l[t] /. sol, {t, 0, 20}]
In[29]:=
ParametricPlot[{x[\[Alpha], t]/.sol, y[\[Alpha], t]/.sol}, {t, 0, 20}, PlotRange -> All]
