Double pendulum simulations

This simulates a trajectory of the double pendulum with (Subscript[l, 1]) = (Sqrt[2]), (Subscript[l, 2]) = 1, (Subscript[m, 1]) = (Subscript[m, 2]) = 1 and initial conditions (Subscript[θ, 1])(0)= (π/4), (Subscript[θ, 2])(0) = π, (Subscript[θ, 1])'(0) = (Subscript[θ, 2])'(0) = 0.

deqns = {(Subscript[m, 1]) (Subscript[x, 1])''[t] == ((((Subscript[λ, 1])[t])/(Subscript[l, 1])))(Subscript[x, 1])[t] - ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[x, 2])[t] - (Subscript[x, 1])[t]),
(Subscript[m, 1]) (Subscript[y, 1])''[t] == ((((Subscript[λ, 1])[t])/(Subscript[l, 1])))(Subscript[y, 1])[t] - ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[y, 2])[t] - (Subscript[y, 1])[t])-(Subscript[m, 1]) g,
(Subscript[m, 2]) (Subscript[x, 2])''[t] ==  ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[x, 2])[t] - (Subscript[x, 1])[t]),
(Subscript[m, 2]) (Subscript[y, 2])''[t] ==  ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[y, 2])[t] - (Subscript[y, 1])[t]) - (Subscript[m, 2]) g};
aeqns = {(Subscript[x, 1])[t]^2 + (Subscript[y, 1])[t]^2 == (Subscript[l, 1])^2, ( (Subscript[x, 2])[t]-(Subscript[x, 1])[t])^2 +((Subscript[y, 2])[t]- (Subscript[y, 1])[t])^2 == (Subscript[l, 2])^2};
ics = {(Subscript[x, 1])[0] ==  1, (Subscript[y, 1])[0] == -1, (Subscript[x, 1])'[0] == 0, (Subscript[y, 1])'[0] == 0, (Subscript[x, 2])[0] == 1 , (Subscript[y, 2])[0] == 0, (Subscript[x, 2])'[0] == 0, (Subscript[y, 2])'[0] == 0};
params = {g -> 9.81, (Subscript[m, 1]) -> 1, (Subscript[m, 2])-> 1, (Subscript[l, 1]) -> Sqrt[2], (Subscript[l, 2]) -> 1};
soldp = First[NDSolve[{deqns, aeqns,ics} /. params,
 {(Subscript[x, 1]), (Subscript[y, 1]), (Subscript[x, 2]), (Subscript[y, 2]),(Subscript[λ, 1]), (Subscript[λ, 2])}, {t, 0, 15},
Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}]];

Animate[
Graphics[{{PointSize[.025], {Red, Point[{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}]},
{Blue, Point[{(Subscript[x, 2])[t], (Subscript[y, 2])[t]}]}, Line[{{0,0}, {(Subscript[x, 1])[t], (Subscript[y, 1])[t]},{(Subscript[x, 2])[t], (Subscript[y, 2])[t]}}]} /. soldp,
{Gray, Line[Map[Function[Evaluate[{(Subscript[x, 2])[#], (Subscript[y, 2])[#]} /. soldp]], Range[0,t,0.025]]]}},
PlotRange -> {{-4,4},{-4,0}}, Axes -> True, Ticks -> False, ImageSize -> 500],
{t,0,15,.025}, SaveDefinitions -> True]

This simulates a trajectory of the double pendulum with (Subscript[l, 1]) = (Sqrt[2]), (Subscript[l, 2]) = 1, (Subscript[m, 1]) = (Subscript[m, 2]) = 1 and initial conditions (Subscript[θ, 1])(0)= (π/4), (Subscript[θ, 2])(0) = π - .03, (Subscript[θ, 1])'(0) = (Subscript[θ, 2])'(0) = 0.

deqns2 = {(Subscript[m, 1]) (Subscript[x, 1])''[t] == ((((Subscript[λ, 1])[t])/(Subscript[l, 1])))(Subscript[x, 1])[t] - ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[x, 2])[t] - (Subscript[x, 1])[t]),
(Subscript[m, 1]) (Subscript[y, 1])''[t] == ((((Subscript[λ, 1])[t])/(Subscript[l, 1])))(Subscript[y, 1])[t] - ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[y, 2])[t] - (Subscript[y, 1])[t])-(Subscript[m, 1]) g,
(Subscript[m, 2]) (Subscript[x, 2])''[t] ==  ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[x, 2])[t] - (Subscript[x, 1])[t]),
(Subscript[m, 2]) (Subscript[y, 2])''[t] ==  ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[y, 2])[t] - (Subscript[y, 1])[t]) - (Subscript[m, 2]) g};
aeqns2 = {(Subscript[x, 1])[t]^2 + (Subscript[y, 1])[t]^2 == (Subscript[l, 1])^2, ( (Subscript[x, 2])[t]-(Subscript[x, 1])[t])^2 +((Subscript[y, 2])[t]- (Subscript[y, 1])[t])^2 == (Subscript[l, 2])^2};
ics2 = {(Subscript[x, 1])[0] ==  1, (Subscript[y, 1])[0] == -1, (Subscript[x, 1])'[0] == 0, (Subscript[y, 1])'[0] == 0, (Subscript[x, 2])[0] == 1 + Sin[.03], (Subscript[y, 2])[0] == -1 + Cos[.03], (Subscript[x, 2])'[0] == 0, (Subscript[y, 2])'[0] == 0};
params = {g -> 9.81, (Subscript[m, 1]) -> 1, (Subscript[m, 2])-> 1, (Subscript[l, 1]) -> Sqrt[2], (Subscript[l, 2]) -> 1};
soldp2 = First[NDSolve[{deqns2, aeqns2,ics2} /. params,
 {(Subscript[x, 1]), (Subscript[y, 1]), (Subscript[x, 2]), (Subscript[y, 2]),(Subscript[λ, 1]), (Subscript[λ, 2])}, {t, 0, 15},
Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}]];

Animate[
Graphics[{{PointSize[.025], {Red, Point[{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}]},
{Blue, Point[{(Subscript[x, 2])[t], (Subscript[y, 2])[t]}]}, Line[{{0,0}, {(Subscript[x, 1])[t], (Subscript[y, 1])[t]},{(Subscript[x, 2])[t], (Subscript[y, 2])[t]}}]} /. soldp2,
{Gray, Line[Map[Function[Evaluate[{(Subscript[x, 2])[#], (Subscript[y, 2])[#]} /. soldp2]], Range[0,t,0.025]]]}},
PlotRange -> {{-4,4},{-4,0}}, Axes -> True, Ticks -> False, ImageSize -> 500],
{t,0,15,.025}, SaveDefinitions -> True]

This simulates a trajectory of the double pendulum with (Subscript[l, 1]) = (Sqrt[2]), (Subscript[l, 2]) = 1, (Subscript[m, 1]) = (Subscript[m, 2]) = 1 and initial conditions (Subscript[θ, 1])(0)= (π/4), (Subscript[θ, 2])(0) = π - .05, (Subscript[θ, 1])'(0) = (Subscript[θ, 2])'(0) = 0.

deqns3 = {(Subscript[m, 1]) (Subscript[x, 1])''[t] == ((((Subscript[λ, 1])[t])/(Subscript[l, 1])))(Subscript[x, 1])[t] - ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[x, 2])[t] - (Subscript[x, 1])[t]),
(Subscript[m, 1]) (Subscript[y, 1])''[t] == ((((Subscript[λ, 1])[t])/(Subscript[l, 1])))(Subscript[y, 1])[t] - ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[y, 2])[t] - (Subscript[y, 1])[t])-(Subscript[m, 1]) g,
(Subscript[m, 2]) (Subscript[x, 2])''[t] ==  ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[x, 2])[t] - (Subscript[x, 1])[t]),
(Subscript[m, 2]) (Subscript[y, 2])''[t] ==  ((((Subscript[λ, 2])[t])/(Subscript[l, 2])))((Subscript[y, 2])[t] - (Subscript[y, 1])[t]) - (Subscript[m, 2]) g};
aeqns3 = {(Subscript[x, 1])[t]^2 + (Subscript[y, 1])[t]^2 == (Subscript[l, 1])^2, ( (Subscript[x, 2])[t]-(Subscript[x, 1])[t])^2 +((Subscript[y, 2])[t]- (Subscript[y, 1])[t])^2 == (Subscript[l, 2])^2};
ics3 = {(Subscript[x, 1])[0] ==  1, (Subscript[y, 1])[0] == -1, (Subscript[x, 1])'[0] == 0, (Subscript[y, 1])'[0] == 0, (Subscript[x, 2])[0] == 1 + Sin[.05], (Subscript[y, 2])[0] == -1 + Cos[.05], (Subscript[x, 2])'[0] == 0, (Subscript[y, 2])'[0] == 0};
params = {g -> 9.81, (Subscript[m, 1]) -> 1, (Subscript[m, 2])-> 1, (Subscript[l, 1]) -> Sqrt[2], (Subscript[l, 2]) -> 1};
soldp3 = First[NDSolve[{deqns3, aeqns3,ics3} /. params,
 {(Subscript[x, 1]), (Subscript[y, 1]), (Subscript[x, 2]), (Subscript[y, 2]),(Subscript[λ, 1]), (Subscript[λ, 2])}, {t, 0, 15},
Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}]];

Animate[
Graphics[{{PointSize[.025], {Red, Point[{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}]},
{Blue, Point[{(Subscript[x, 2])[t], (Subscript[y, 2])[t]}]}, Line[{{0,0}, {(Subscript[x, 1])[t], (Subscript[y, 1])[t]},{(Subscript[x, 2])[t], (Subscript[y, 2])[t]}}]} /. soldp3,
{Gray, Line[Map[Function[Evaluate[{(Subscript[x, 2])[#], (Subscript[y, 2])[#]} /. soldp3]], Range[0,t,0.025]]]}},
PlotRange -> {{-4,4},{-4,0}}, Axes -> True, Ticks -> False, ImageSize -> 500],
{t,0,15,.025}, SaveDefinitions -> True]

Commands for graphics in the lecture notes

ParametricPlot[
Evaluate[{{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}, {(Subscript[x, 2])[t], (Subscript[y, 2])[t]}} /. soldp], {t, 0, 15}, PlotStyle -> {{Gray}, Black}, ImageSize -> Medium, PlotLegends -> {"Mass 1", "Mass 2"}]

ParametricPlot[
Evaluate[{{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}, {(Subscript[x, 2])[t], (Subscript[y, 2])[t]}} /. soldp], {t, 0, 15}, PlotStyle -> {Red, Blue}, ImageSize -> Medium, PlotLegends -> {"Mass 1", "Mass 2"}]

ParametricPlot[
Evaluate[{{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}, {(Subscript[x, 2])[t], (Subscript[y, 2])[t]}} /. soldp], {t, 0, 15}, PlotStyle -> {{Gray}, Black}, ImageSize -> Medium, PlotLegends -> {"Mass 1", "Mass 2"}]

ParametricPlot[
Evaluate[{{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}, {(Subscript[x, 2])[t], (Subscript[y, 2])[t]}} /. soldp], {t, 0, 15}, PlotStyle -> {Red, Blue}, ImageSize -> Medium, PlotLegends -> {"Mass 1", "Mass 2"}]

ParametricPlot[
Evaluate[{{(Subscript[x, 1])[t], (Subscript[y, 1])[t]}, {(Subscript[x, 2])[t], (Subscript[y, 2])[t]}} /. soldp], {t, 0, 15}, PlotStyle -> {{Gray}, Black}, ImageSize -> Medium, PlotLegends -> {"Mass 1", "Mass 2"}]