P3945 三体问题

题意

题目描述

给定 \(n\) 个天体的位矢 \((r_{ix},r_{iy},r_{iz})\)、初速度 \((v_{ix},v_{iy},v_{iz})\) 与质量 \(m_i\)，万有引力常数 \(G\) 取 \(6.67408\times 10^{-11}\)，求 \(t\) 时刻所有天体的位矢（天体视为质点，单位都为 SI 单位）。

为了简化计算，假设时间的流逝是不连续的，每 \(0.01\) 秒刷新一次（下文中，也称 \(0.01\) 秒称为一单位时间）；并且在一单位时间内，视天体的运动为匀速直线运动。

输入格式

第一行一个正整数 \(n\) 和一个非负实数 \(t\)（\(3\le n\le30,0\le t\le100\)），含义见题目描述。
接下来 \(n\) 行，每行 \(7\) 个实数 \(r_{ix},r_{iy},r_{iz},m_i,v_{ix},v_{iy},v_{iz}\)（\(-100\le r_{ix},r_{iy},r_{iz}\le 100\)，\(m_i\) 在 long long 范围内），含意见题目描述。

输出格式

共 \(n\) 行，每行 \(3\) 个实数 \(r'_{ix},r'_{iy},r'_{iz}\)，表示 \(t\) 时刻第 \(i\) 个天体的位矢。
本题开启 Special Judge，标准答案保留 \(12\) 位小数，当你的答案与标准答案的相对误差不超过 \(0.5\%\) 时视为正确。

思路

首先列出本题目背景下需要用到的物理公式：

\(F=\frac{Gm_1m_2}{r^2}\)，\(G\) 为万有引力常数，\(m_1\)、\(m_2\) 为两天体的质量，\(r\) 为两天体之间的距离。
\(a=\frac{F}{m}\)，\(F\) 为天体受到的力，\(m\) 为天体的质量。
\(v=v_0+at\)，\(v\) 为末速度，\(v_0\) 为初速度，\(a\) 为加速度，\(t\) 为经过的时间。
\(r=r_0+vt\)，\(r\) 为末位矢，\(r_0\) 为初位矢，\(v\) 为末速度，\(t\) 为经过的时间（现实世界中应当是 \(r=r_0+\frac{1}{2}at^2\)，但由于本题中时间流逝不连续，且在一单位时间内天体的运动视为匀速直线运动，因此修改为上式）。

然后就可以开始模拟了：

在一单位时间内，对于每一天体 \(i\)，计算其他每一天体对其的万有引力所产生的加速度的矢量和；
依照上述公式更新每一天体的位矢以及速度。
重复上述步骤，直至到达时刻 \(t\)。

实现时需要注意的几个细节：

由于相对误差不得超过 \(0.5\%\)，因此浮点数类型采用 long double；
计算天体 \(i\) 对天体 \(j\) 的万有引力所产生的加速度时，可以先计算加速度 \(a\) 的大小，在计算其在 \(x\) 轴、\(y\) 轴、\(z\) 轴上的分量。具体内容详见代码实现。
将每一天体天体的加速度都计算好后，再更新每一天体的位矢和速度。

程序

AC 记录

#include<cstdlib>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<queue>
#include<cstdio>
#include<iostream>
#include<vector>
#include<map>
#include<cmath>
#include<iomanip>
#include<string>
#include<stack>
#define re register
#define ll long long
#define ull unsigned long long
#define vl __int128
#define ld long double
#define LL 2e18
#define INT 1e9
#define INF 0x3f3f3f3f
#define lb(x) (x&(-x))
#ifdef __linux__
#define gc getchar_unlocked
#define pc putchar_unlocked
#else
#define gc _getchar_nolock
#define pc _putchar_nolock
#endif
int T=1;
using namespace std;
inline bool blank(const char x){return !(x^32)||!(x^10)||!(x^13)||!(x^9);}
template<typename Tp>inline void read(Tp &x){x=0;re bool z=true;re char a=gc();for(;!isdigit(a);a=gc())if(a=='-')z=false;for(;isdigit(a);a=gc())x=(x<<1)+(x<<3)+(a^48);x=(z?x:~x+1);}
inline void read(double &x){x=0.0;re bool z=true;re double y=0.1;re char a=gc();for(;!isdigit(a);a=gc())if(a=='-')z=false;for(;isdigit(a);a=gc())x=x*10+(a^48);if(a!='.')return x=z?x:-x,void();for(a=gc();isdigit(a);a=gc(),y/=10)x+=y*(a^48);x=(z?x:-x);}
inline void read(ld &x){x=0.0;re bool z=true;re ld y=0.1;re char a=gc();for(;!isdigit(a);a=gc())if(a=='-')z=false;for(;isdigit(a);a=gc())x=x*10+(a^48);if(a!='.')return x=z?x:-x,void();for(a=gc();isdigit(a);a=gc(),y/=10)x+=y*(a^48);x=(z?x:-x);}
inline void read(char &x){for(x=gc();blank(x)&&(x^-1);x=gc());}
inline void read(char *x){re char a=gc();for(;blank(a)&&(a^-1);a=gc());for(;!blank(a)&&(a^-1);a=gc())*x++=a;*x=0;}
inline void read(string &x){x="";re char a=gc();for(;blank(a)&&(a^-1);a=gc());for(;!blank(a)&&(a^-1);a=gc())x+=a;}
template<typename T,typename ...Tp>inline void read(T &x,Tp &...y){read(x),read(y...);}
template<typename T>inline void read(T *begin,T *end){re T *i;if(begin<end)for(i=begin;i<end;++i)read(*i);else for(i=begin-1;i>=end;--i)read(*i);}
template<typename Tp>inline void write(Tp x){if(!x)return pc(48),void();if(x<0)pc('-'),x=~x+1;re int len=0;re char tmp[64];for(;x;x/=10)tmp[++len]=x%10+48;while(len)pc(tmp[len--]);}
inline void write(const double x){re int a=6;re double b=x,c=b;if(b<0)pc('-'),b=-b,c=-c;re double y=5*powl(10,-a-1);b+=y,c+=y;re int len=0;re char tmp[64];if(b<1)pc(48);else for(;b>=1;b/=10)tmp[++len]=floor(b)-floor(b/10)*10+48;while(len)pc(tmp[len--]);pc('.');for(c*=10;a;a--,c*=10)pc(floor(c)-floor(c/10)*10+48);}
inline void write(const ld x){re int a=6;re ld b=x,c=b;if(b<0)pc('-'),b=-b,c=-c;re ld y=5*powl(10,-a-1);b+=y,c+=y;re int len=0;re char tmp[64];if(b<1)pc(48);else for(;b>=1;b/=10)tmp[++len]=floor(b)-floor(b/10)*10+48;while(len)pc(tmp[len--]);pc('.');for(c*=10;a;a--,c*=10)pc(floor(c)-floor(c/10)*10+48);}
inline void write(const pair<int,double>x){re int a=x.first;if(a<7){re double b=x.second,c=b;if(b<0)pc('-'),b=-b,c=-c;re double y=5*powl(10,-a-1);b+=y,c+=y;re int len=0;re char tmp[64];if(b<1)pc(48);else for(;b>=1;b/=10)tmp[++len]=floor(b)-floor(b/10)*10+48;while(len)pc(tmp[len--]);a&&(pc('.'));for(c*=10;a;a--,c*=10)pc(floor(c)-floor(c/10)*10+48);}else printf("%.*lf",a,x.second);}
inline void write(const pair<int,ld>x){re int a=x.first;if(a<7){re ld b=x.second,c=b;if(b<0)pc('-'),b=-b,c=-c;re ld y=5*powl(10,-a-1);b+=y,c+=y;re int len=0;re char tmp[64];if(b<1)pc(48);else for(;b>=1;b/=10)tmp[++len]=floor(b)-floor(b/10)*10+48;while(len)pc(tmp[len--]);a&&(pc('.'));for(c*=10;a;a--,c*=10)pc(floor(c)-floor(c/10)*10+48);}else printf("%.*Lf",a,x.second);}
inline void write(const char x){pc(x);}
inline void write(const bool x){pc(x?49:48);}
inline void write(char *x){fputs(x,stdout);}
inline void write(const char *x){fputs(x,stdout);}
inline void write(const string &x){fputs(x.c_str(),stdout);}
template<typename T,typename ...Tp> inline void write(T x,Tp ...y){write(x),write(y...);}
template<typename T>inline void write(T *begin,T *end,const char c=' '){re T *i;for(i=begin;i<end;++i)write(*i,c);}
template<typename T>inline void init(T *begin,T *end,const T& val=T()){re T* i;for(i=begin;i<end;++i)*i=val;}
template<typename T>inline T max(T *begin,T *end){re T *ans,*i;for(i=begin;i<end;++i)if(i==begin||*ans<*i)ans=i;return *ans;}
template<typename T>inline T min(T *begin,T *end){re T *ans,*i;for(i=begin;i<end;++i)if(i==begin||*i<*ans)ans=i;return *ans;}
template<typename T>inline T calc_sum(T *begin,T *end,const T& val=T()){re T ans=val,*i;for(i=begin;i<end;++i)ans+=*i;return ans;}
template<typename T>inline bool is_equal(T *begin,T *end,const T& val=T()){re T *i;for(i=begin;i<end;++i)if(*i!=val)return false;return true;} //模板，无需在意。 

ll mod=0;
const int MAX=30;
const int N=MAX+10;
const ld G=6.67408e-11;
const ld dt=0.01;
//#define DEBUG
//#define more_text
struct planet{
	ld rx,ry,rz; //位矢 
	ld vx,vy,vz; //速度 
	ld ax,ay,az; //加速度 
	ld m;        //质量 
};
inline void read(planet &p){read(p.rx,p.ry,p.rz,p.m,p.vx,p.vy,p.vz);}
inline void write(planet p){write(make_pair(12,p.rx),' ',make_pair(12,p.ry),' ',make_pair(12,p.rz));}

int n;ld t;planet p[N];

void solve(int step){
	read(n,t);
	read(p+1,p+n+1);
	
	while(t>0){
		for(int i=1;i<=n;++i){
			p[i].ax=p[i].ay=p[i].az=0;
			
			for(int j=1;j<=n;++j){
				if(i==j)continue;
				
				ld dx=p[j].rx-p[i].rx,
				   dy=p[j].ry-p[i].ry,
				   dz=p[j].rz-p[i].rz; //计算由 p[i] 的位置引向 p[j] 的位置的矢量在 x 轴、y 轴、z 轴上的分量。 
				
				ld d=sqrt(dx*dx+dy*dy+dz*dz); //计算 p[i] 与 p[j] 之间的距离 。 
				
				ld a=G*p[j].m/(d*d); //a=F/p[i].m
				                     // =(G*p[i].m*p[j].m*p[j].m/(d*d))/p[i].m
				                     // =G*p[j].m/(d*d)
				p[i].ax+=a*dx/d;
				p[i].ay+=a*dy/d;
				p[i].az+=a*dz/d; //计算加速度在  x 轴、y 轴、z 轴上的分量。
			}
		}
		
		for(int i=1;i<=n;++i){
			p[i].vx+=p[i].ax*dt;
			p[i].vy+=p[i].ay*dt;
			p[i].vz+=p[i].az*dt; //计算速度在  x 轴、y 轴、z 轴上的分量。
			
			p[i].rx+=p[i].vx*dt;
			p[i].ry+=p[i].vy*dt;
			p[i].rz+=p[i].vz*dt; //计算位矢在  x 轴、y 轴、z 轴上的分量。
		}
		
		t-=dt;
	}
	
	write(p+1,p+n+1,'\n');
}
int main(){ //模板，无需在意。 
	#ifdef DEBUG
	freopen("test.in","r",stdin);freopen("test.out","w",stdout);
	#endif
	#ifdef more_text
	read(T);
	#endif 
	for(int i=0;i<T;++i)solve(i);
	#ifdef DEBUG
	fclose(stdin);fclose(stdout);
	#endif
	return 0;
}
/*
Input:
3 100
0 10 0 10000000 0.006207480877613 0 0
8.660254037844 -5 0 10000000 -0.003103740438807 -0.00537583613352 0
-8.660254037844 -5 0 10000000 -0.003103740438807 0.00537583613352 0
Output:
0.620349511786 9.980741705470 0.000000000000
8.333401109655 -5.527609289167 0.000000000000
-8.953750621441 -4.453132416303 0.000000000000
Outline:

*/

posted @ 2025-03-09 09:53 LXcjh4998 阅读(35) 评论(0) 收藏举报

刷新页面返回顶部

LXcjh4998

P3945 三体问题

题意

题目描述

输入格式

输出格式

思路

程序

公告