天遇---地球在宇宙流浪，人类在地球流浪---天体力学的几个图像切片

查看原文

其他

天遇---地球在宇宙流浪，人类在地球流浪---天体力学的几个图像切片

Original 洞穴之外理念世界的影子 2021-06-23

公众号：理念世界的影子

文不可无观点，观点不可无论据。

转载请注明出处

太阳系是稳定的吗？确切地讲，答案还不知道，但是由这个问题已经引出了许多深刻的结果，它们可能比对初始问题的解答更为重要。---莫泽

流浪地球票房大卖，本号也借个势。按著名的六神磊磊的话说，咱这是姨妈精(《杠精、屁精和姨妈精》)。本篇是上个姨妈精《永元搞定范冰冰---FBB、Matlab与航天》的姐妹篇。

引子

只是这里没有氦闪，没有红巨星，没有行星发动机。但是这里有三体、有引力弹弓、有洛希极限，这里是天体力学。

有个问题是：假设太阳系几大行星永不消亡，仅考虑它们间万有引力作用下的轨道运动，地球会不会永远稳定地围绕太阳旋转？

康德说：有两种东西，我对它们的思考越是深沉和持久，它们在我心灵中唤起的惊奇和敬畏就会日新月异，不断增长，这就是我头上的星空和心中的道德定律。

冯友兰说：在道德伦理之外还有更高的价值，可以称之为超道德价值。人不满足于现实世界而追求超越现实世界，这是人类内心深处的一种渴望。譬如为什么有星星？星星为什么会动，它是怎么动的？。。。这是孩童常挂在嘴边的问题，我们追求和欣赏星空的运行规律和永恒的世界。

这个世界太过神奇，为满足超越现实世界的追求，西方曾求助于宗教，但它的权威在科学前进的历程中不断被削弱。无奈对待太阳系，科学也还是不给力，最后一位科学通才庞加莱甩出了一句听起来特别牛逼的话：

繁星无法超越

不要被这句话的气势震慑，这句话其实没有太多微言大义，它的含义一如字面，阐述的是一个事实：太阳系轨道运动会不会稳定这个问题，天纵英才的庞加莱穷尽一生，回答不了。

地球的流浪---庞加莱初窥混沌

三体人在是否能见到明天的太阳的恐惧中，度过了一次又一次的毁灭与重生。明天的太阳会不会照常升起，这是个古老的问题。

如果太阳系中的星体不发生碰撞并且其中的行星不能从太阳系中逃逸出去，那么太阳系就被称作稳定的。

在牛顿之前，我们不可能对太阳系的稳定性进行严格的研究。

经过第谷艰苦卓绝的观察和记录，以及开普勒从大量数据中反演出的三大定律“为天空立法”，人们已知知道地球绕太阳呈椭圆旋转，但在牛顿之前，没有人通过数学原理进行精确论证。

牛顿，1687，《自然哲学的数学原理》，第三编《宇宙体系》

命题7：一切物体都有受一种引力的吸引，该引力正比于物体各自所含的质量。

命题8：两个互相吸引的球体，如果球体内到球心距离相等处的物质是相似的，则其中一个球体的重量与另一个的重量之比反比于它们球心之间的距离的平方。

当存在两个天体时，利用万有引力定律可以解出两个天体的运动。

图二体运动

看起来很复杂，换个角度，站在第一个天体上看相对位置，则另一个质点的轨道一定是个椭圆，抛物线，双曲线的一支或者直线(数学上统称为圆锥曲线)。

图二体运动相对位置

将两个天体扩展到多个是个很自然的问题，在数学上，这个问题叫：N体问题。其中最简单的是三体，就是刘慈欣小说下的《三体》。

但牛顿、欧拉和其他人在试图理解更多天体时碰到了困难。

在没有计算机前，一般采用摄动法（级数展开，忽略次要部分计算）计算近似解。如1846年，法国教师勒威耶根据天王星理论轨道与真实轨道差别，计算出影响天王星运动的第八颗行星轨道，9月23日，德国天文学家在这个位置发现了海王星。下图为摄动的示意图，分别为绝对空间和相对空间下的图像，轨道已经不再那么规律。

图摄动轨道示意图(从恒星看)

图摄动轨道示意图(从第一颗行星看)

有没有看到，被摄动的轨道有内有外。在整个16和17世纪，天文学家发现，木星正在缓慢地螺旋向内运动，而土星正在逐渐向外，如果这一趋势继续进行几万年，太阳系就会陷入危机。虽然付出了巨大努力，但牛顿无法用他的理论解释正在发生的现象。他承认：如果我没有弄错的话，同时考虑所有对（行星）运动的影响并且使用简单计算即可得的精确定律来定义这些运动超出了人类思维所能及的范围。

还好，拉普拉斯解决了这个问题，他证明了这两个轨道存在一个周期为数千年的振荡，木星和土星的轨道在千年后还会回来，而不会一边倒。

拉普拉斯第一个把牛顿的万有引力定律应用到了整个太阳系，太阳系稳定性问题的研究始于拉普拉斯，他得到了一个乐观结果，但远非最终的结果（他的结果是如果考虑离心率的一阶近似，行星半长轴不含久期项，即随时间累积不发散）。但这不妨碍拉普拉斯的乐观主义，在拿破仑大帝问拉普拉斯，他的著作《天体力学》中为什么没有上帝，拉普拉斯回答：陛下，我不需要上帝这个假设。

好景不长，继拉普拉斯、拉格朗日、泊松给出太阳系稳定性的较为乐观的结果后，1877年哈雷特的结果则不那么好了，他证明了行星轨道长轴的值关于质量的三阶幂级数逼近中确实存在久期项，即行星轨道的长轴在时间变化时不一定是稳定的(但不能坐实，如果久期项能以某种巧妙方式抵消，则仍然可能是稳定的)。

鉴于太阳系稳定性、n体问题解析解问题迟迟得不到解决，同时又这么充满了超道德价值，一时成为巴黎时髦沙龙的话题，也曾引起伏尔泰和曼坦侬夫人等非科学人士的兴趣。1885年瑞典与挪威国王奥斯卡二世被说服设立了一个奖项，以庆祝1889年1月21日国王的六十寿辰。推荐的理由是促进天体力学与数学的研究，提高国王的威望，提高瑞典数学界以及他的杂志的影响力。国王欣然接受这个建议，事实上，国王是明智的，后来人知道他，大半因为他设立的这个奖。

关于这个奖的竞赛公告第一段是要解决这个问题：“具有任意多个质点的系统，其中任何两点间的作用力满足牛顿定律，在任意两个点不发生碰撞的情况下，试给出每个点的坐标，这个坐标可以以时间的某个已知函数作为变量的级数表示，并且对于所有的取值，该级数是一致收敛的。”

庞加莱出手了，他首先证明了，试图写出三体问题一般化的显式表达式和积分的定量方法是不可能的。为解决这个问题，他创立了微分方程定性理论。定量方法只能解决某些问题，人们还应当定性地和几何地思考。

图庞加莱

他的理论和证明最后写成了《天体力学新方法》这本书，在最后一卷第397节开头时，庞加莱用了这样一段描述：...双向渐进解，这两条曲线中的任何一条都不会再次与自身相交，但是它以一种极其复杂的方式弯曲着回到自身附近，从而无限次地穿过格子结构中的所有网眼...

庞加莱对此发现感到震惊，他提出，这样的图形因过于复杂而很难画出来。事实上，这些相交曲线会形成网状构图，并不断穿过网眼，形成更细密更复杂纷乱的网眼，因此系统将陷入混乱无序之中。这就是后来科学家所定义的“混沌”（chaos）的雏形。庞加莱因此成为窥见确定性系统产生不确定性的第一人，他对自己的发现非常震惊，并很难接受，因为这与他所处时代的哲学是背道而驰的。

根据参赛规则，只能通过一段题记来确认身份，庞加莱选择了“繁星无法超越”。

庞加莱说：“这一问题将为未来几代的数学家提供取之不竭的成果来源，如此多的工作需要去完成，即使是向前迈几步，也会占用我的余生。现在至少我懂得，没有人能够独立地解决这个问题。我们这代人不可能完全理解它。”

今天绝大多数人应该理解不了庞加莱的证明，实际上，即使是受过不少科学训练的笔者，甚至连《天遇》这本有趣的科普书都看不懂。但今天，我们有了计算机，那些复杂的结论，我们反而都能接受了。

图 "混沌"轨道

注:由于本文只是简单地采用Runge-Kutta数值积分算法进行轨道计算，此方法存在数值耗散，无法进行长时间轨道动力学模拟，因此无法给出图形，要进行长时间轨道动力学模拟，需要使用辛积分。

1989年，法国精度管理局的Jacques Laskar发表了他对太阳系超过2亿年的数值积分结果，他沿线运用了拉普拉斯平衡方程而非完整的动量方程。Lasker的工作显示，地球轨道是混沌的，现在对地球轨道15米那么小的测量误差就会导致仅仅1亿年后其轨道的不可测量。(更多信息见：太阳系是稳定的吗？http://bolide.lamost.org/articles/article245.htm)

没有尽头的流浪---夏志宏的超级引力弹弓

既然三体问题没有一般解，那有没有几个特殊的构型，能给出已知解呢？数学家们对此事很感兴趣。

1767年，欧拉给出，三个天体共线，边上两颗围绕中间一个旋转，由于对称性，这个解是显而易见的。

图 L2点

1772年，拉格朗日给出，三星组成三角形，围绕三角形中心旋转。

图 L4点

他们给出的就是五个拉格朗日点，这里用到了一个蜕化的假设，即第三个天体的质量为无穷小，它仅受其他两个天体万有引力作用，而自身对其他天体无作用。这个称为限制性三体问题。

拉格朗日点现在已取得广泛的应用，如嫦娥四号的中继星鹊桥就运行在L2点附近的Halo轨道上。自然天体中也有在拉格朗日点的，详细可阅读“解读《三体》7. 三体运动：欧拉、拉格朗日特解与三体运动的新特解 ”，本文不再抄一遍了。

然后就是多年没有进展，直到1970年，Michel Hénon和Roger A. Broucke利用计算机进行模拟，发现了更多结果，并将之归结为3族：

欧拉.拉格朗日族；
Broucke-Hadjidemetriou-Hénon 族；
8字形族。

下图为一个8字形轨道。

图 8字形轨道

2013年，研究获得了一个突破性进展，Milovan Šuvakov and Veljko Dmitrašinović在《物理评论快报》上发表文章，宣布发现13族新的特解。

2017年，上海交通大学的廖世俊教授利用计算机，又发现了669种，相关初始参数可在http://numericaltank.sjtu.edu.cn/three-body/three-body-movies.htm网站查到，并进行模拟，看看那些奇迹般的形状。

图复杂轨道1---蝴蝶型

图复杂轨道2

在寻找稳定解的同时，也开展了寻找不稳定解的探索。

从牛顿运动定律看，万有引力与距离平方反比成正比，数学上两天体相遇时，万有引力为无穷大（物理上是撞上），这在数学上称为奇点。数学家们关心，对于N体问题，有没有可能存在其他非碰撞的奇点。只能说，数学家们都对这些纯粹的事物极其感兴趣。

1895年，法国数学家潘勒韦提出了这个问题，他还证明了对于三体问题，存在的奇点都是碰撞，但对于n大于3的情况是不是这样？谁都不知道？潘勒韦猜想这个结果是存在的（想起费马大定理了吗？这些猜想是不是都有异曲同工之妙？），他心目中的图像是：质点剧烈地振动，接近碰撞，接着在另一次碰撞更为接近之前分散开来，上述过程将在有限时间内不断重复。

1966年，塞拜海伊做了一个数值试验，考虑质量比为3：4：5的三个质点，分别位于直角三角形的顶点，假设这些质量从零速释放，结果会怎么样？

图三体问题中的毕达哥拉斯问题

它们会直接全部撞上吗？未必？因为两两对撞上，会碰到第三个质点的引力而擦肩而过。这个现象只能用计算机进行计算。

图毕达哥拉斯三角形弹弓

结果令人吃惊！在经过一段时间复杂运动过程后，其中两个质点形成了二体系统紧密地互相绕转，而第三个质点以很高的速度被排斥出去！

之所以是这样，是因为中间有个时刻(t=41附近)，三个质点非常靠近，按照牛顿万有引力定律，两个质点间的作用力由于它们相互靠近而变得非常大，这就是"引力弹弓"效应(slingshot effect)。

引力弹弓其实就是一个二体问题，从中心天体看，结果就是一个圆锥曲线轨道，这就造成了两个结果，一是在靠近焦点时速度加快，二是经过中心天体后速度转向，此时被弹弓甩摆出去的天体绝对速度就叠加上了中心天体的速度，可以起到加速和减速作用。

图引力弹弓示意图(红色为此天体速度值曲线)

只是计算机模拟不能代替数学证明，计算机有助于直观地认识，但真正的理解还得依赖分析。中间经过冯.塞佩尔、麦吉、马瑟、杰弗、萨瑞的铺垫性工作后，1987年，中国学者夏志宏终于证明了潘勒韦猜想，在这座建造已逾90年的金字塔的顶端，放上了最后一块石头。

图夏志宏(网上有现在的清晰照，但笔者喜欢书上这张年轻时的帅帅照)

这次，夏志宏构造的是空间五体问题，尽管增加了变量数目，但他通过设定某些对称性简化了分析。

图夏志宏构造的非碰撞奇点的例子

在夏志宏构造的方案中，m5在穿过m3、m4时，三者正好接近于三体碰撞，从而被加速，它上升了一段距离后受到强烈吸引向下运行，之后m3、m4相互分离，m5快速向下经过m1和m2，再重复这个过程。上述步骤往复循环，m5每次都更加接近二体系统，并且获得更大的加速度，而且运动在有限时间内达到无穷远。这就是潘勒韦猜测存在的奇点。笔者没有找论文或其他介绍，不知道相应设置的初值，这里就不放仿真图了。

这里的方案定性听起来很简单，但为什么每组二体问题的质点会越来越靠近，而且不会碰撞，需要严格证明。夏志宏在每一步都巧妙地运用碰撞点的稳定流形和不稳定流形结构排除“坏的”解的集合，即排除在下次相遇时将导致碰撞的解以及m5在错误的时间接近二体系统的解。除了正确“瞄准”m5，二体系统的角速度需要在每一步都减少一点，从而与m5相遇时会更加接近。所有过程都需要非常精确的计算。在无穷多次重复相遇过程之后，留下了由解的初值条件构成的康托尔集，它们将导致非碰撞奇点。证明不同于计算机仿真，在证明中夏志宏克服了大量的困难，运用了许多数学技巧，引入了许多新的思想。总之，夏氏定理(Xia’s theorem)位于了这座大厦的顶端。

这就是夏志宏和他的超级引力弹弓，一个瞬间将天体加速到无穷远处的引力弹弓。谁说中国人性世俗，讲究实用？在超道德关怀上，我们也曾经走在了世界的巅峰。

流浪的代价---洛希极限

在引力弹弓中，特别是超级引力弹弓，需要天体间的距离越近越好，但由于天体存在半径，在距离小于半径之和时，两个天体会产生碰撞。

而且在碰撞之前，由于天体A对天体B两侧的吸引力不同，当吸引力大于天体B自身的引力时，天体B会被撕碎，这个撕碎的距离就被称为洛希极限。

不多说，上图，由于matlab程序计算力有限，粒子少了点，看起来有点丑，因此还从网上扒拉了一个漂亮一点的效果图。

图洛希极限示意图

图洛希极限效果图

人类的流浪---确定性的终结

19世纪末人类的精神是积极进取和乐观的，人们相信可以用人类的聪明才智解决任何问题。在拉普拉斯的数学中，如果知道宇宙中每一个粒子的形状，就能彻底地预言它们的整个未来，不需要上帝。

经典科学强调有序和稳定性，20世纪物理学诞生了相对论、量子力学和混沌。现在，我们在观测的所有层次上都看到了涨落、不稳定性、多种选择和有限可预测性，像混沌这样的思想已变得相当流行，影响着从宇宙学到经济学，实际上所有科学领域乃至社会领域的思想。现代科学里，时间是一种错觉、上帝在掷骰子、繁星也无法超越。

身处17世纪的牛顿，面对如此广阔的宇宙，他真的感到了造物的伟大，自我的渺小。“我不知道我可以向世界呈现什么，但是对于我自己来说，我似乎只是像一个在海岸边玩耍的孩子，不时为发现比寻常更为美丽的一块卵石或一片贝壳而沾沾自喜，至于展现在我面前的浩瀚的真理海洋，却全然没有发现。”

参考文献

弗洛林.迪亚库、菲利普.霍尔姆斯，王兰宇译，天遇---混沌与稳定性的起源，上海科技教育出版社，2001。
太阳系是稳定的吗？http://bolide.lamost.org/articles/article245.htm
三体周期轨道解动画http://numericaltank.sjtu.edu.cn/three-body/three-body-movies.htm
解读《三体》7. 三体运动：欧拉、拉格朗日特解与三体运动的新特解

本文使用程序

本文使用程序，保存后一键运行即可

function celestialbodymotion
% 天体运动求解
GM=1; % 万有引力常数，采用归一化的1(物理值6.67e-11)简化计算
PS.colors={'k' 'b' 'r' 'g' 'b' 'k' 'm' 'c'}; % 每个天体显示颜色
PS.linestyles={'-'}; % 每个天体自身轨迹线型
PS.markers={'none'}; % 每个天体标记
PS.frontmarkers={'o'}; % 动画中天体标记
PS.linewidthes={1 1 1 1 1}; % 轨迹粗细
icenter=0; % 绘图中心，0为采用绝对中心，其他数值代表采用此天体为中心
giffilename='star.gif'; % 输出动画名称
type=input(['1: 二体轨道\n'...
'2: 海王星轨道示意\n'...
'3: 混沌\n'...
'4: 拉格朗日点\n'...
'5: 8字形\n'...
'6: 复杂形状\n'...
'7: 毕达哥斯拉三角形\n'...
'8: 引力弹弓\n'...
'9: 洛希极限\n']);
% 定义天体质量、位置(三维坐标)、速度(三维坐标)。如果是限制性三体问题，则可将此天体质量设为0
switch(type)
case 1 % 二体轨道
Celestialbodies=[...
1e2 0 0 0 0 0 0
3 1 0 0 0 sqrt(GM*1e2/1) 0]; tt=5;icenter=0;
case 2 % 海王星示意
Celestialbodies=[...
1e2 0 0 0 0 0 0
3 1 0 0 0 sqrt(GM*1e2/1) 0
3 2 0 0 0 sqrt(GM*1e2/2) 0]; tt=5;icenter=2;
case 3 % 混沌示意图
v=[0. 0.3]; tt=30;
Celestialbodies=[...
1 -1 0 0 v 0
1 1 0 0 v.*[1 -1] 0
0 0 0 0 -2*v 0
];
case 4 % 欧拉解(L1,L2,L3)，拉格朗日解(L4,L5)
nL=input('第几个拉格朗日点?');
m1=5; m2=1; a=m2/(m1+m2); b=1-a;
r1=[-a 0 0]; r2=[b 0 0];
if(nL==1 || nL==2 || nL==3)
coeff=@(s1, s2) [1 3 3 1-s2 0 0]-a*[0 1 2 1+s1-s2 2*s1 s1];
ss=[-1 1; 1 1; -1 -1]; s=ss(nL,:);
u=roots(coeff(s(1), s(2)));
u=u(imag(u)==0);
r3=[u+b 0 0];
elseif(nL==4)
r3=[(b-a)/2 sqrt(3)/2 0];
elseif(nL==5)
r3=[(b-a)/2 -sqrt(3)/2 0];
end
omega=[0 0 sqrt(GM*(m1+m2)/(a+b)^3)];
v1=cross(omega, r1);
v2=cross(omega, r2);
v3=cross(omega, r3);
Celestialbodies=[...
m1 r1 v1
m2 r2 v2
0 r3 v3]; tt=4;icenter=[0 ];
case 5 % 8字形
m1=1; m2=1; m3=1;
r1=[-0.97000436, 0.24308753 0]; r3=-r1;
r2=[0 0 0];
v1=[0.4662036850, 0.4323657300 0]; v3=v1;
v2=[-0.93240737, -0.86473146 0];
Celestialbodies=[...
m1 r1 v1
m2 r2 v2
m3 r3 v3]; tt=10;icenter=0;
case 6 % 复杂形状
v=[0.6150407229 0.5226158545]; tt=40; % 蝴蝶型
% v=[0.4121028725 0.2833837497]; tt=40;
Celestialbodies=[...
1 -1 0 0 v 0
1 1 0 0 v 0
1 0 0 0 -2*v 0
]; icenter=0;
case 7 % 毕达哥斯拉三角形
Celestialbodies=[...
3 0 4 0 0 0 0
4 -3 0 0 0 0 0
5 0 0 0 0 0 0]; tt=100;icenter=0;
case 8 % 引力弹弓
Celestialbodies=[...
100 0 0 0 1 0 0
1 -1 -1 0 10 0 0]; tt=0.3;icenter=0;
case 9 % 洛希极限
PS.colors={'k'}; PS.frontmarkers=['o' repmat({'+'}, 1, 1000)]; PS.linestyles={'none'};
buf=linspace(-0.1, 0.1, 10); [x y]=meshgrid(buf, buf); x=x(:); y=y(:); O=x*0;
index=find(x.^2+y.^2<=0.01); x=x(index); y=y(index); O=O(index);
Celestialbodies=[...
100 0 0 0 1 0 0
1e-5/length(O)+O -1+x -1+y O 10+O O O]; tt=0.3; icenter=1;
end
%% 计算
% 初始值写成向量
n_bodies=size(Celestialbodies,1); % 天体数目
posvec0=reshape(Celestialbodies(:, 2:7), 1, []); % 积分初始位置和速度
options = odeset('RelTol',1e-8,'AbsTol', 1e-8*ones(1, n_bodies*6));
[TT YY]=ode45(@f_gravity, [0 tt], posvec0, options); % 积分
%% 绘图
% 复制标记，避免绘图时溢出
fns=fieldnames(PS);
for i=1:length(fns)
PS.(fns{i})=repmat(PS.(fns{i}), 1, n_bodies);
end
GifImCount=1; % 动画输出标记
if(icenter==0)
centerx=0*TT; centery=0*TT;
elseif(isnumeric(icenter))
centerx=mean(YY(:, icenter), 2); centery=mean(YY(:, icenter+n_bodies), 2);
end
% 绘制起始位置，以及运行轨迹
close all; figure('unit', 'centimeter', 'position', [0 0 14 12]); hold on;
for i=1:n_bodies
hm(i)=plot(YY(1, i)-centerx(1), YY(1, i+n_bodies)-centery(1), 'color', PS.colors{i}, 'marker', PS.frontmarkers{i});
plot(YY(:, i)-centerx, YY(:, i+n_bodies)-centery, 'linestyle', PS.linestyles{i}, 'color', PS.colors{i}, 'marker', PS.markers{i}, 'linewidth', PS.linewidthes{i});
end
set(gcf, 'color', 'w'); axis off; axis equal;
%% 制作动画并输出到gif
ht=title(sprintf('%.2f秒', 0)); % 显示时间
nfigs=100; nn=floor(size(YY, 1)/nfigs); if(nn==0) nn=1; nfigs=size(YY,1); end % 时间切片
for ifig = 1:nfigs
set(ht, 'string', sprintf('%.2f秒', TT(ifig*nn)));
for i=1:n_bodies % 设置新的位置
set(hm(i), 'xdata', YY(ifig*nn, i)-centerx(ifig*nn), 'ydata', YY(ifig*nn, i+n_bodies)-centery(ifig*nn));
end
frame(ifig) = getframe(gcf); % 结果输出至动画
im = frame2im(frame(ifig));
[imind,cm] = rgb2ind(im,256);
if GifImCount == 1
imwrite(imind,cm,giffilename,'gif','DelayTime',1/25, 'Loopcount',inf);
else
imwrite(imind,cm,giffilename,'gif','DelayTime',1/25, 'WriteMode','append');
end
GifImCount = GifImCount + 1;
end
movie(gcf, frame, 1, 25);
%% 万有引力方程，计算过程。
% 此处只是简单地采用Runge-Kutta数值积分算法进行轨道计算，此方法存在数值耗散，无法进行长时间轨道动力学模拟，因此无法给出图形，要进行长时间轨道动力学模拟，需要使用辛积分。
function d_posvec=f_gravity(t, posvec)
posvec=posvec(:);
% y: x1-xn, y1-yn, z1-zn, vx1-n, vy1-n ,vz1-n, 6nx1
n=length(posvec)/6; % 天体数量
d_posvec(1:3*n, 1)=posvec(3*n+1:end); % dx=vx
for i=1:n
accel=[0; 0; 0];
pos_self=posvec(i+[0 n 2*n]);
for j=1:n
pos_other=posvec(j+[0 n 2*n]);
m_other=Celestialbodies(j, 1);
dir=pos_other-pos_self; dis=norm(dir);
if(i~=j)
accel=accel+GM*m_other*dir/dis.^3;
end
end
d_posvec(i+[3*n 4*n 5*n])=accel;
end
end
end