Gravitational planar three-body problem

Global initial conditions and parameters

Tmax=4.5;mxlimit=-0.5;xlimit=1.5;mylimit=-0.5;ylimit=1.1;
m1=1.0;m2=2.6;m3=1;(*no gravitational constant G=1*)

dx10=0.002;(*difference in initial coordinate*)
x10=0;y10=1;
x20=1;y20=0;
x30=1;y30=1;

vx10=0.0;vy10=-0.31;
vx20=-0.1;vy20=0.08;
vx30=-0.1;vy30=-0.11;

One set of initial conditions

Clear[t,tmax]
Manipulate[
integrate=NDSolve[{
x1'[t]==vx1[t],y1'[t]==vy1[t],
x2'[t]==vx2[t],y2'[t]==vy2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m2 vx2'[t]==((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m2 vy2'[t]==((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vx3'[t]==((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vy3'[t]==((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),
x1[0]==x10,y1[0]==y10,
x2[0]==x20,y2[0]==y20,
x3[0]==x30,y3[0]==y30,
vx1[0]==vx10,vy1[0]==vy10,
vx2[0]==vx20,vy2[0]==vy20,
vx3[0]==vx30,vy3[0]==vy30},
{x1,x2,x3,y1,y2,y3,vx1,vx2,vx3,vy1,vy2,vy3},{t,0,tmax},MaxSteps->Infinity,PrecisionGoal->9];(*solved for positions and velocities*)
initpos1=Graphics[{PointSize[0.02 m1],Blue,Point[{x10,y10}]}];
initpos2=Graphics[{PointSize[0.02 m2],Black,Point[{x20,y20}]}];
initpos3=Graphics[{PointSize[0.02 m3],Red,Point[{x30,y30}]}];
mass1=ParametricPlot[Evaluate[{x1[t],y1[t]}/.integrate],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotLabel->Style["Mass 1",20],PlotRange->All,ImageSize->200,LabelStyle->18,Frame-> True,AspectRatio->1,PerformanceGoal->"Quality"];
mass2=ParametricPlot[Evaluate[{x2[t],y2[t]}/.integrate],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotLabel->Style["Mass 2",20],PlotRange->All,PlotStyle->Black,ImageSize->200,LabelStyle->18,Frame-> True,AspectRatio->1,PerformanceGoal->"Quality"];
mass3=ParametricPlot[Evaluate[{x3[t],y3[t]}/.integrate],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotLabel->Style["Mass 3",20],PlotRange->All,PlotStyle->Red,ImageSize->200,LabelStyle->18,Frame-> True,AspectRatio->1,PerformanceGoal->"Quality"];
traj=Show[{mass1,initpos1,mass2,initpos2,mass3,initpos3},PlotLabel->Style["Positions of all masses together",20],ImageSize->500,LabelStyle->18,AspectRatio->1];
Grid[{{traj,SpanFromLeft},
{Show[mass1,initpos1],Show[mass2,initpos2],Show[mass3,initpos3] }},Frame-> All],
{{tmax,0.15 Tmax,"Time"},0.001,Tmax,Appearance->"Labeled"},SaveDefinitions->True,TrackedSymbols:>{tmax}]

One set of initial conditions, animation

Clear[t,tmax]
Animate[
integrate=NDSolve[{
x1'[t]==vx1[t],y1'[t]==vy1[t],
x2'[t]==vx2[t],y2'[t]==vy2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m2 vx2'[t]==((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m2 vy2'[t]==((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vx3'[t]==((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vy3'[t]==((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),
x1[0]==x10,y1[0]==y10,
x2[0]==x20,y2[0]==y20,
x3[0]==x30,y3[0]==y30,
vx1[0]==vx10,vy1[0]==vy10,
vx2[0]==vx20,vy2[0]==vy20,
vx3[0]==vx30,vy3[0]==vy30},
{x1,x2,x3,y1,y2,y3},{t,0,tmax},MaxSteps->Infinity,PrecisionGoal->9];(*solved for positions only*)
pointm1=Graphics[{PointSize[0.02 m1],Blue,Point[Dynamic[Evaluate[{x1[tmax],y1[tmax]}/.integrate]]]},
PlotRange->{{mxlimit,xlimit},{mylimit,ylimit}},Axes->True,FrameLabel-> {x[t],y[t]},Frame->True];
pointm2=Graphics[{PointSize[0.02 m2],Black,Point[Dynamic[Evaluate[{x2[tmax],y2[tmax]}/.integrate]]]},
PlotRange->{{mxlimit,xlimit},{mylimit,ylimit}},Axes->True,FrameLabel-> {x[t],y[t]},Frame->True];
pointm3=Graphics[{PointSize[0.02 m3],Red,Point[Dynamic[Evaluate[{x3[tmax],y3[tmax]}/.integrate]]]},
PlotRange->{{mxlimit,xlimit},{mylimit,ylimit}},Axes->True,FrameLabel-> {x[t],y[t]},Frame->True];
Show[{pointm1,pointm2,pointm3},PlotLabel->Style[ "Position of masses",20],ImageSize->500,LabelStyle-> 18],
{tmax,0,Tmax},AnimationRunning->False,AnimationRate->0.3,SaveDefinitions->True,TrackedSymbols:>{tmax}]

Animation with tails

Clear[t,tmax]
Animate[
integrate=NDSolve[{
x1'[t]==vx1[t],y1'[t]==vy1[t],
x2'[t]==vx2[t],y2'[t]==vy2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m2 vx2'[t]==((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m2 vy2'[t]==((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vx3'[t]==((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vy3'[t]==((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),
x1[0]==x10,y1[0]==y10,
x2[0]==x20,y2[0]==y20,
x3[0]==x30,y3[0]==y30,
vx1[0]==vx10,vy1[0]==vy10,
vx2[0]==vx20,vy2[0]==vy20,
vx3[0]==vx30,vy3[0]==vy30},
{x1,x2,x3,y1,y2,y3},{t,0,tmax},MaxSteps->Infinity,PrecisionGoal->9];(*solved for positions only*)
mass1=ParametricPlot[Evaluate[{x1[t],y1[t]}/.integrate],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotRange->{{mxlimit,xlimit},{mylimit,ylimit}},PlotStyle->{Blue,Thickness[0.003]},ColorFunction->Function[{x,y,w},Opacity[w,Blue]],PlotPoints->150,Frame->True,PerformanceGoal->"Quality"];
mass2=ParametricPlot[Evaluate[{x2[t],y2[t]}/.integrate],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotRange->{{mxlimit,xlimit},{mylimit,ylimit}},PlotStyle->{Black,Thickness[0.003]},ColorFunction->Function[{x,y,w},Opacity[w,Black]],PlotPoints->150,Frame->True,PerformanceGoal->"Quality"];
mass3=ParametricPlot[Evaluate[{x3[t],y3[t]}/.integrate],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotRange->{{mxlimit,xlimit},{mylimit,ylimit}},PlotStyle->{Red,Thickness[0.003]},ColorFunction->Function[{x,y,w},Opacity[w,Red]],PlotPoints->150,Frame->True,PerformanceGoal->"Quality"];
pointm1=Graphics[{PointSize[0.02 m1],Blue,Point[Dynamic[Evaluate[{x1[tmax],y1[tmax]}/.integrate]]]}];
pointm2=Graphics[{PointSize[0.02 m2],Black,Point[Dynamic[Evaluate[{x2[tmax],y2[tmax]}/.integrate]]]}];
pointm3=Graphics[{PointSize[0.02 m3],Red,Point[Dynamic[Evaluate[{x3[tmax],y3[tmax]}/.integrate]]]}];
Show[{mass1,mass2,mass3,pointm1,pointm2,pointm3},PlotLabel->Style[ "Position of masses",20],ImageSize->500,LabelStyle->18],
{tmax,0+$MachineEpsilon,Tmax},AnimationRunning->False,AnimationRate->0.3,SaveDefinitions->True,TrackedSymbols:>{tmax}]

Two sets of initial condition (demonstration of chaos)

Clear[t,tmax]
Manipulate[
integrate1=NDSolve[{
x1'[t]==vx1[t],y1'[t]==vy1[t],
x2'[t]==vx2[t],y2'[t]==vy2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m2 vx2'[t]==((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m2 vy2'[t]==((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vx3'[t]==((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vy3'[t]==((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),
x1[0]==x10,y1[0]==y10,
x2[0]==x20,y2[0]==y20,
x3[0]==x30,y3[0]==y30,
vx1[0]==vx10,vy1[0]==vy10,
vx2[0]==vx20,vy2[0]==vy20,
vx3[0]==vx30,vy3[0]==vy30},
{x1,x2,x3,y1,y2,y3},{t,0,tmax},MaxSteps->Infinity,PrecisionGoal->15];(*solved only for positions*)
integrate2=NDSolve[{
x1'[t]==vx1[t],y1'[t]==vy1[t],
x2'[t]==vx2[t],y2'[t]==vy2[t],
x3'[t]==vx3[t],y3'[t]==vy3[t],
m1 vx1'[t]==-((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m1 vy1'[t]==-((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2))),m2 vx2'[t]==((m1 m2(x1[t]-x2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m2 vy2'[t]==((m1 m2(y1[t]-y2[t]))/((((((x1[t]-x2[t]))^2)+(((y1[t]-y2[t]))^2)))^(3/2)))-((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vx3'[t]==((m1 m3(x1[t]-x3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(x2[t]-x3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),m3 vy3'[t]==((m1 m3(y1[t]-y3[t]))/((((((x1[t]-x3[t]))^2)+(((y1[t]-y3[t]))^2)))^(3/2)))+((m2 m3(y2[t]-y3[t]))/((((((x2[t]-x3[t]))^2)+(((y2[t]-y3[t]))^2)))^(3/2))),
x1[0]==x10+dx10,y1[0]==y10,
x2[0]==x20,y2[0]==y20,
x3[0]==x30,y3[0]==y30,
vx1[0]==vx10,vy1[0]==vy10,
vx2[0]==vx20,vy2[0]==vy20,
vx3[0]==vx30,vy3[0]==vy30},
{x1,x2,x3,y1,y2,y3},{t,0,tmax},MaxSteps->Infinity,PrecisionGoal->15];(*solved only for positions*)

initpos1=Graphics[{PointSize[0.02 m1],Blue,Point[{x10,y10}]}];
initpos2=Graphics[{PointSize[0.02 m2],Black,Point[{x20,y20}]}];
initpos3=Graphics[{PointSize[0.02 m3],Red,Point[{x30,y30}]}];
mass1=ParametricPlot[Evaluate[{x1[t],y1[t]}/.integrate1],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotLabel->Style["Mass 1",20],PlotRange->All,ImageSize->200,LabelStyle->18,Frame->True,AspectRatio->1,PerformanceGoal->"Quality"];
dmass1=ParametricPlot[Evaluate[{x1[t],y1[t]}/.integrate2],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotStyle->Dashed,PlotRange->All,ImageSize->200,LabelStyle->18,Frame->True,AspectRatio->1,PerformanceGoal->"Quality"];
mass2=ParametricPlot[Evaluate[{x2[t],y2[t]}/.integrate1],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotLabel->Style["Mass 2",20],PlotRange->All,PlotStyle->Black,ImageSize->200,LabelStyle->18,Frame->True,AspectRatio->1,PerformanceGoal->"Quality"];
dmass2=ParametricPlot[Evaluate[{x2[t],y2[t]}/.integrate2],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotRange->All,PlotStyle->{Dashed,Black},ImageSize->200,LabelStyle->18,Frame->True,AspectRatio->1,PerformanceGoal->"Quality"];
mass3=ParametricPlot[Evaluate[{x3[t],y3[t]}/.integrate1],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotLabel->Style["Mass 3",20],PlotRange->All,PlotStyle->Red,ImageSize->200,LabelStyle->18,Frame->True,AspectRatio->1,PerformanceGoal->"Quality"];
dmass3=ParametricPlot[Evaluate[{x3[t],y3[t]}/.integrate2],{t,0,tmax},FrameLabel->{x[t],y[t]},PlotRange->All,PlotStyle->{Dashed,Red},ImageSize->200,LabelStyle->18,Frame->True,AspectRatio->1,PerformanceGoal->"Quality"];
traj=Show[{mass1,dmass1,initpos1,mass2,dmass2,initpos2,mass3,dmass3,initpos3},PlotLabel-> Style["Positions of all masses together",20],ImageSize->500,LabelStyle->18,AspectRatio->1];
Grid[{{traj,SpanFromLeft},
{Show[mass1,dmass1,initpos1],Show[mass2,dmass2,initpos2],Show[mass3,dmass3,initpos3]}},Frame-> All],
{{tmax,0.25 Tmax,"Time"},0.001,Tmax,Appearance->"Labeled"},SaveDefinitions->True,TrackedSymbols:>{tmax}]

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