用python模拟三体运动.docx

用python模拟三体运动

第一篇：用python模拟三体运动用python模拟三体运动知乎不让传动图，大家可以移步我的博客wend.blog.ustc.edu.cn,上面有动图大四就要有大四的样子，最近闲来无事，想起来以前一直想做的一件事：模拟三体运动。正好可以练一练python。小试牛刀先从模拟二维正方形中，弹性小球的运动轨迹入手，熟悉相关操作。物理情景很简单，就是弹性小球位于正方形空间中，有某初速度，未与边界碰撞时进行匀速直线运动，与边界碰撞时发生弹性碰撞。模拟思路很简单，在t时间，位于r处，取小时间微元dt，认为小球在此微元内进行匀速直线运动，计算t+dt时刻的位置和速度，继续迭代。代码如下importnumpyasnpimportmatplotlib.pyplotaspltimportmatplotlib.animationasanimationdefsimData():dt=0.0001d1=0.25d2=0.25v1=12v2=5whileTrue:d1+=v1*dtd2+=v2*dtif(d1>1)or(d11)or(d2无限模拟下去，可以发现如果小球的初始条件不是那么特殊，小球最终将遍历正方形中所有点。这里吐槽一下python的动画模块，太不直观，帮助简直不知所云，在用动画展示数据这方面，mathematica就做的非常好。知乎不能发动图，所以我将图片存放在请输入提取码访问密码b684可以下载模拟三体为了简单，假设三个物体质量相同，有着相同的引力系数，三个物体的在t时刻的位置r、速度v都保存下来，取小时间微元dt，计算每个球所受的引力，然后计算加速度（按照xyz三个分量计算），推出t+dt时刻的位置和速度，继续迭代。代码如下importnumpyasnpimportmatplotlib.pyplotaspltimportmpl_toolkits.mplot3d.axes3dasp3importmatplotlib.animationasanimationx1=[]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=1000000foriinrange(5000000):d11+=v11*dtd12+=v12*dtd21+=v21*dtd22+=v22*dtd13+=v13*dtd23+=v23*dtd31+=v31*dtd32+=v32*dtd33+=v33*dtx1.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)+.1r13=pow(np.sqrt(pow(d11-d31,2)+pow(d12-d32,2)+pow(d13-d33,2)),3)+.1r23=pow(np.sqrt(pow(d31-d21,2)+pow(d32-d22,2)+pow(d33-d23,2)),3)+.1a121=g*(d21-d11)/r12a122=g*(d22-d12)/r12a123=g*(d23-d13)/r12a211=-a121a212=-a122a213=-a123a131=g*(d31-d11)/r13a132=g*(d32-d12)/r13a133=g*(d33-d13)/r13a311=-a131a312=-a132a313=-a133a321=g*(d21-d31)/r23a322=g*(d22-d32)/r23a323=g*(d23-d33)/r23a231=-a321a232=-a322a233=-a323v11=v11+(a121+a131)*dtv12=v12+(a122+a132)*dtv13=v13+(a123+a133)*dtv31=v31+(a321+a311)*dtv32=v32+(a322+a312)*dtv33=v33+(a323+a313)*dtv21=v21+(a211+a231)*dtv22=v22+(a212+a232)*dtv23=v23+(a213+a233)*dt'''fig=plt.figure()ax=fig.gca(projection='3d')ax.plot(x2,y2,z2)ax