Chaos in double mathematical pendulum

One set of initial conditions

tmax=20;
m1=1;
m2=2;
g=9.80665;

l1=1;
l2=0.9;

init1=π;
init2=π+1;
initpri1=1;
initpri2=-2;

integration=NDSolve[{
(m1+m2) l1 θ1''[t]+m2 l2 θ2''[t] Cos[θ1[t]-θ2[t]]+m2 l2 θ2'[t]^2 Sin[θ1[t]-θ2[t]]+g(m1+m2)Sin[θ1[t]]==0,m2 l2 θ2''[t]+m2 l1 θ1''[t]Cos[θ1[t]-θ2[t]]-m2 l1 θ1'[t]^2 Sin[θ1[t]-θ2[t]]+m2 g Sin[θ2[t]]==0,
θ1[0]==init1,
θ2[0]==init2,
θ1'[0]==initpri1,
θ2'[0]==initpri2},
{θ1,θ2},{t,0,tmax},MaxSteps->Infinity, PrecisionGoal->10];
plot1=ParametricPlot[Evaluate[{θ1[t],θ2[t]}/.integration],{t,0,tmax},FrameLabel->{θ1[t],θ2[t]},ImageSize->200,LabelStyle->18,Frame->True];
plot2=Plot[Evaluate[θ1[t]/.integration],{t,0,tmax},AxesLabel-> {t,θ1[t]},PlotRange-> All,ImageSize->260,LabelStyle->18];
plot3=Plot[Evaluate[θ2[t]/.integration],{t,0,tmax},AxesLabel-> {t,θ2[t]},PlotRange-> All,ImageSize->260,LabelStyle->18];
Row[{plot1,plot2,plot3}]
(*PROJECTIONS*)
plo1=ParametricPlot[Evaluate[{Sin[θ1[t]],Sin[θ2[t]]}/.integration],{t,0,tmax},Frame-> True, FrameLabel->{Sin[θ1[t]],Sin[θ2[t]]},ImageSize->240,LabelStyle->18];
plo2=Plot[Evaluate[Sin[θ1[t]]/.integration],{t,0,tmax},AxesLabel-> {t,Sin[θ1[t]]},PlotRange-> All,ImageSize->260,LabelStyle->18];
plo3=Plot[Evaluate[Sin[θ2[t]]/.integration],{t,0,tmax},AxesLabel-> {t,Sin[θ2[t]]},PlotRange-> All,ImageSize->260,LabelStyle->18];
Row[{plo1,plo2,plo3}]
pplo1=ParametricPlot[Evaluate[{
Sin[θ1[t]],
D[Sin[θ1[t]],t]}/.integration],{t,0,tmax},
FrameLabel-> {"sin θ1","(∂_t)sin θ1"},PlotLabel->Style["Phase plot of sin θ1",20],PlotRange-> All,ImageSize->330,AspectRatio->1,LabelStyle->18,Frame->True];
pplo2=ParametricPlot[Evaluate[{
Sin[θ2[t]],
D[Sin[θ2[t]],t]}/.integration],{t,0,tmax},
FrameLabel-> {"sin θ2","(∂_t)sin θ2"},PlotLabel->Style["Phase plot of sin θ2",20],PlotRange-> All,ImageSize->330,AspectRatio->1,LabelStyle->18,Frame->True];
Row[{pplo1,pplo2}]

Two sets of initial conditions (chaos demonstration)

tmax=20;
m1=1;
m2=2;
g=9.80665;

l1=1;
l2=0.9;

init1=π;
init2=π+1;
initpri1=1;
initpri2=-2;

dinit1=0.00001;
dinit2=0.0;
dinitpri1=0.0;
dinitpri2=0.0;

integration1=NDSolve[{
(m1+m2) l1 θ1''[t]+m2 l2 θ2''[t] Cos[θ1[t]-θ2[t]]+m2 l2 θ2'[t]^2 Sin[θ1[t]-θ2[t]]+g(m1+m2)Sin[θ1[t]]==0,m2 l2 θ2''[t]+m2 l1 θ1''[t]Cos[θ1[t]-θ2[t]]-m2 l1 θ1'[t]^2 Sin[θ1[t]-θ2[t]]+m2 g Sin[θ2[t]]==0,
θ1[0]==init1,
θ2[0]==init2,
θ1'[0]==initpri1,
θ2'[0]==initpri2},
{θ1,θ2},{t,0,tmax},MaxSteps->Infinity, PrecisionGoal->10];
integration2=NDSolve[{
(m1+m2) l1 θ1''[t]+m2 l2 θ2''[t] Cos[θ1[t]-θ2[t]]+m2 l2 θ2'[t]^2 Sin[θ1[t]-θ2[t]]+g(m1+m2)Sin[θ1[t]]==0,m2 l2 θ2''[t]+m2 l1 θ1''[t]Cos[θ1[t]-θ2[t]]-m2 l1 θ1'[t]^2 Sin[θ1[t]-θ2[t]]+m2 g Sin[θ2[t]]==0,
θ1[0]==init1+dinit1,
θ2[0]==init2+dinit2,
θ1'[0]==initpri1+dinitpri1,
θ2'[0]==initpri2+dinitpri2},
{θ1,θ2},{t,0,tmax},MaxSteps->Infinity, PrecisionGoal->10];
plot1=ParametricPlot[Evaluate[{θ1[t],θ2[t]}/.integration1],{t,0,tmax},FrameLabel->{θ1[t],θ2[t]},ImageSize->260,LabelStyle->18,Frame->True];
plot12=ParametricPlot[Evaluate[{θ1[t],θ2[t]}/.integration2],{t,0,tmax},FrameLabel->{θ1[t],θ2[t]},ImageSize->260,PlotStyle->Red,LabelStyle->18,Frame->True];
plot2=Plot[Evaluate[θ1[t]/.integration1],{t,0,tmax},AxesLabel-> {t,θ1[t]},PlotRange-> All,ImageSize->260,LabelStyle->18];
plot22=Plot[Evaluate[θ1[t]/.integration2],{t,0,tmax},AxesLabel-> {t,θ1[t]},PlotStyle->Red,PlotRange-> All,ImageSize->220];
plot3=Plot[Evaluate[θ2[t]/.integration1],{t,0,tmax},AxesLabel-> {t,θ2[t]},PlotRange-> All,ImageSize->260,LabelStyle->18];
plot32=Plot[Evaluate[θ2[t]/.integration2],{t,0,tmax},AxesLabel-> {t,θ2[t]},PlotStyle->Red,PlotRange-> All,ImageSize->260,LabelStyle->18];
Row[{Show[{plot1,plot12}],Show[{plot2,plot22}],Show[{plot3,plot32}]}]
(*PROJECTIONS*)
plo1=ParametricPlot[Evaluate[{Sin[θ1[t]],Sin[θ2[t]]}/.integration1],{t,0,tmax},FrameLabel->{Sin[θ1[t]],Sin[θ2[t]]},ImageSize->260,LabelStyle->18,Frame->True];
plo12=ParametricPlot[Evaluate[{Sin[θ1[t]],Sin[θ2[t]]}/.integration2],{t,0,tmax},FrameLabel->{Sin[θ1[t]],Sin[θ2[t]]},PlotStyle->{Dashed,Red},ImageSize->260,LabelStyle->18,Frame->True];
plo2=Plot[Evaluate[Sin[θ1[t]]/.integration1],{t,0,tmax},AxesLabel-> {t,Sin[θ1[t]]},PlotRange-> All,ImageSize->260,LabelStyle->18];
plo22=Plot[Evaluate[Sin[θ1[t]]/.integration2],{t,0,tmax},AxesLabel-> {t,Sin[θ1[t]]},PlotStyle->Red,PlotRange-> All,ImageSize->260,LabelStyle->18];
plo3=Plot[Evaluate[Sin[θ2[t]]/.integration1],{t,0,tmax},AxesLabel-> {t,Sin[θ2[t]]},PlotRange-> All,ImageSize->260,LabelStyle->18];
plo32=Plot[Evaluate[Sin[θ2[t]]/.integration2],{t,0,tmax},AxesLabel-> {t,Sin[θ2[t]]},PlotStyle->Red,PlotRange-> All,ImageSize->260,LabelStyle->18];
Row[{Show[{plo1,plo12}],Show[{plo2,plo22}],Show[{plo3,plo32}]}]
pplo1=ParametricPlot[Evaluate[{
Sin[θ1[t]],
D[Sin[θ1[t]],t]}/.integration1],{t,0,tmax},
FrameLabel-> {"sin θ1","(∂_t)sin θ1"},PlotLabel->Style["Phase plot of sin θ1",20],PlotRange-> All,ImageSize->330,AspectRatio->1,LabelStyle->18,Frame->True];
dpplo1=ParametricPlot[Evaluate[{
Sin[θ1[t]],
D[Sin[θ1[t]],t]}/.integration2],{t,0,tmax},
FrameLabel-> {"sin θ1","(∂_t)sin θ1"},PlotRange-> All,ImageSize->330,AspectRatio->1,PlotStyle->{Red,Dashed},LabelStyle->18,Frame->True];
pplo2=ParametricPlot[Evaluate[{
Sin[θ2[t]],
D[Sin[θ2[t]],t]}/.integration1],{t,0,tmax},
FrameLabel-> {"sin θ2","(∂_t)sin θ2"},PlotLabel->Style["Phase plot of sin θ2",20],PlotRange-> All,ImageSize->330,AspectRatio->1,LabelStyle->18,Frame->True];
dpplo2=ParametricPlot[Evaluate[{
Sin[θ2[t]],
D[Sin[θ2[t]],t]}/.integration2],{t,0,tmax},
FrameLabel-> {"sin θ2","(∂_t)sin θ2"},PlotRange-> All,ImageSize->330,AspectRatio->1,PlotStyle->{Red,Dashed},LabelStyle->18,Frame->True];
Row[{Show[{pplo1,dpplo1}],Show[{pplo2,dpplo2}]}]

Creation of this notebook was supported by IT Academy program of Information Technology Foundation for Education (HITSA).