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

Die kartesischen Koordinaten ausgedrueckt in ϕ und ψ:

xM = L*Sin[ϕ[t]];

yM = -L*Cos[ϕ[t]];

xm = xM + l*Sin[ϕ[t] + ψ[t]];

ym = yM - l*Cos[ϕ[t] + ψ[t]];

Kinetische und potentielle Energie:

TM = 1/2*M*L^2*ϕ'[t]^2;

Tm = 1/2*m*(D[xm, t]^2 + D[ym, t]^2);

VM = M*g*yM;

Vm = m*g*ym;

Die Lagrange Funktion:

Lag[ϕ[t] _, ϕ'[t] _, ψ[t] _, ψ'[t] _] := 
        TM + Tm - (VM + Vm);

Die beiden Bewegungsgleichungen:

equ1 = 
        D[D[Lag[ϕ[t] _, ϕ'[t] _, ψ[t] _, ψ'[
                      t] _], ϕ'[t]], t] - 
            D[Lag[ϕ[t] _, ϕ'[t] _, ψ[t] _, ψ'[
                    t] _], ϕ[t]] == 0;

equ2 = 
        D[D[Lag[ϕ[t] _, ϕ'[t] _, ψ[t] _, ψ'[
                      t] _], ψ'[t]], t] - 
            D[Lag[ϕ[t] _, ϕ'[t] _, ψ[t] _, ψ'[
                    t] _], ψ[t]] == 0;

Bestimmung der numerischen Loesung mit NDSolve:

Loesung1 = NDSolve[{equ1, 
              equ2, ϕ[0] == Pi, ψ[0] == 
                10^(-6), ϕ'[0] == 
                0, ψ'[0] == 0}, {ϕ, ψ}, {t, 0, 20}, 
            MaxSteps -> 10000][[1]];

Loesung2 = NDSolve[{equ1, 
              equ2, ϕ[0] == 
                Pi - 10^(-6), ψ[0] == 
                10^(-6), ϕ'[0] == 
                0, ψ'[0] == 0}, {ϕ, ψ}, {t, 0, 20}, 
            MaxSteps -> 10000][[1]];

Die Plots:

Plot[Evaluate[ϕ[t] /. Loesung1], {t, 0, 20}];

Plot[Evaluate[ψ[t] /. Loesung1], {t, 0, 20}];

Plot[Evaluate[ϕ[t] /. Loesung2], {t, 0, 20}];

Plot[Evaluate[ψ[t] /. Loesung2], {t, 0, 20}];

Die Differenzen der Winkel zu minimal verschiedenen Anfangsbedingungen:

Plot[
        Evaluate[(ϕ[t] /. Loesung1) - (ϕ[t] /. 
                Loesung2)], {t, 0, 20}, 
        PlotRange -> All];

Plot[
        Evaluate[(ψ[t] /. Loesung1) - (ψ[t] /. 
                Loesung2)], {t, 0, 20}, 
        PlotRange -> All];

Man sieht,  
    dass sie kleine Unterschiede in den Anfangsbedingungen schon nach kurzer Zeit sehr stark auswirken!

Ein paar nette Graphiken,  
    die Mathematica auch annimieren kann:  ->  Cell  ->  Animate Selected Graphics

Needs["Mechanics`DoublePendulum`"]

movie1 = 
        Table[DoublePendulumPlot[{ϕ[t] /. Loesung1, ψ[t] /. 
                Loesung1}], {t, 0, 20,  .1}];

movie2 = 
        Table[DoublePendulumPlot[{ϕ[t] /. Loesung2, ψ[t] /. 
                Loesung2}], {t, 0, 20,  .1}];

Man kann es auch als animated gif exportieren.

Export["dp1.gif", movie1, "GIF"]

Export["dp2.gif", movie2, "GIF"]