1. 代码
现在为了把之前的代码延伸到三体系统，需要给常数增加一些东西——增加第三体的质量、位置和速率向量。把第三恒星的质量视作和太阳的质量等同。
#Mass of the Third Starm3=1.0 #Third Star#Position of the Third Starr3=[0,1,0] #mr3=sci.array(r3,dtype='float64')#Velocity of the Third Starv3=[0,-0.01,0]v3=sci.array(v3,dtype='float64')
需要更新代码中质心和质心速率的公式。#Update COM formular_com=(m1*r1+m2*r2+m3*r3)/(m1+m2+m3)#Update velocity of COM formulav_com=(m1*v1+m2*v2+m3*v3)/(m1+m2+m3)
对一个三体系统来说，需要修改运动方程使之包括另一物体施加的额外引力。因此，需要在RHS上，对问题中每一对物体施加力的其他物体增加一个力项。在三体系统的情况下，一个物体会受到其余两个物体施加的力的影响并因此在RHS上出现两个力项。数学上可表示为：
为在代码中反映这些变化，需要为odeint求解器创建一个新函数。
def ThreeBodyEquations(w,t,G,m1,m2,m3): r1=w[:3] r2=w[3:6] r3=w[6:9] v1=w[9:12] v2=w[12:15] v3=w[15:18] r12=sci.linalg.norm(r2-r1) r13=sci.linalg.norm(r3-r1) r23=sci.linalg.norm(r3-r2) dv1bydt=K1*m2*(r2-r1)/r12**3+K1*m3*(r3-r1)/r13**3 dv2bydt=K1*m1*(r1-r2)/r12**3+K1*m3*(r3-r2)/r23**3 dv3bydt=K1*m1*(r1-r3)/r13**3+K1*m2*(r2-r3)/r23**3 dr1bydt=K2*v1 dr2bydt=K2*v2 dr3bydt=K2*v3 r12_derivs=sci.concatenate((dr1bydt,dr2bydt)) r_derivs=sci.concatenate((r12_derivs,dr3bydt)) v12_derivs=sci.concatenate((dv1bydt,dv2bydt)) v_derivs=sci.concatenate((v12_derivs,dv3bydt)) derivs=sci.concatenate((r_derivs,v_derivs)) return derivs
最后，调用odeint函数并向其提供上述函数连同初始条件。#Package initial parametersinit_params=sci.array([r1,r2,r3,v1,v2,v3]) #Initial parametersinit_params=init_params.flatten() #Flatten to make 1D arraytime_span=sci.linspace(0,20,500) #20 orbital periods and 500 points#Run the ODE solverimport scipy.integratethree_body_sol=sci.integrate.odeint(ThreeBodyEquations,init_params,time_span,args=(G,m1,m2,m3))
给你写了一个求质心的函数，代码如下：（因为函数中使用到求平方根的函数sqrt，所以请包含math.h头文件）
#include <math.h>
POINT ZX(int X1,int Y1,int X2,int Y2,int X3,int Y3) //参数分别为三角形的三个坐标点
{float L1,L2,L3,N //L1,L2,L3分别代表三条边的长,(N用来作交换用)
POINT PN//用来表示质心的坐标
L1=sqrt((X1-X2)*(X1-X2)+(Y1-Y2)*(Y1-Y2)) //分别求出三条边的长
L2=sqrt((X1-X3)*(X1-X3)+(Y1-Y3)*(Y1-Y3))
L3=sqrt((X3-X2)*(X3-X2)+(Y3-Y2)*(Y3-Y2))
if (L1<L2) //如果L2比L1大,就把两个数交换
{N=L1
L1=L2
L2=N}
if (L1<L3) //如果L3比L1大,就把两个数交换
{N=L1
L1=L3
L3=N}
//经过两轮的比较和交换,可以确保L1是三条边中最大的一条
if (L1>=(L2+L3))//如果最大边大于等于两条小条的和,则三点构不成一个三角形
{PN.x=0xffffffff //设置一个错误值
PN.y=0xffffffff
return PN}//让函数返回错误值,这样调用函数之后就可以作出相应的判断
/*如果通过判断符合三角形的条件,求质心,质心就是重心,公式很简单.
就是X=(X1+X2+X3)/3Y=(Y1+Y2+Y3)/3,如果要证明有点长,这里就不说.
你可以自己试着证明一下,或百度一下*/
PN.x=(X1+X2+X3)/3
PN.y=(Y1+Y2+Y3)/3
return PN}
一维数据也是一样， 就是数组里面只有一个元素就行了，二维数据在画面上是平面，一维数据就是一根直线上的分布而已
你的数据改成[ [84],[78],[89]。。 ]这样的格式就可以用大部分算法跑了
