Clear[“Global`*”];

(*定义参数*)
{g, l1, l2, m1, m2, tt} = {9.8, 1, 1, 1, 3, 100};

(*定义拉格朗日方程*)
L = (m1 + m2)/2 l1^2 (a'[t])^2 + m2/2 l2^2 (b'[t])^2 + 
   m2*l1*l2*a'[t]*b'[t]*Cos[a[t] – b[t]] + (m1 + m2)*g*l1*Cos[a[t]] + 
   m2*g*l2*Cos[b[t]];

(*定义微分方程以及合适的边界条件*)
funcs = {D[L, a[t]] == D[D[L, a'[t]], t], 
   D[L, b[t]] == D[D[L, b'[t]], t], a[0] == 3 Pi/4, b[0] == Pi/6, 
   a'[0] == 0, b'[0] == 0};

(*微分方程的100精度数字解*)
sov = NDSolve[funcs, {a, b}, {t, 0, tt}, PrecisionGoal -> 100][[1]];

(*随时间变化的角度*)
a = a /. sov;
b = b /. sov;

(*整个系统的动量，机械能*)
T[t_] := (m1 + m2)/2 l1^2 ((D[a[x], x]) /. x -> t)^2 + 
   m2/2 l2^2 ((D[b[x], x]) /. x -> t)^2 + 
   m2*l1*l2*((D[a[x], x]) /. x -> t)*((D[b[x], x]) /. x -> t)*
    Cos[a[t] – b[t]];
V[t_] := -(m1 + m2)*g*l1*Cos[a[t]] – m2*g*l2*Cos[b[t]];
Energy[t_] := T[t] + V[t];

(*m1和m2的坐标*)
p1[t_] := {l1*Sin[a[t]], -l1*Cos[a[t]]};
p2[t_] := {l1*Sin[a[t]] + l2*Sin[b[t]], -l1*Cos[a[t]] – l2*Cos[b[t]]};

(*可视化模块*)
result = Animate[{Show[
    Graphics[Circle[{0, 0}, 1], PlotRange -> {{-2, 2}, {-2, 2}}, 
     AxesOrigin -> {0, 0}, Axes -> True, Frame -> False, 
     ImageSize -> Medium], 
    Graphics[{Line[{{0, 0}, p1[u], p2[u]}], Disk[p1[u], m1/15], 
      Disk[p2[u], m2/15]}], 
    If[u < 3, 
     ParametricPlot[p2[t], {t, 0, u}, 
      ColorFunction -> 
       Function[{x, y, u}, f = Hue[0.58, 0.58, 0.8, u]]], 
     ParametricPlot[p2[t], {t, u – 3, u}, 
      ColorFunction -> 
       Function[{x, y, u}, f = Hue[0.58, 0.58, 0.8, u]]]]], 
   ParametricPlot[{{t, T[t]}, {t, V[t]}, {t, Energy[t]}}, {t, 0, u}, 
    ImageSize -> Medium, AspectRatio -> 0.75, 
    PlotRange -> {{0, u}, {-25, 30}}]}, {u, 0.01, tt, 0.01}, 
  AnimationDirection -> Forward, RefreshRate -> 180, 
  DefaultDuration -> tt]