问题描述
我正在尝试使用Runge Kutta 4集成算法编写代码来解决n体问题。 我正在使用质量在0和1之间均匀分布的两个物体测试代码，位置遵循与1 / r ^ 2成比例的密度定律分布，并且速度按照麦克斯韦-玻耳兹曼分布进行分配。我曾尝试针对不同的tmax集成系统，但是我在曲线图中得到了轨道，但我无法找出问题所在。 任何帮助将不胜感激。
这是代码：
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
typedef struct {
double x,y,z;
}vector;
double *masse;
vector *pos,*vel,*forza;
int main(int argc,char** argv){
double epsilon,dt,tmax,t;
double dx,dy,dz,dist,invdist,invdist3;
int N,i,l,m,n,s,j;
double a,b;
vector *k1,*k2,*k3,*k4,*w1,*w2,*w3,*w4,*pos1,*vel1;
if(argc!=5) {
fprintf(stdout,"Il programma prende in input il softening,il passo d'integrazione,il tempo massimo d'integrazione e il numero di corpi del sistema\n",argv[0]);
exit(1);
}
epsilon=strtod(argv[1],NULL);
dt=strtod(argv[2],NULL);
tmax=strtod(argv[3],NULL);
N=strtod(argv[4],NULL);
FILE* fp=fopen("Cond_ini.out","r");
if(fp==NULL){
perror("Errore: file non trovato\n");
exit(1);
}
masse=(double*)malloc(N*sizeof(double));
pos=(vector*)malloc(N*sizeof(vector));
vel=(vector*)malloc(N*sizeof(vector));
forza=(vector*)malloc(N*sizeof(vector));
k1=(vector*)malloc(N*sizeof(vector));
k2=(vector*)malloc(N*sizeof(vector));
k3=(vector*)malloc(N*sizeof(vector));
k4=(vector*)malloc(N*sizeof(vector));
w1=(vector*)malloc(N*sizeof(vector));
w2=(vector*)malloc(N*sizeof(vector));
w3=(vector*)malloc(N*sizeof(vector));
w4=(vector*)malloc(N*sizeof(vector));
pos1=(vector*)malloc(N*sizeof(vector));
vel1=(vector*)malloc(N*sizeof(vector));
for(i=0;i<N;i++){
fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf",&masse[i],&pos[i].x,&pos[i].y,&pos[i].z,&vel[i].x,&vel[i].y,&vel[i].z);
}
fclose(fp);
printf("Condizioni iniziali:\n");
for(l=0;l<N;L++){
printf("%lf %lf %lf %lf %lf %lf %lf\n",masse[l],pos[l].x,pos[l].y,pos[l].z,vel[l].x,vel[l].y,vel[l].z);
}
for(t=0;t<tmax;t+=dt){
for(m=0;m<N;m++){
for(n=0;n<N;n++){
if(m!=n){
k1[n].x=dt*vel[n].x;
k1[n].y=dt*vel[n].y;
k1[n].z=dt*vel[n].z;
dx=pos[n].x-pos[m].x;
dy=pos[n].y-pos[m].y;
dz=pos[n].z-pos[m].z;
dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
invdist=(1/sqrt(dist));
invdist3=(invdist*invdist*invdist);
forza[n].x=dx*invdist3*masse[n];
forza[n].y=dy*invdist3*masse[n];
forza[n].z=dz*invdist3*masse[n];
w1[n].x=dt*forza[n].x;
w1[n].y=dt*forza[n].y;
w1[n].z=dt*forza[n].z;
pos1[n].x=pos[n].x+(0.5*k1[n].x);
pos1[n].y=pos[n].y+(0.5*k1[n].y);
pos1[n].z=pos[n].z+(0.5*k1[n].z);
vel1[n].x=vel[n].x+(0.5*w1[n].x);
vel1[n].y=vel[n].y+(0.5*w1[n].y);
vel1[n].z=vel[n].z+(0.5*w1[n].z);
k2[n].x=dt*(vel[n].x+(0.5*w1[n].x));
k2[n].y=dt*(vel[n].y+(0.5*w1[n].y));
k2[n].z=dt*(vel[n].z+(0.5*w1[n].z));
dx=pos1[n].x-pos[m].x;
dy=pos1[n].y-pos[m].y;
dz=pos1[n].z-pos[m].z;
dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
invdist=(1/sqrt(dist));
invdist3=(invdist*invdist*invdist);
forza[n].x=dx*invdist3*masse[n];
forza[n].y=dy*invdist3*masse[n];
forza[n].z=dz*invdist3*masse[n];
w2[n].x=dt*forza[n].x;
w2[n].y=dt*forza[n].y;
w2[n].z=dt*forza[n].z;
pos1[n].x=pos[n].x+(0.5*k2[n].x);
pos1[n].y=pos[n].y+(0.5*k2[n].y);
pos1[n].z=pos[n].z+(0.5*k2[n].z);
vel1[n].x=vel[n].x+(0.5*w2[n].x);
vel1[n].y=vel[n].y+(0.5*w2[n].y);
vel1[n].z=vel[n].z+(0.5*w2[n].z);
k3[n].x=dt*(vel[n].x+(0.5*w2[n].x));
k3[n].y=dt*(vel[n].y+(0.5*w2[n].y));
k3[n].z=dt*(vel[n].z+(0.5*w2[n].z));
dx=pos1[n].x-pos[m].x;
dy=pos1[n].y-pos[m].y;
dz=pos1[n].z-pos[m].z;
dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
invdist=(1/sqrt(dist));
invdist3=(invdist*invdist*invdist);
forza[n].x=dx*invdist3*masse[n];
forza[n].y=dy*invdist3*masse[n];
forza[n].z=dy*invdist3*masse[n];
w3[n].x=dt*forza[n].x;
w3[n].y=dt*forza[n].y;
w3[n].z=dt*forza[n].z;
pos1[n].x=pos[n].x+(k3[n].x);
pos1[n].y=pos[n].y+(k3[n].y);
pos1[n].z=pos[n].z+(k3[n].z);
vel1[n].x=vel[n].x+(w3[n].x);
vel1[n].y=vel[n].y+(w3[n].y);
vel1[n].z=vel[n].z+(w3[n].z);
k4[n].x=dt*(vel[n].x+w3[n].x);
k4[n].y=dt*(vel[n].y+w3[n].y);
k4[n].z=dt*(vel[n].z+w3[n].z);
dx=pos1[n].x-pos[m].x;
dy=pos1[n].y-pos[m].y;
dz=pos1[n].z-pos[m].z;
dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
invdist=(1/sqrt(dist));
invdist3=(invdist*invdist*invdist);
forza[n].x=dx*invdist3*masse[n];
forza[n].y=dy*invdist3*masse[n];
forza[n].z=dy*invdist3*masse[n];
w4[n].x=dt*forza[n].x;
w4[n].y=dt*forza[n].y;
w4[n].z=dt*forza[n].z;
a=k1[n].x+(2*k2[n].x)+(2*k3[n].x)+k4[n].x;
a=a/6;
pos1[n].x=pos[n].x+a;
a=k1[n].y+(2*k2[n].y)+(2*k3[n].y)+k4[n].y;
a=a/6;
pos1[n].y=pos[n].y+a;
a=k1[n].z+(2*k2[n].z)+(2*k3[n].z)+k4[n].z;
a=a/6;
pos1[n].z=pos[n].z+a;
b=w1[n].x+(2*w2[n].x)+(2*w3[n].x)+w4[n].x;
b=b/6;
vel1[n].x=vel[n].x+a;
b=w1[n].y+(2*w2[n].y)+(2*w3[n].y)+w4[n].y;
b=b/6;
vel1[n].y=vel[n].y+a;
b=w1[n].z+(2*w2[n].z)+(2*w3[n].z)+w4[n].z;
b=b/6;
vel1[n].z=vel[n].z+a;
}
}
for(j=0;j<N;j++) {
forza[j].x=0;
forza[j].y=0;
forza[j].z=0;
}
}
for(i=0;i<N;i++){
pos[i].x=pos1[i].x;
pos[i].y=pos1[i].y;
pos[i].z=pos1[i].z;
vel[i].x=vel1[i].x;
vel[i].y=vel1[i].y;
vel[i].z=vel1[i].z;
printf("%lf %lf %lf %lf %lf %lf %lf\n",t,pos[i].x,pos[i].y,pos[i].z,vel[i].x,vel[i].y,vel[i].z);
/*forza[i].x=0;
forza[i].y=0;
forza[i].z=0;*/
}
}
}
这是不同tmax时的轨道图：
解决方法
浏览代码时，第一个错误是算法错误。循环的顺序必须是时间步长，然后是RK4阶段，其内是恒星体相互作用项/作用力的总和。这尤其意味着您需要分别对Runge-Kutta方法的每个阶段分别计算交互作用。
进行从位置数组填充力数组的单独过程以避免复制粘贴编辑错误是很有意义的。
有关这种记录良好的测试问题，请参阅IVP测试套件https://archimede.dm.uniba.it/~testset/testsetivpsolvers/?page_id=26的Pleiades系统