洛谷 P3945 三体问题

洛谷 P3945 三体问题

在物竞dalao的帮助下（简化下？）终于A了此题，于是在他的提议下来喷出题人。

接下来看题。

题意分析

模拟三维空间中\(n\)个星体的运动，求\(Ts\)后\(n\)个星体的坐标。使用微元法，\(dt=0.01s\)。

\(n\le 30,0\le T\le 100\)。

精度要求：\(0.5\%\)。

Solution

以下变量大多都是矢量，就不打箭头了

首先列出我们在小学一年级就学过的万有引力公式

\[F=\frac{GMm}{r^2} \]

其中\(M,m\)分别是两颗星球的质量，\(r\)是距离。

那么在每时每刻，用\(n\)个星体的质量和位置就可以算出每对星球之间的相互作用力。

设\(F_i\)表示第\(i\)颗星球的受力，则

\[F_i=\sum_{j=1,j\not =i}^n \frac{Gm_im_j}{r^2}=Gm_i\cdot\sum_{j=1,j\not =i}^n \frac{m_j}{r^2} \]

那么，根据牛顿第二定律\(F=am\)，这一刻的加速度为

\[a_i=\frac{F_i}{m_i}=G\sum_{j=1;j\not =i}^n \frac{m_j}{r^2} \]

注意这里前面有个\(\times m_i\)，这里有个\(\div m_i\)，在程序中可以不用计算了，以保证精度。

然后再用我们在幼儿园就学过的基本运动学公式

\[x=v_0t+\frac12at^2 \\ v=at \]

当您高高兴兴的算到了这里，觉得这道题配不起它是道蓝题的时候——

欸！我怎么才\(38pts\)？

难道是要先更新\(v\)再更新\(x\)？（即，用\(t+dt\)时刻的速度\(v+a\cdot dt\)当作\(t\)时刻的初速度去计算位移\(x\)）

欸！用错误的做法反而分更高了？！（\(40pts\)）

再开个\(\texttt{long double}\)试试——

欸！才\(50pts\)？难道要手写实数高精？

看看样例，居然样例跑的都这么吃劲（差距肉眼可见！）

当然也不是控制输出多少位小数的问题，因为错误只会出在\(0.005\)中。

于是我请\(LHQing\)来帮我调这一道题（他以前写过一个物理模拟器），在10min后他过了。

原因是，此题\(std\)认为在\(dt\)中星球做匀速直线运动，速度只会在每个\(dt\)之后改变,即

\[v=at \\ x=vt \]

只需要这两个公式，不需要考虑星球在这段\(dt\)中的加速度！

那么，从某种意义上，\(std\)的精度岂不是比我的要劣……

Code

三维矢量

struct vector
{
	long double x,y,z;
	vector(long double xx=0,long double yy=0,long double zz=0)
	{
		x=xx,y=yy,z=zz;
	}
	void clear()
	{
		x=y=z=0;
	}
	void in()
	{
		cin>>x>>y>>z;
	}
	void out()
	{
		printf("%Lf %Lf %Lf",x,y,z);
	}
	vector operator+(const vector& b)const
	{
		return vector(x+b.x,y+b.y,z+b.z);
	}
	vector operator+=(const vector& b)
	{
		x+=b.x;
		y+=b.y;
		z+=b.z;
		return *this;
	}
	vector operator*(const long double& b)
	{
		return vector(x*b,y*b,z*b);
	}
	vector operator*=(const long double& b)
	{
		x*=b;
		y*=b;
		z*=b;
		return *this;
	}
	vector operator/(const long double& b)const
	{
		return vector(x/b,y/b,z/b);
	}
}F[100];

星球

struct Star
{
	vector pos,v;
	long double m;
	IL void in()
	{
		pos.in();
		cin>>m;
		v.in();
	}
	IL void out()
	{
		pos.out();
	}
	
}star[100];

计算

for(;T>0;T-=dt)
{
	for(int i=1;i<=n;i++) F[i].clear();
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
		{
			if(i==j) continue;
			r=(star[i].pos.x-star[j].pos.x)*(star[i].pos.x-star[j].pos.x)+(star[i].pos.y-star[j].pos.y)*(star[i].pos.y-star[j].pos.y)+(star[i].pos.z-star[j].pos.z)*(star[i].pos.z-star[j].pos.z);
			f=G*star[j].m/r;
			F[i]+=vector(f*(star[j].pos.x-star[i].pos.x)/sqrt(r),f*(star[j].pos.y-star[i].pos.y)/sqrt(r),f*(star[j].pos.z-star[i].pos.z)/sqrt(r));
		}
	}
	for(int i=1;i<=n;i++)
	{
		star[i].v+=F[i]*dt;
		star[i].pos+=(star[i].v)*dt;
	}
}

完整代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<limits.h>
#define IL inline
#define re register
#define LL long long
#define ULL unsigned long long
#ifdef TH
#define debug printf("Now is %d\n",__LINE__);
#else
#define debug
#endif
using namespace std;
const long double G = 6.67408e-11;
const long double dt = 1e-2;
int n;
long double T;
struct vector
{
	long double x,y,z;
	vector(long double xx=0,long double yy=0,long double zz=0)
	{
		x=xx,y=yy,z=zz;
	}
	void clear()
	{
		x=y=z=0;
	}
	void in()
	{
		cin>>x>>y>>z;
	}
	void out()
	{
		printf("%Lf %Lf %Lf",x,y,z);
	}
	vector operator+(const vector& b)const
	{
		return vector(x+b.x,y+b.y,z+b.z);
	}
	vector operator+=(const vector& b)
	{
		x+=b.x;
		y+=b.y;
		z+=b.z;
		return *this;
	}
	vector operator*(const long double& b)
	{
		return vector(x*b,y*b,z*b);
	}
	vector operator*=(const long double& b)
	{
		x*=b;
		y*=b;
		z*=b;
		return *this;
	}
	vector operator/(const long double& b)const
	{
		return vector(x/b,y/b,z/b);
	}
}F[100];
struct Star
{
	vector pos,v;
	long double m;
	IL void in()
	{
		pos.in();
		cin>>m;
		v.in();
	}
	IL void out()
	{
		pos.out();
	}
	
}star[100];
int main()
{
	cin>>n>>T;
	for(int i=1;i<=n;i++) star[i].in();
	long double f,r;
	for(;T>0;T-=dt)
	{
		for(int i=1;i<=n;i++) F[i].clear();
		for(int i=1;i<=n;i++)
		{
			for(int j=1;j<=n;j++)
			{
				if(i==j) continue;
				r=(star[i].pos.x-star[j].pos.x)*(star[i].pos.x-star[j].pos.x)+(star[i].pos.y-star[j].pos.y)*(star[i].pos.y-star[j].pos.y)+(star[i].pos.z-star[j].pos.z)*(star[i].pos.z-star[j].pos.z);
				f=G*star[j].m/r;
				F[i]+=vector(f*(star[j].pos.x-star[i].pos.x)/sqrt(r),f*(star[j].pos.y-star[i].pos.y)/sqrt(r),f*(star[j].pos.z-star[i].pos.z)/sqrt(r));
			}
		}
		for(int i=1;i<=n;i++)
		{
			star[i].v+=F[i]*dt;
			star[i].pos+=(star[i].v)*dt;
		}
	}
	for(int i=1;i<=n;i++)
	{
		star[i].out();
		cout<<endl;
	}
	return 0;
}

注意第\(95\)行，\(r\)的定义是\(r^2\)，不先开方是为了精度。

以及玄学的物理计算在第\(102,103\)行

喷喷喷