Return to Answer

Needs["VariationalMethods`"] Clear[s1, s2, ϕ1, ϕ2, t, g, m1, m2]; variables = {ϕ1[t], ϕ2[t]}; r1 = s1 {Sin[ϕ1[t]], -Cos[ϕ1[t]]}; r2 = r1 + s2 {Sin[ϕ2[t]], -Cos[ϕ2[t]]}; lagrangian = m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - g {0, 1}.(m1 r1 + m2 r2); eqs = EulerEquations[lagrangian, variables, t]; s1 = 1; s2 = 1/2; m1 = 2; m2 = 1; g = 9.81; tMax = 10; initial = {ϕ1[0] == Pi/2, ϕ2[0] == Pi, ϕ1'[0] == Pi/2, ϕ2'[0] == 0}; sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]]; nFrames = 100; frames = Table[Graphics[{Gray, Thick, Line[{{0, 0}, r1, r2}], Darker[Blue], Disk[{0, 0}, .1], Darker[Red], Disk[r1, m1/10], Disk[r2, m2/10]} /. sol, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}, Background -> Lighter[Cyan]], {t, 0, tMax, tMax/nFrames}]; ListAnimate[frames] 

Note that the differential equations we're solving are never written down by hand anywhere. They are automatically generated by EulerEquations.

Needs["VariationalMethods`"] Clear[s1, s2, ϕ1, ϕ2, t, g, m1, m2]; variables = {ϕ1[t], ϕ2[t]}; r1 = s1 {Sin[ϕ1[t]], -Cos[ϕ1[t]]}; r2 = r1 + s2 {Sin[ϕ2[t]], -Cos[ϕ2[t]]}; lagrangian = m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - g {0, 1}.(m1 r1 + m2 r2); eqs = EulerEquations[lagrangian, variables, t]; s1 = 1; s2 = 1/2; m1 = 2; m2 = 1; g = 9.81; tMax = 10; initial = {ϕ1[0] == Pi/2, ϕ2[0] == Pi, ϕ1'[0] == Pi/2, ϕ2'[0] == 0}; sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]]; nFrames = 100; frames = Table[Graphics[{Gray, Thick, Line[{{0, 0}, r1, r2}], Darker[Blue], Disk[{0, 0}, .1], Darker[Red], Disk[r1, m1/10], Disk[r2, m2/10]} /. sol, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}, Background -> Lighter[Cyan]], {t, tMax/nFrames, tMax, tMax/nFrames}]; ListAnimate[frames] 

Note that the differential equations we're solving are never written down by hand anywhere. They are automatically generated by EulerEquations.

Edit

As requested in the comment, the solution above can also be re-used to plot a trace of the trajectory together with the pendulum graphic:

traceFrames = Map[ParametricPlot[r2 /. sol, {t, 0, #}, PlotStyle -> Directive[Thick, Lighter[Purple]], PlotPoints -> Floor[(2 nFrames #)/tMax], PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}, Background -> Lighter[Cyan]] &, Range[tMax/nFrames, tMax, tMax/nFrames]]; ListAnimate[MapThread[Show, {traceFrames, frames}]] 

Here I re-use the list of frames from above and Show it together with a new list traceFrames in which the ParametricPlot of the end point r2 is shown up to the intermediate time corresponding to the instants used in frame. Since ParametricPlot uses interpolation, we have to make sure that the plots for long times don't lose details, by increasing the number of PlotPoints with time.

Needs["VariationalMethods`"] Clear[s1, s2, ϕ1, ϕ2, t, g, m1, m2]; variables = {ϕ1[t], ϕ2[t]}; r1 = s1 {Sin[ϕ1[t]], -Cos[ϕ1[t]]}; r2 = r1 + s2 {Sin[ϕ2[t]], -Cos[ϕ2[t]]}; lagrangian = m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - g {0, 1}.(m1 r1 + m2 r2); eqs = EulerEquations[lagrangian, variables, t]; s1 = 1; s2 = 1/2; m1 = 2; m2 = 1; g = 9.81; tMax = 10; initial = {ϕ1[0] == Pi/2, ϕ2[0] == Pi, ϕ1'[ 0] == Pi/2, ϕ2'[0] == 0};  sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]]; nFrames = 100; frames = Table[Graphics[{Gray, Thick, Line[{{0, 0}, r1, r2}], Darker[Blue], Disk[{0, 0}, .1], Darker[Red], Disk[r1, m1/10], Disk[r2, m2/10]} /. sol, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}], {t, 0, tMax,  tMax/nFrames}]; 

Needs["VariationalMethods`"] Clear[s1, s2, ϕ1, ϕ2, t, g, m1, m2]; variables = {ϕ1[t], ϕ2[t]}; r1 = s1 {Sin[ϕ1[t]], -Cos[ϕ1[t]]}; r2 = r1 + s2 {Sin[ϕ2[t]], -Cos[ϕ2[t]]}; lagrangian = m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - g {0, 1}.(m1 r1 + m2 r2); eqs = EulerEquations[lagrangian, variables, t]; s1 = 1; s2 = 1/2; m1 = 2; m2 = 1; g = 9.81; tMax = 10; initial = {ϕ1[0] == Pi/2, ϕ2[0] == Pi,  ϕ1'[0] == Pi/2, ϕ2'[0] == 0}; sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]]; nFrames = 100; frames = Table[Graphics[{Gray, Thick, Line[{{0, 0}, r1, r2}], Darker[Blue], Disk[{0, 0}, .1], Darker[Red], Disk[r1, m1/10], Disk[r2, m2/10]} /. sol, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}, Background -> Lighter[Cyan]],  {t, 0, tMax, tMax/nFrames}]; ListAnimate[frames] 

Note that the differential equations we're solving are never written down by hand anywhere. They are automatically generated by EulerEquations.

Edit: more general solution

To get more out of this animation, it may be of interest to change the simulation so that it can handle different masses, and different pendulum lengths. But this means you need to understand how the differential equation must be modified. For many constrained mechanics problems, including the double pendulum, the Lagrange formalism is the most efficient way to set up the equations of motion. Mathematica has a VariationalMethods package that helps to automate most of the steps. This is extremely useful for solving mechanics problems because it lets you focus on the physics while taking over the mathematical work of getting from the physics to the equations of motion.

Here is how it works. The following few lines contain all the steps from defining the variables to establishing and then solving the equations of motion, with arbitrary parameters for the pendulum geometry. Here s1, s2 are the lengths of the two rods connecting the masses m1, m2. Most of the following equations can be done entirely in terms of the vectorial positions r1 and r2 without thinking about their components too much. The actual variables we solve for are ϕ1, ϕ2. To show the effect of doubling one of the masses and halving one of the pendulum rods, I then set those variables before calling NDSolve:

Needs["VariationalMethods`"] Clear[s1, s2, ϕ1, ϕ2, t, g, m1, m2]; variables = {ϕ1[t], ϕ2[t]}; r1 = s1 {Sin[ϕ1[t]], -Cos[ϕ1[t]]}; r2 = r1 + s2 {Sin[ϕ2[t]], -Cos[ϕ2[t]]}; lagrangian = m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - g {0, 1}.(m1 r1 + m2 r2); eqs = EulerEquations[lagrangian, variables, t]; s1 = 1; s2 = 1/2; m1 = 2; m2 = 1; g = 9.81; tMax = 10; initial = {ϕ1[0] == Pi/2, ϕ2[0] == Pi, ϕ1'[ 0] == Pi/2, ϕ2'[0] == 0}; sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]]; nFrames = 100; frames = Table[Graphics[{Gray, Thick, Line[{{0, 0}, r1, r2}], Darker[Blue], Disk[{0, 0}, .1], Darker[Red], Disk[r1, m1/10], Disk[r2, m2/10]} /. sol, PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}], {t, 0, tMax, tMax/nFrames}];