二重振り子のカオス（アニメーション）

Clear["Global'*"];

○初期値の設定
ｌ：腕の長さ　ｍ：おもりの質量　ｇ：重力加速度

l1=1.0;
l2=1.0;
m1=1;
m2=1;
g=9.8;

l3=1.0;
l4=1.0;
m3=1;
m4=1;

θshoki：θの初期値　θshokiD：θの導関数

Θ1shoki=0.9Pi;
Θ1shokiD =0;
Θ2shoki=0.9Pi;
Θ2shokiD=0;

Θ3shoki=0.9Pi;
Θ3shokiD =0;
Θ4shoki=0.9000000001Pi;                   (*θ２の初期値と1/(10^10)のずれを設定した。*)
Θ4shokiD=0;

○運動方程式の解

１つ目の系に関する運動方程式

solution2 = NDSolve[
{(m1+m2)*g*Sin[Θ1[t]]+m2*l2*Sin[Θ1[t]-Θ2[t]]*D[Θ2[t],{t,1}]*D[Θ2[t],{t,1}]+(m1+m2)*l1*D[Θ1[t],{t,2}]+m2*l2*D[Θ2[t],{t,2}]*Cos[Θ1[t]-Θ2[t]]==0,
g*Sin[Θ2[t]]-l1*D[Θ1[t],{t,1}]*D[Θ1[t],{t,1}]*Sin[Θ1[t]-Θ2[t]]+l1*D[Θ1[t],{t,2}]*Cos[Θ1[t]-Θ2[t]]+l2*D[Θ2[t],{t,2}]==0,
Θ1[0]==Θ1shoki,Θ1'[0]==Θ1shokiD,Θ2[0]==Θ2shoki,Θ2'[0]==Θ2shokiD}
,
{Θ1,Θ2}
,
{t,0,40}
,
MaxSteps->∞
]

２つ目の系に関する運動方程式

solution3 = NDSolve[
{(m3+m4)*g*Sin[Θ3[t]]+m4*l4*Sin[Θ3[t]-Θ4[t]]*D[Θ4[t],{t,1}]*D[Θ4[t],{t,1}]+(m3+m4)*l3*D[Θ3[t],{t,2}]+m4*l4*D[Θ4[t],{t,2}]*Cos[Θ3[t]-Θ4[t]]==0,
g*Sin[Θ4[t]]-l1*D[Θ3[t],{t,1}]*D[Θ3[t],{t,1}]*Sin[Θ3[t]-Θ4[t]]+l1*D[Θ3[t],{t,2}]*Cos[Θ3[t]-Θ4[t]]+l2*D[Θ4[t],{t,2}]==0,
Θ3[0]==Θ3shoki,Θ3'[0]==Θ3shokiD,Θ4[0]==Θ4shoki,Θ4'[0]==Θ4shokiD}
,
{Θ3,Θ4}
,
{t,0,40}
,
MaxSteps->∞
]

○アニメーションの実装

Animate[
Show[
ListLinePlot[{
{0,0}/.solution2[[1]],

{l1 Sin[Θ1[t]],-l1 Cos[Θ1[t]]}/.solution2[[1]],
{l1 Sin[Θ1[t]]+l2 Sin[Θ2[t]],-l1 Cos[Θ1[t]]-l2 Cos[Θ2[t]]}/.solution2[[1]]
},
PlotRange->{{-2,2},{-2,2}},PlotStyle->Hue[.1],AspectRatio->1]
,
ListLinePlot[{
{0,0}/.solution3[[1]],

{l3 Sin[Θ3[t]],-l3 Cos[Θ3[t]]}/.solution3[[1]],
{l3 Sin[Θ3[t]]+l4 Sin[Θ4[t]],-l3 Cos[Θ3[t]]-l4 Cos[Θ4[t]]}/.solution3[[1]]
},
PlotRange->{{-2,2},{-2,2}},PlotStyle->Hue[.6] ,AspectRatio->1]
,
ListPlot[{
{l1 Sin[Θ1[t]],-l2 Cos[Θ1[t]]}/.solution2[[1]],
{l1 Sin[Θ1[t]]+l2 Sin[Θ2[t]],-l1 Cos[Θ1[t]]-l2 Cos[Θ2[t]]}/.solution2[[1]]},
 PlotStyle->PointSize[0.05],PlotStyle->Hue[.1]]
,
ListPlot[{
{l3 Sin[Θ3[t]],-l4 Cos[Θ3[t]]}/.solution3[[1]],
{l3 Sin[Θ3[t]]+l4 Sin[Θ4[t]],-l3 Cos[Θ3[t]]-l4 Cos[Θ4[t]]}/.solution3[[1]]},
 PlotStyle->PointSize[0.05],PlotStyle->Hue[.6]]
]
,
{t,0,40,0.005}
,
DefaultDuration->40
,
AnimationRunning->False
]