Numerische Berechnung der Bewegung eines ebenen Doppelpendels (mit Lagrange 2. Art).

Die x-Achse geht von links nach rechts (links negativ, rechts positiv), die y-Achse von unten nach oben (unten negativ, oben positiv). Der Nullpunkt (x = 0, y = 0) befindet sich am Aufhängepunkt des ersten (oberen) Pendels.

Das erste Pendel hat die feste Pendellänge L und Masse M, das zweite die feste Pendellänge l und Masse m.
ϕ ist der Winkel zwischen erstem Pendel und negativer y-Achse, ψ ist der Winkel zwischen erstem und zweitem Pendel.

\(\(\(Needs["Mechanics`DoublePendulum`"]\)\(
\)\(
\)
    \)\)
x1 := L Sin[ϕ[t]]
y1 := -L Cos[ϕ[
            t]]

x2 := l Sin[ϕ[t] + ψ[t]]  +  x1
y2 := -l Cos[ϕ[t] + ψ[t]] + 
        y1

T1 := ((M )/2) (D[x1, t]^2 + 
            D[y1, t]^2)
V1 := M g y1

T2 := 
      (m/2)  (D[x2, t]^2 + D[y2, t]^2)
V2 := 
      m g y2  


Lagrange = 
      T1 - V1 + T2 - V2 // 
        FullSimplify

gleichung1 = 
      D[Lagrange, ϕ'[t], t] - D[Lagrange, ϕ[t]] == 0 // 
        FullSimplify

gleichung2 = 
      D[Lagrange, ψ'[t], t] - D[Lagrange, ψ[t]] == 0 // 
        FullSimplify

L := 1
l := 1
M := 1
m := 1
g := 9.81

    
Zeitspanne := 20

    
loesung1 = NDSolve[{gleichung1, 
            gleichung2, ϕ[0] ==  π, ψ[0] ==  
              (10^(-6)), ϕ'[0] == 
              0, ψ'[0] ==  0}, {ϕ[t], ψ[t]}, {t, 0, 
            Zeitspanne}, 
          MaxSteps -> 
            Zeitspanne 200][[1]]

    
loesung2 = NDSolve[{gleichung1, 
            gleichung2, ϕ[0] ==  π - (10^(-6)), ψ[0] ==  
              (10^(-6)), ϕ'[0] == 0, ψ'[0] ==  
              0}, {ϕ[t], ψ[t]}, {t, 0, Zeitspanne}, 
          MaxSteps -> Zeitspanne 200][[1]]

Plot[ϕ[t] /. loesung1, {t, 0, Zeitspanne}, 
      AxesLabel -> {t, ϕ}]
Plot[ψ[t] /. loesung1, {t, 0, Zeitspanne}, 
      AxesLabel -> {t, ψ}]

Plot[(ϕ[t] /. loesung1)  -  (ϕ[t]  /. 
            loesung2), {t, 0, Zeitspanne}, 
      AxesLabel -> {t, ϕ}]
Plot[(ψ[t] /. loesung1) - (ψ[t] /. loesung2), {t, 0, 
        Zeitspanne}, AxesLabel -> {t, ψ}]

Table[
      DoublePendulumPlot[{ϕ[t] /. loesung1, ψ[t] /. loesung1}, 
        LengthRatio -> L/l, Mass1 -> M, Mass2 -> m, 
        AspectRatio ->   1], {t, 0, Zeitspanne, 0.1}]