三体问题Mathematica数值解求解

本文参考百度经验如何用Mathematica计算三体问题？，进行了三体问题在三维空间的数值求解，给出了核心部分的计算代码。

注意我下面这些初值设定来自一篇2013年的一篇论文（感兴趣的可以私信我），在这种设定下能构成闭合的三体周期轨道。

以下这个部分是求解动力学方程的数值计算部分：

{{r1xo,r1yo,r1zo,v1xo,v1yo,v1zo,
r2xo,r2yo,r2zo,v2xo,v2yo,v2zo,
r3xo,r3yo,r3zo,v3xo,v3yo,v3zo}}
=NDSolve[{
(*第一组*)
r1x'[t]==v1x[t],r1y'[t]==v1y[t],r1z'[t]==v1z[t],
m1*v1x'[t]==-G*m2*m1*(r1x[t]-r2x[t])/lowR21[t]-G*m3*m1*(r1x[t]-r3x[t])/lowR31[t],
m1*v1y'[t]==-G*m2*m1*(r1y[t]-r2y[t])/lowR21[t]-G*m3*m1*(r1y[t]-r3y[t])/lowR31[t],
m1*v1z'[t]==-G*m2*m1*(r1z[t]-r2z[t])/lowR21[t]-G*m3*m1*(r1z[t]-r3z[t])/lowR31[t],
(*第二组*)
r2x'[t]==v2x[t],r2y'[t]==v2y[t],r2z'[t]==v2z[t],
m2*v2x'[t]==-G*m3*m2*(r2x[t]-r3x[t])/lowR32[t]-G*m1*m2*(r2x[t]-r1x[t])/lowR21[t],
m2*v2y'[t]==-G*m3*m2*(r2y[t]-r3y[t])/lowR32[t]-G*m1*m2*(r2y[t]-r1y[t])/lowR21[t],
m2*v2z'[t]==-G*m3*m2*(r2z[t]-r3z[t])/lowR32[t]-G*m1*m2*(r2z[t]-r1z[t])/lowR21[t],
(*第三组*)
r3x'[t]==v3x[t],r3y'[t]==v3y[t],r3z'[t]==v3z[t],
m3*v3x'[t]==-G*m3*m1*(r3x[t]-r1x[t])/lowR31[t]-G*m3*m2*(r3x[t]-r2x[t])/lowR32[t],
m3*v3y'[t]==-G*m3*m1*(r3y[t]-r1y[t])/lowR31[t]-G*m3*m2*(r3y[t]-r2y[t])/lowR32[t],
m3*v3z'[t]==-G*m3*m1*(r3z[t]-r1z[t])/lowR31[t]-G*m3*m2*(r3z[t]-r2z[t])/lowR32[t],
(*初值第1组*)
r1x[0]==Part[R10,1],r1y[0]==Part[R10,2],r1z[0]==Part[R10,3],v1x[0]==Part[V10,1],v1y[0]==Part[V10,2],v1z[0]==Part[V10,3],
(*初值第2组*)
r2x[0]==Part[R20,1],r2y[0]==Part[R20,2],r2z[0]==Part[R20,3],v2x[0]==Part[V20,1],v2y[0]==Part[V20,2],v2z[0]==Part[V20,3],
(*初值第3组*)
r3x[0]==Part[R30,1],r3y[0]==Part[R30,2],r3z[0]==Part[R30,3],
v3x[0]==Part[V30,1],v3y[0]==Part[V30,2],v3z[0]==Part[V30,3]},
{r1x,r1y,r1z,v1x,v1y,v1z,r2x,r2y,r2z,v2x,v2y,v2z,r3x,r3y,r3z,v3x,v3y,v3z},
{t,0,30}]

绘图

ParametricPlot3D[
{
{Evaluate[r1x[t] /. r1xo], Evaluate[r1y[t] /. r1yo],    Evaluate[r1z[t] /. r1zo]}, 
{Evaluate[r2x[t] /. r2xo],    Evaluate[r2y[t] /. r2yo],    Evaluate[r2z[t] /. r2zo]}, 
{Evaluate[r3x[t] /. r3xo],    Evaluate[r3y[t] /. r3yo], Evaluate[r3z[t] /. r3zo]}
}, 
{t, 0, 30}, 
PlotStyle -> {Red, Green, Blue}, PlotPoints -> 1000
]

下面是我最终导出的动图，如何导出动图仍然参考前面那篇百度经验.