用python模拟三体运动_三体运动python

用python模拟三体运动由刀豆文库小编整理，希望给你工作、学习、生活带来方便，猜你可能喜欢“三体运动python”。

用python模拟三体运动

知乎不让传动图，大家可以移步我的博客wend.blog.ustc.edu.cn,上面有动图大四就要有大四的样子，最近闲来无事，想起来以前一直想做的一件事：模拟三体运动。正好可以练一练python。小试牛刀先从模拟二维正方形中，弹性小球的运动轨迹入手，熟悉相关操作。物理情景很简单，就是弹性小球位于正方形空间中，有某初速度，未与边界碰撞时进行匀速直线运动，与边界碰撞时发生弹性碰撞。模拟思路很简单，在t时间，位于r处，取小时间微元dt，认为小球在此微元内进行匀速直线运动，计算t+dt时刻的位置和速度，继续迭代。代码如下import numpy as npimport matplotlib.pyplot as pltimport matplotlib.animation as animation def simData(): dt = 0.0001 d1 = 0.25 d2 = 0.25 v1 = 12 v2 = 5 while True: d1+=v1*dt d2+=v2*dt if(d1>1)or(d11)or(d2无限模拟下去，可以发现如果小球的初始条件不是那么特殊，小球最终将遍历正方形中所有点。这里吐槽一下python的动画模块，太不直观，帮助简直不知所云，在用动画展示数据这方面，mathematica就做的非常好。知乎不能发动图，所以我将图片存放在请输入提取码 访问密码 b684 可以下载

模拟三体为了简单，假设三个物体质量相同，有着相同的引力系数，三个物体的在t时刻的位置r、速度v都保存下来，取小时间微元dt，计算每个球所受的引力，然后计算加速度（按照xyz三个分量计算），推出t+dt时刻的位置和速度，继续迭代。

代码如下import numpy as npimport matplotlib.pyplot as pltimport mpl_toolkits.mplot3d.axes3d as p3import matplotlib.animation as animationx1=[]x2=[]y1=[]y2=[]z1=[]z2=[]x3=[]y3=[]z3=[]dt=0.00002d11 =-10d12 = 0d13=0v11 =-23v12 = 20v13=0d21 = 100d22 = 0d23=5v21 = 14v22 = 38v23=0d31=0d32=100d33=45v31=0v32=-30v33=34g=1000000for i in range(5000000): d11+=v11*dt d12+=v12*dt d21+=v21*dt d22+=v22*dt d13+=v13*dt d23+=v23*dt d31+=v31*dt d32+=v32*dt d33+=v33*dt x1.append(d11)y1.append(d12)z1.append(d13)x2.append(d21)y2.append(d22)z2.append(d23)x3.append(d31)y3.append(d32)z3.append(d33)r12=pow(np.sqrt(pow(d11-d21,2)+pow(d12-d22,2)+pow(d13-d23,2)),3)+.1 r13=pow(np.sqrt(pow(d11-d31,2)+pow(d12-d32,2)+pow(d13-d33,2)),3)+.1 r23=pow(np.sqrt(pow(d31-d21,2)+pow(d32-d22,2)+pow(d33-d23,2)),3)+.1 a121=g*(d21-d11)/r12 a122=g*(d22-d12)/r12 a123=g*(d23-d13)/r12 a211=-a121 a212=-a122 a213=-a123 a131=g*(d31-d11)/r13 a132=g*(d32-d12)/r13 a133=g*(d33-d13)/r13 a311=-a131 a312=-a132 a313=-a133 a321=g*(d21-d31)/r23 a322=g*(d22-d32)/r23 a323=g*(d23-d33)/r23 a231=-a321 a232=-a322 a233=-a323 v11=v11+(a121+a131)*dt v12=v12+(a122+a132)*dt v13=v13+(a123+a133)*dt v31=v31+(a321+a311)*dt v32=v32+(a322+a312)*dt v33=v33+(a323+a313)*dt v21=v21+(a211+a231)*dt v22=v22+(a212+a232)*dt v23=v23+(a213+a233)*dt'''fig = plt.figure()ax = fig.gca(projection='3d')ax.plot(x2, y2, z2)ax2 = fig.gca(projection='3d')ax2.plot(x1, y1, z1)ax3 = fig.gca(projection='3d')ax3.plot(x3, y3, z3)#ani = animation.FuncAnimation(fig, blit=False,interval=0.01, repeat=True)plt.show()'''输出动画的代码如下import matplotlib.pyplot as pltimport mpl_toolkits.mplot3d.axes3d as p3import matplotlib.animation as animationdef update_lines(num, dataLines, lines): for line, data in zip(lines, dataLines): # NOTE: there is no.set_data()for 3 dim data...line.set_data(data[0:2, :num])line.set_3d_properties(data[2, :num])return lines# Attaching 3D axis to the figurefig = plt.figure()ax = p3.Axes3D(fig)# Fifty lines of random 3-D lines#data = [Gen_RandLine(25, 3)for index in range(1)]data=[np.array([x1,y1,z1])[:,0:1000000:1000],np.array([x2,y2,z2])[:,0:1000000:1000],np.array([x3,y3,z3])[:,0:1000000:1000]]# Creating fifty line objects.# NOTE: Can't pa empty arrays into 3d version of plot()lines = [ax.plot(dat[0, 0:1], dat[1, 0:1], dat[2, 0:1])[0] for dat in data]# Setting the axes propertiesax.set_xlim3d([-200,200])ax.set_xlabel('X')ax.set_ylim3d([-200,200])ax.set_ylabel('Y')ax.set_zlim3d([-200,200])ax.set_zlabel('Z')ax.set_title('3D Test')# Creating the Animation objectline_ani = animation.FuncAnimation(fig, update_lines, fargs=(data, lines),interval=.1, blit=False)line_ani.save('3body2.gif')plt.show()由于gif能保存的信息太少，只能保存前段的运动情况，所以最后的运动情况给出静态图。最后的样子如图所示可以发现在初期运动相当没有规律，几乎陷入混沌，但是一段时间以后，有一个物体远离其他两个，剩下的两个纠缠在一起进行二体运动比翼齐飞，我试过不少初始条件，发现最后几乎都会发生这种情况，可能这就是三个质量差不多的物体的最终归宿吧。两个人相濡以沫，小三只好相忘于江湖。写到这里我突然发现，上周的GRE考试由个题大意是宇宙中没有质量相近的三星系统，但是却被观测到了，要你用逻辑解释这个问题。呵呵呵，直接用数学物理解释了。如果早早写程序的话，也许可以当场把模拟扔到GRE出题人身边。。

等有时间继续尝试一个超大质量物体和两个小质量物体的运动情况。另一组混乱之治镇楼。可以看到这组虽然初始的时候运动极其混乱，但是最后仍然是双星相伴，另外一个远离。

另外补充一句，如果随机产生初始条件的话，最终形成几个行星绕着一个大质量恒星转的概率基本为0，大多数情况是星体们擦肩而过，永不回头。可见我们生活在如此稳定的太阳系是多么的幸运。