P3945 三体问题

题目传送门

题目注意事项：

本题万有引力常数 $G$ 取 $6.67408 \times 10^{-11}$。

做这道题所需的知识：

1：距离公式

在三维空间内有两个点，设第一个点的坐标为 $(x_1,y_1,z_1)$，第二个点的坐标为 $(x_2,y_2,z_2)$，那两点之间的距离 distance 为：

$\operatorname{distance}=\sqrt{(x_2-x_1)^2+(y_2-y_1)^2+(z_2-z_1)^2}$。

2：算加速度的

设 $F$ 表示力的大小，$M_1,M_2$ 表示两个物体的质量，$r$ 表示两个物体之间的距离，$G$ 是万有引力常数，有 $F=G \times \dfrac{M_1 \times M_2}{r^2}$。

而且还有设 $F$ 表示力的大小，$M$ 表示物体的质量，$a$ 表示加速度 $F=M \times a$。

联立可得 $M·a=G \times \dfrac{M_1 \times M_2}{r^2}$。

化简可得 $a=G \times \dfrac{M}{r^2}$。

3：匀变速直线运动的公式

位移公式 $x=v_0t+\frac{1}{2}at^2$。

速度公式 $v=at$。

然后就可以开始构思代码了

先定义一个结构体来存天体。

struct stars{
	long double x,y,z;//位置 
	long double xv,yv,zv;//三个方向上的速度
	long double nowx,nowy,nowz;//现在的速度改变量 
	long double changex,changey,changez;//计算时的改变量
	long double tmp;
	long double mass;//质量 
}m[MAXN];

再解决一下改变量的计算。

m[i].changex=m[j].x-m[i].x; m[i].changey=m[j].y-m[i].y; m[i].changez=m[j].z-m[i].z;
//三个方向上的改变量 

m[i].tmp=sqrt(m[i].changex*m[i].changex + m[i].changey*m[i].changey + m[i].changez*m[i].changez);

m[i].nowx += m[j].mass*G/(m[i].tmp*m[i].tmp)/*距离*/ * m[i].changex/(m[i].tmp);
m[i].nowy += m[j].mass*G/(m[i].tmp*m[i].tmp)/*距离*/ * m[i].changey/(m[i].tmp);
m[i].nowz += m[j].mass*G/(m[i].tmp*m[i].tmp)/*距离*/ * m[i].changez/(m[i].tmp);
//引力和距离平方成反比,再乘一个change/tmp

最后再解决单位时间内改变的问题。

m[i].x+=t0*(m[i].xv+=t0*m[i].nowx);
m[i].y+=t0*(m[i].yv+=t0*m[i].nowy);
m[i].z+=t0*(m[i].zv+=t0*m[i].nowz);
//单位时间内的改变量 

再保留一个小数位数。

cout<<fixed<<setprecision(12);

最后就是完整代码。

#include <bits/stdc++.h>
#define G 6.67408e-11
//重力加速度 
using namespace std;

const long double t0=0.01;//单位时间 
const int MAXN=30+10;//最大值 

int n;//天体个数 
long double tt;//时间 

struct stars{
	long double x,y,z;//位置 
	long double xv,yv,zv;//三个方向上的速度
	long double nowx,nowy,nowz;//现在的速度改变量 
	long double changex,changey,changez;//计算时的改变量
	long double tmp;
	long double mass;//质量 
}m[MAXN];

int main(){
	cin>>n>>tt;
	for(int i=1;i<=n;i++){
		cin>>m[i].x>>m[i].y>>m[i].z;//现在的位置 
		cin>>m[i].mass;//质量 
		cin>>m[i].xv>>m[i].yv>>m[i].zv;
	}
	while(tt>0){
		tt-=t0;
		for(int i=1;i<=n;i++){
			m[i].nowx=m[i].nowy=m[i].nowz=0;
			for(int j=1;j<=n;j++){
				if(i!=j){
					m[i].changex=m[j].x-m[i].x; m[i].changey=m[j].y-m[i].y; m[i].changez=m[j].z-m[i].z;
					//三个方向上的改变量 
					m[i].tmp=sqrt(m[i].changex*m[i].changex + m[i].changey*m[i].changey + m[i].changez*m[i].changez);
					
					m[i].nowx += m[j].mass*G/(m[i].tmp*m[i].tmp)/*距离*/ * m[i].changex/(m[i].tmp);
					m[i].nowy += m[j].mass*G/(m[i].tmp*m[i].tmp)/*距离*/ * m[i].changey/(m[i].tmp);
					m[i].nowz += m[j].mass*G/(m[i].tmp*m[i].tmp)/*距离*/ * m[i].changez/(m[i].tmp);
					//引力和距离平方成反比,再乘一个change/tmp
				}
			}
		}
		for(int i=1;i<=n;i++){
			m[i].x+=t0*(m[i].xv+=t0*m[i].nowx);
			m[i].y+=t0*(m[i].yv+=t0*m[i].nowy);
			m[i].z+=t0*(m[i].zv+=t0*m[i].nowz);
			//单位时间内的改变量 
		}
	}
	cout<<fixed<<setprecision(12);
	for(int i=1;i<=n;i++) cout<<m[i].x<<' '<<m[i].y<<' '<<m[i].z<<endl;
	return 0;
}