Double Pendulum
Model the motion of a double pendulum in Cartesian coordinates.
deqns = {Subscript[m, 1]x1''[t] == ((λ1[t]/Subscript[l, 1])) x1[t] - ((λ2[t]/Subscript[l, 2]))(x2[t] - x1[t]), Subscript[m, 1]y1''[t] == ((λ1[t]/Subscript[l, 1])) y1[t] - ((λ2[t]/Subscript[l, 2]))(y2[t] - y1[t]) - Subscript[m, 1] g,
Subscript[m, 2]x2''[t] == ((λ2[t]/Subscript[l, 2]))(x2[t] - x1[t]),
Subscript[m, 2]y2''[t] == ((λ2[t]/Subscript[l, 2])) (y2[t] - y1[t]) - Subscript[m, 2]g};aeqns = {x1[t]^2 + y1[t]^2 == Subscript[l, 1]^2, (x2[t] - x1[t])^2 + (y2[t] - y1[t])^2 == Subscript[l, 2]^2};ics = {x1[0] == 1, y1[0] == 0, y1'[0] == 0, x2[0] == 1, y2[0] == -1, y2'[0] == 0};params = {g -> 9.81, Subscript[m, 1] -> 1, Subscript[m, 2] -> 1, Subscript[l, 1] -> 1, Subscript[l, 2] -> 1};soldp = First[NDSolve[{deqns, aeqns, ics} /. params, {x1, y1, x2, y2, λ1, λ2}, {t, 0, 15}, Method -> {"IndexReduction" -> {"Pantelides", "ConstraintMethod" -> "Projection"}}]];Animate[Graphics[{{PointSize[.025], {Red, Point[{x1[t], y1[t]}]}, {Blue, Point[{x2[t], y2[t]}]}, Line[{{0, 0}, {x1[t], y1[t]}, {x2[t], y2[t]}}]} /. soldp,
{Gray, Line[Map[Function[Evaluate[{x2[#], y2[#]} /. soldp]], Range[0, t, 0.025]]]}}, PlotRange -> {{-2, 2}, {-2, 0}}, Axes -> True, Ticks -> False, ImageSize -> 500], {t, 0, 10, .025}, SaveDefinitions -> True]