This sample program is developed to show the power of Mathematica OOP as shown in other my OOP projects.

It is well known that the double pendulum motion becomes a chaos with time development. Also, in the time development of a pendulum, we can observe another kind of chaos caused by the initial condition. A very small fluctuation of the initial condition affect the time development deeply. In this program case, angle perturbation less than 10^-12 can have an effect on the time development.

To observe this phenomenon this program traces simultaneously a several tens of double pendulums each has random initial angle difference. The pendulums are represented by the instance constructed from the OOP Lagrange equation of motion class. The Mathematica OOP can represent these number of double pendulums motion in time development with a simple and an effective way.

Setup for global parameters and OOP class

{g = 9.8, m = 1, r1 = 1, r2 = 0.5, time = 50};

case[nam_] := 
  Module[{\[Theta]1, \[Theta]2, ans, T1, T2, V1, V2, t, L, lkeq, 
    initcond},

   initialize[nam[th1_, th2_]] ^:= (
     (* Lagrangian setup *)
     T1 = 1/2 m*r1^2*\[Theta]1'[t]^2;
     V1 = -m*g*r1*Cos[\[Theta]1[t]];
     T2 = 
      1/2 m*(r1^2*\[Theta]1'[t]^2 + r2^2*\[Theta]2'[t]^2 + 
         2 r1*r2*\[Theta]1'[t]*\[Theta]2'[t]*r1*r2*
          Cos[\[Theta]1[t] - \[Theta]2[t]]);
     V2 = -m*g*(r1*Cos[\[Theta]1[t]] + r2*Cos[\[Theta]2[t]]);
     L = T1 + T2 - (V1 + V2);

     (* Lagrange equation of motion *)

     lkeq = {D[D[L, \[Theta]1'[t]], t] - D[L, \[Theta]1[t]] == 0,
       D[D[L, \[Theta]2'[t]], t] - D[L, \[Theta]2[t]] == 0};
     initcond = {
       \[Theta]1[0] == th1,
       \[Theta]2[0] == th2,
       \[Theta]1'[0] == 0,
       \[Theta]2'[0] == 0};

     (* Numerical solve of equation *)

     ans = NDSolve[{lkeq, initcond}, {\[Theta]1, \[Theta]2}, {t, 0, 
         time}, MaxSteps -> Infinity, PrecisionGoal -> \[Infinity]][[
       1]];
     );

   pendulum[nam[tr_]] ^:= (
     (* Pendulum graphics return *)

     Graphics[{Line[{{0, 
          0}, {r1*Sin[\[Theta]1[tr]], -r1*Cos[\[Theta]1[tr]]} /. 
          ans}], Line[{{r1*Sin[\[Theta]1[tr]], -r1*
            Cos[\[Theta]1[tr]]}, {(r1*Sin[\[Theta]1[tr]] + 
             r2*Sin[\[Theta]2[tr]]), (-r1*Cos[\[Theta]1[tr]] - 
             r2*Cos[\[Theta]2[tr]])}} /. ans]}, 
      PlotRange -> {{-1.8, 1.8}, {-1.8, 1.8}}]
     )
   ];

Setup initial conditions and construct instances, and display of results

{pendulums = 50, angle1 = 4 Pi/3, angle2 = 4 Pi/3, butterfly = 10^-12};

objectList = Table[Unique[], {pendulums}];
Map[case[#] &, objectList];
Map[initialize[#[angle1, angle2 + butterfly*RandomReal[{-1, 1}]]] &, 
  objectList];

Animate[Show[
  Map[pendulum[#[tr]] &, objectList]
  ], {tr, 0, time}, AnimationRepetitions -> 1, AnimationRate -> 1]

Enjoy a sudden divergence.