Gravitatives Zwei- und Dreikörperproblem animiert

Anleitung: Zuerst die Zelle ausführen, dann die Animationen auswählen und mit  ▸  starten!
Die Begriffe “Sonne”, “Mond” und “Erde” sind nur der Veranschaulichung halber gewählt - die Animationen zeigen die gravitative Dynamik, die sich aus den Newtonschen Bewegungsgleichungen mit den jeweils gewählten Massen und Anfangsbedingungen ergeben.
Gravitationskonstante G=1.

(* Sonne und Erde *)

m1[zweik]=10;m2[zweik]=3;
y2ini[zweik]=-1;x2strichini[zweik]=(Sqrt[m1[zweik]/Abs[y2ini[zweik]]])-0.4;
r12[zweik]=(Sqrt[(((x1[zweik][t]-x2[zweik][t]))^2)+(((y1[zweik][t]-y2[zweik][t]))^2)]);
(* Print["Lösung des Differentialgleichungs-Systems:"]; *)
solution2K=NDSolve[{
x1[zweik]''[t]==m2[zweik]((x2[zweik][t]-x1[zweik][t])/((r12[zweik])^3)),
y1[zweik]''[t]==m2[zweik]((y2[zweik][t]-y1[zweik][t])/((r12[zweik])^3)),
x2[zweik]''[t]==m1[zweik]((x1[zweik][t]-x2[zweik][t])/((r12[zweik])^3)),
y2[zweik]''[t]==m1[zweik]((y1[zweik][t]-y2[zweik][t])/((r12[zweik])^3)),
x1[zweik][0]==0,y1[zweik][0]==0,x1[zweik]'[0]==0,y1[zweik]'[0]==0,
x2[zweik][0]==0,y2[zweik][0]==y2ini[zweik],x2[zweik]'[0]==x2strichini[zweik],y2[zweik]'[0]==0
},{x1[zweik][t],y1[zweik][t],x2[zweik][t],y2[zweik][t]},{t,0,10}];
X1[zweik][t_]=x1[zweik][t]/.solution2K[[1]];
X2[zweik][t_]=x2[zweik][t]/.solution2K[[1]];
Y1[zweik][t_]=y1[zweik][t]/.solution2K[[1]];
Y2[zweik][t_]=y2[zweik][t]/.solution2K[[1]];
mmX[zweik][t_]=((m1[zweik] X1[zweik][t]+m2[zweik] X2[zweik][t])/(m1[zweik]+m2[zweik]));mmY[zweik][t_]=((m1[zweik] Y1[zweik][t]+m2[zweik] Y2[zweik][t])/(m1[zweik]+m2[zweik]));
run2K:=Manipulate[Graphics[{
{Orange,Disk[{X1[zweik][t],Y1[zweik][t]},.06]},
{Blue,Disk[{X2[zweik][t],Y2[zweik][t]},.03]},If[massenmittelpunkt,{Red,{Line[{{mmX[zweik][t]-0.04,mmY[zweik][t]},{mmX[zweik][t]+0.04,mmY[zweik][t]}}],Line[{{mmX[zweik][t],mmY[zweik][t]-0.04},{mmX[zweik][t],mmY[zweik][t]+0.04}}]
}},{}]
},PlotRange->{{-0.85,6.75},{-1.29,1.15}},Axes->achsen,ImageSize->700],{{achsen,False,"Achsen"},{True,False}},{{massenmittelpunkt,False,"Massenmittelpunkt"},{True,False}},{t,0,10,0.01,AnimationRate->0.5,Appearance->"Open"},ControlPlacement->{Bottom,Bottom,Top}];

(* ------------------------------------------------------------------------------------------ *)
(* Sonne, Mond und Erde *)

m1[dreik1]=10;m2[dreik1]=3;m3[dreik1]=1;
y2ini[dreik1]=-1;x2strichini[dreik1]=(Sqrt[m1[dreik1]/Abs[y2ini[dreik1]]])-0.4;
y3ini[dreik1]=-1.2;x3strichini[dreik1]=(Sqrt[m2[dreik1]/Abs[y3ini[dreik1]-y2ini[dreik1]]])+3;
r12[dreik1]=(Sqrt[(((x1[dreik1][t]-x2[dreik1][t]))^2)+(((y1[dreik1][t]-y2[dreik1][t]))^2)]);
r13[dreik1]=(Sqrt[(((x1[dreik1][t]-x3[dreik1][t]))^2)+(((y1[dreik1][t]-y3[dreik1][t]))^2)]);
r23[dreik1]=(Sqrt[(((x2[dreik1][t]-x3[dreik1][t]))^2)+(((y2[dreik1][t]-y3[dreik1][t]))^2)]);
(* Print["Lösung des Differentialgleichungs-Systems:"]; *)
solution3K[dreik1]=NDSolve[{
x1[dreik1]''[t]==m2[dreik1]((x2[dreik1][t]-x1[dreik1][t])/((r12[dreik1])^3))+m3[dreik1]((x3[dreik1][t]-x1[dreik1][t])/((r13[dreik1])^3)),
y1[dreik1]''[t]==m2[dreik1]((y2[dreik1][t]-y1[dreik1][t])/((r12[dreik1])^3))+m3[dreik1]((y3[dreik1][t]-y1[dreik1][t])/((r13[dreik1])^3)),
x2[dreik1]''[t]==m1[dreik1]((x1[dreik1][t]-x2[dreik1][t])/((r12[dreik1])^3))+m3[dreik1]((x3[dreik1][t]-x2[dreik1][t])/((r23[dreik1])^3)),
y2[dreik1]''[t]==m1[dreik1]((y1[dreik1][t]-y2[dreik1][t])/((r12[dreik1])^3))+m3[dreik1]((y3[dreik1][t]-y2[dreik1][t])/((r23[dreik1])^3)),
x3[dreik1]''[t]==m1[dreik1]((x1[dreik1][t]-x3[dreik1][t])/((r13[dreik1])^3))+m2[dreik1]((x2[dreik1][t]-x3[dreik1][t])/((r23[dreik1])^3)),
y3[dreik1]''[t]==m1[dreik1]((y1[dreik1][t]-y3[dreik1][t])/((r13[dreik1])^3))+m2[dreik1]((y2[dreik1][t]-y3[dreik1][t])/((r23[dreik1])^3)),

x1[dreik1][0]==0,y1[dreik1][0]==0,x1[dreik1]'[0]==0,y1[dreik1]'[0]==0,
x2[dreik1][0]==0,y2[dreik1][0]==y2ini[dreik1],x2[dreik1]'[0]==x2strichini[dreik1],y2[dreik1]'[0]==0,
x3[dreik1][0]==0,y3[dreik1][0]==y3ini[dreik1],x3[dreik1]'[0]==x3strichini[dreik1],y3[dreik1]'[0]==0
},{x1[dreik1][t],y1[dreik1][t],x2[dreik1][t],y2[dreik1][t],x3[dreik1][t],y3[dreik1][t]},{t,0,10},MaxSteps->50000,MaxStepFraction->1/100,PrecisionGoal->10];
X1[dreik1][t_]=x1[dreik1][t]/.solution3K[dreik1][[1]];
X2[dreik1][t_]=x2[dreik1][t]/.solution3K[dreik1][[1]];
X3[dreik1][t_]=x3[dreik1][t]/.solution3K[dreik1][[1]];
Y1[dreik1][t_]=y1[dreik1][t]/.solution3K[dreik1][[1]];
Y2[dreik1][t_]=y2[dreik1][t]/.solution3K[dreik1][[1]];
Y3[dreik1][t_]=y3[dreik1][t]/.solution3K[dreik1][[1]];
mmX[dreik1][t_]=((m1[dreik1] X1[dreik1][t]+m2[dreik1] X2[dreik1][t]+m3[dreik1] X3[dreik1][t])/(m1[dreik1]+m2[dreik1]+m3[dreik1]));mmY[dreik1][t_]=((m1[dreik1] Y1[dreik1][t]+m2[dreik1] Y2[dreik1][t]+m3[dreik1] Y3[dreik1][t])/(m1[dreik1]+m2[dreik1]+m3[dreik1]));
run3K[dreik1]:=Manipulate[Graphics[{
{Orange,Disk[{X1[dreik1][t],Y1[dreik1][t]},.1]},
{Blue,Disk[{X2[dreik1][t],Y2[dreik1][t]},.04]},
{RGBColor[0,0.9,0],Disk[{X3[dreik1][t],Y3[dreik1][t]},.03]},
If[massenmittelpunkt,{Red,{Line[{{mmX[dreik1][t]-0.07,mmY[dreik1][t]},{mmX[dreik1][t]+0.07,mmY[dreik1][t]}}],Line[{{mmX[dreik1][t],mmY[dreik1][t]-0.07},{mmX[dreik1][t],mmY[dreik1][t]+0.07}}]
}},{}]
},PlotRange->{{-0.85,11.3},{-2.29,2.15}},Axes->achsen,ImageSize->700],{{achsen,False,"Achsen"},{True,False}},{{massenmittelpunkt,False,"Massenmittelpunkt"},{True,False}},{t,0,10,0.001,AnimationRate->0.25,Appearance->"Open"},ControlPlacement->{Bottom,Bottom,Top}];

(* ------------------------------------------------------------------------------------------ *)
(* Sonne und Erde stehlen einander den Mond *)

m1[dreik2]=10;m2[dreik2]=3;m3[dreik2]=1;
y2ini[dreik2]=-1;x2strichini[dreik2]=(Sqrt[m1[dreik2]/Abs[y2ini[dreik2]]])-0.4;
y3ini[dreik2]=-1.29;x3strichini[dreik2]=(Sqrt[m2[dreik2]/Abs[y3ini[dreik2]-y2ini[dreik2]]])+3;
r12[dreik2]=(Sqrt[(((x1[dreik2][t]-x2[dreik2][t]))^2)+(((y1[dreik2][t]-y2[dreik2][t]))^2)]);
r13[dreik2]=(Sqrt[(((x1[dreik2][t]-x3[dreik2][t]))^2)+(((y1[dreik2][t]-y3[dreik2][t]))^2)]);
r23[dreik2]=(Sqrt[(((x2[dreik2][t]-x3[dreik2][t]))^2)+(((y2[dreik2][t]-y3[dreik2][t]))^2)]);
(* Print["Lösung des Differentialgleichungs-Systems:"]; *)
solution3K[dreik2]=NDSolve[{
x1[dreik2]''[t]==m2[dreik2]((x2[dreik2][t]-x1[dreik2][t])/((r12[dreik2])^3))+m3[dreik2]((x3[dreik2][t]-x1[dreik2][t])/((r13[dreik2])^3)),
y1[dreik2]''[t]==m2[dreik2]((y2[dreik2][t]-y1[dreik2][t])/((r12[dreik2])^3))+m3[dreik2]((y3[dreik2][t]-y1[dreik2][t])/((r13[dreik2])^3)),
x2[dreik2]''[t]==m1[dreik2]((x1[dreik2][t]-x2[dreik2][t])/((r12[dreik2])^3))+m3[dreik2]((x3[dreik2][t]-x2[dreik2][t])/((r23[dreik2])^3)),
y2[dreik2]''[t]==m1[dreik2]((y1[dreik2][t]-y2[dreik2][t])/((r12[dreik2])^3))+m3[dreik2]((y3[dreik2][t]-y2[dreik2][t])/((r23[dreik2])^3)),
x3[dreik2]''[t]==m1[dreik2]((x1[dreik2][t]-x3[dreik2][t])/((r13[dreik2])^3))+m2[dreik2]((x2[dreik2][t]-x3[dreik2][t])/((r23[dreik2])^3)),
y3[dreik2]''[t]==m1[dreik2]((y1[dreik2][t]-y3[dreik2][t])/((r13[dreik2])^3))+m2[dreik2]((y2[dreik2][t]-y3[dreik2][t])/((r23[dreik2])^3)),

x1[dreik2][0]==0,y1[dreik2][0]==0,x1[dreik2]'[0]==0,y1[dreik2]'[0]==0,
x2[dreik2][0]==0,y2[dreik2][0]==y2ini[dreik2],x2[dreik2]'[0]==x2strichini[dreik2],y2[dreik2]'[0]==0,
x3[dreik2][0]==0,y3[dreik2][0]==y3ini[dreik2],x3[dreik2]'[0]==x3strichini[dreik2],y3[dreik2]'[0]==0
},{x1[dreik2][t],y1[dreik2][t],x2[dreik2][t],y2[dreik2][t],x3[dreik2][t],y3[dreik2][t]},{t,0,15},MaxSteps->50000,MaxStepFraction->1/100,PrecisionGoal->10];
X1[dreik2][t_]=x1[dreik2][t]/.solution3K[dreik2][[1]];
X2[dreik2][t_]=x2[dreik2][t]/.solution3K[dreik2][[1]];
X3[dreik2][t_]=x3[dreik2][t]/.solution3K[dreik2][[1]];
Y1[dreik2][t_]=y1[dreik2][t]/.solution3K[dreik2][[1]];
Y2[dreik2][t_]=y2[dreik2][t]/.solution3K[dreik2][[1]];
Y3[dreik2][t_]=y3[dreik2][t]/.solution3K[dreik2][[1]];
mmX[dreik2][t_]=((m1[dreik2] X1[dreik2][t]+m2[dreik2] X2[dreik2][t]+m3[dreik2] X3[dreik2][t])/(m1[dreik2]+m2[dreik2]+m3[dreik2]));mmY[dreik2][t_]=((m1[dreik2] Y1[dreik2][t]+m2[dreik2] Y2[dreik2][t]+m3[dreik2] Y3[dreik2][t])/(m1[dreik2]+m2[dreik2]+m3[dreik2]));
run3K[dreik2]:=Manipulate[Graphics[{
{Orange,Disk[{X1[dreik2][t],Y1[dreik2][t]},.15]},
{Blue,Disk[{X2[dreik2][t],Y2[dreik2][t]},.06]},
{RGBColor[0,0.9,0],Disk[{X3[dreik2][t],Y3[dreik2][t]},.05]},
If[massenmittelpunkt,{Red,{Line[{{mmX[dreik2][t]-0.09,mmY[dreik2][t]},{mmX[dreik2][t]+0.09,mmY[dreik2][t]}}],Line[{{mmX[dreik2][t],mmY[dreik2][t]-0.09},{mmX[dreik2][t],mmY[dreik2][t]+0.09}}]
}},{}]
},PlotRange->{{-0.85,17.3},{-3.49,3.15}},Axes->achsen,ImageSize->700],{{achsen,False,"Achsen"},{True,False}},{{massenmittelpunkt,False,"Massenmittelpunkt"},{True,False}},{t,0,15,0.001,AnimationRate->0.25,Appearance->"Open"},ControlPlacement->{Bottom,Bottom,Top}];

(* ------------------------------------------------------------------------------------------ *)
(* Der Mond verabschiedet sich *)

m1[dreik3]=10;m2[dreik3]=3;m3[dreik3]=1;
y2ini[dreik3]=-1;x2strichini[dreik3]=(Sqrt[m1[dreik3]/Abs[y2ini[dreik3]]])-0.4;
y3ini[dreik3]=-1.5;x3strichini[dreik3]=(Sqrt[m2[dreik3]/Abs[y3ini[dreik3]-y2ini[dreik3]]])+2;
r12[dreik3]=(Sqrt[(((x1[dreik3][t]-x2[dreik3][t]))^2)+(((y1[dreik3][t]-y2[dreik3][t]))^2)]);
r13[dreik3]=(Sqrt[(((x1[dreik3][t]-x3[dreik3][t]))^2)+(((y1[dreik3][t]-y3[dreik3][t]))^2)]);
r23[dreik3]=(Sqrt[(((x2[dreik3][t]-x3[dreik3][t]))^2)+(((y2[dreik3][t]-y3[dreik3][t]))^2)]);
(* Print["Lösung des Differentialgleichungs-Systems:"]; *)
solution3K[dreik3]=NDSolve[{
x1[dreik3]''[t]==m2[dreik3]((x2[dreik3][t]-x1[dreik3][t])/((r12[dreik3])^3))+m3[dreik3]((x3[dreik3][t]-x1[dreik3][t])/((r13[dreik3])^3)),
y1[dreik3]''[t]==m2[dreik3]((y2[dreik3][t]-y1[dreik3][t])/((r12[dreik3])^3))+m3[dreik3]((y3[dreik3][t]-y1[dreik3][t])/((r13[dreik3])^3)),
x2[dreik3]''[t]==m1[dreik3]((x1[dreik3][t]-x2[dreik3][t])/((r12[dreik3])^3))+m3[dreik3]((x3[dreik3][t]-x2[dreik3][t])/((r23[dreik3])^3)),
y2[dreik3]''[t]==m1[dreik3]((y1[dreik3][t]-y2[dreik3][t])/((r12[dreik3])^3))+m3[dreik3]((y3[dreik3][t]-y2[dreik3][t])/((r23[dreik3])^3)),
x3[dreik3]''[t]==m1[dreik3]((x1[dreik3][t]-x3[dreik3][t])/((r13[dreik3])^3))+m2[dreik3]((x2[dreik3][t]-x3[dreik3][t])/((r23[dreik3])^3)),
y3[dreik3]''[t]==m1[dreik3]((y1[dreik3][t]-y3[dreik3][t])/((r13[dreik3])^3))+m2[dreik3]((y2[dreik3][t]-y3[dreik3][t])/((r23[dreik3])^3)),

x1[dreik3][0]==0,y1[dreik3][0]==0,x1[dreik3]'[0]==0,y1[dreik3]'[0]==0,
x2[dreik3][0]==0,y2[dreik3][0]==y2ini[dreik3],x2[dreik3]'[0]==x2strichini[dreik3],y2[dreik3]'[0]==0,
x3[dreik3][0]==0,y3[dreik3][0]==y3ini[dreik3],x3[dreik3]'[0]==x3strichini[dreik3],y3[dreik3]'[0]==0
},{x1[dreik3][t],y1[dreik3][t],x2[dreik3][t],y2[dreik3][t],x3[dreik3][t],y3[dreik3][t]},{t,0,5},MaxSteps->50000,MaxStepFraction->1/100,PrecisionGoal->10];
X1[dreik3][t_]=x1[dreik3][t]/.solution3K[dreik3][[1]];
X2[dreik3][t_]=x2[dreik3][t]/.solution3K[dreik3][[1]];
X3[dreik3][t_]=x3[dreik3][t]/.solution3K[dreik3][[1]];
Y1[dreik3][t_]=y1[dreik3][t]/.solution3K[dreik3][[1]];
Y2[dreik3][t_]=y2[dreik3][t]/.solution3K[dreik3][[1]];
Y3[dreik3][t_]=y3[dreik3][t]/.solution3K[dreik3][[1]];
mmX[dreik3][t_]=((m1[dreik3] X1[dreik3][t]+m2[dreik3] X2[dreik3][t]+m3[dreik3] X3[dreik3][t])/(m1[dreik3]+m2[dreik3]+m3[dreik3]));mmY[dreik3][t_]=((m1[dreik3] Y1[dreik3][t]+m2[dreik3] Y2[dreik3][t]+m3[dreik3] Y3[dreik3][t])/(m1[dreik3]+m2[dreik3]+m3[dreik3]));
run3K[dreik3]:=Manipulate[Graphics[{
{Orange,Disk[{X1[dreik3][t],Y1[dreik3][t]},.07]},
{Blue,Disk[{X2[dreik3][t],Y2[dreik3][t]},.025]},
{RGBColor[0,0.9,0],Disk[{X3[dreik3][t],Y3[dreik3][t]},.02]},
If[massenmittelpunkt,{Red,{Line[{{mmX[dreik3][t]-0.05,mmY[dreik3][t]},{mmX[dreik3][t]+0.05,mmY[dreik3][t]}}],Line[{{mmX[dreik3][t],mmY[dreik3][t]-0.05},{mmX[dreik3][t],mmY[dreik3][t]+0.05}}]
}},{}]
},PlotRange->{{-0.95,5.3},{-1.59,0.81}},Axes->achsen,ImageSize->700],{{achsen,False,"Achsen"},{True,False}},{{massenmittelpunkt,False,"Massenmittelpunkt"},{True,False}},{t,0,5,0.001,AnimationRate->0.25,Appearance->"Open"},ControlPlacement->{Bottom,Bottom,Top}];

(* ------------------------------------------------------------------------------------------ *)

TabView[{
"Sonne und Erde"->Labeled[run2K,Row[{"\nMassen: ",m1[zweik], " und ",m2[zweik]}]],
"Sonne, Mond und Erde"->Labeled[run3K[dreik1],Row[{"Massen: ",m1[dreik1], ", ", m2[dreik1]," und ",m3[dreik1]}]],
"Sonne und Erde stehlen einander den Mond"->Labeled[run3K[dreik2],Row[{"Massen: ",m1[dreik2], ", ", m2[dreik2]," und ",m3[dreik2]}]],
"Der Mond verabschiedet sich"->Labeled[run3K[dreik3],Row[{"Massen: ",m1[dreik3], ", ", m2[dreik3]," und ",m3[dreik3]}]]}]