c - C中的N体模拟
问题描述
我正在尝试使用 Runge Kutta 4 集成算法编写代码来解决 n 体问题。我正在使用质量均匀分布在 0 和 1 之间的两个物体测试代码,位置分布遵循与 1/r^2 成比例的密度定律,速度分布为 Maxwell-Boltzmann 分布。我试图为不同的 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 测试套件的 Pleiades 系统https://archimede.dm.uniba.it/~testset/testsetivpsolvers/?page_id=26
参照。还有如何使用函数更新 c++ 类成员?我在其中发布了实现这些原则/更正的代码,用于使用 Verlet 方法(和 C++)进行模拟的测试问题。
推荐阅读
- spring-boot - 我应该如何更改 JUNIT 中布尔变量的值?
- openssl - 从 base64 字符串中的证书获取可分辨名称(DN)
- python - 无法导入 QUERY_TERMS
- python - 如何使用python将超链接添加到ppt中的图像
- sigmoid - 为什么 Relu 能解决梯度消失问题?
- javascript - 带有工具提示的 FusionCharts
- angular - 在 Angular 8 中避免使用 Mangle 类名称
- sql-server - 为什么 SQL Server 2019 Developer 安装失败?
- mule - Mule 4中黑白foreach,并行foreach和批处理的主要区别是什么
- javascript - 如何有目的地地址的方向(谷歌地图API)