推广了n体模拟中的力
generalizing the force in n-body simulation
本文关键字：模拟 更新时间：2023-10-16
在用c++编程n体问题时，我在推广引力方面遇到了一个问题。如果我试图重写我的程序，使其易于适应任意数量的物体，我在试图推广力时遇到了问题。
我写了以下代码：
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#define h 1000.0
#define G 6.67384*pow(10.0,-11)
using namespace std;
class particle{
public:
double kx1,kx2,kx3,kx4, kv1, kv2, kv3, kv4;
double ky1, ky2, ky3, ky4, kvy1, kvy2, kvy3, kvy4;
double x,y,vx,vy,m;
double dist(particle aap){
double dx = x - aap.x;
double dy = y - aap.y;
return sqrt(pow(dx,2.0)+pow(dy,2.0));
}
double g(double x1, double y1,particle aap){
return G*aap.m*(aap.x-x1)/pow(dist(aap),3.0);
}
double p(double x1, double y1, particle aap){
return G*aap.m*(aap.y-y1)/pow(dist(aap),3.0);
}
void update(){ //zet het object 1 stap vooruit
x = x + (1/6.0)*(kx1+2*kx2+2*kx3+kx4);
vx = vx + (1/6.0)*(kv1+2*kv2+2*kv3+kv4);
y = y + (1/6.0)*(ky1+2*ky2+2*ky3+ky4);
vy = vy + (1/6.0)*(kvy1+2*kvy2+2*kvy3+kvy4);
}
void create(double x1, double y1, double vx1, double vy1, double m1){
x = x1;
y = y1;
vx = vx1;
vy = vy1;
m =m1;
}
bool operator ==(particle &other){
if(x == other.x && y == other.y && vx == other.vx && vy == other.vy){
return true;
}
}
};
particle zon, maan, aarde;
void set(){
zon.create(1, 1, -2, 1, 2*pow(10.0,30));
aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx + -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
}
double xforce(double x1, double y1, particle aap){ //kracht in de x-richting
particle bodies[] = {zon, aarde, maan};
double fx;
for (int i = 0; i < 3; i++){
if (bodies[i].x == aap.x && bodies[i].y == aap.y && bodies[i].vx == aap.vx && bodies[i].vy == aap.vy ){;}
else{
fx += aap.g(x1,y1,bodies[i]);
}
}
return fx;
}
double yforce(double x1, double y1, particle aap){ //kracht in de y-richting
particle bodies[] = {zon, aarde, maan};
double fy;
for (int i = 0; i <= 3; i++){
if (bodies[i].x == aap.x && bodies[i].y == aap.y && bodies[i].vx == aap.vx && bodies[i].vy == aap.vy) {;}
else{
fy += aap.p(x1,y1,bodies[i]);
}
}
return fy;
}
void corr(particle& body){
body.kx1 = h*body.vx;
body.kv1 = h*xforce(body.x, body.y, body);
body.ky1 = h*body.vy;
body.kvy1 = h*yforce(body.x, body.y, body);
body.kx2 = h*(body.vx + 0.5*body.kv1);
body.kv2 = h*xforce(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, body);
body.ky2 = h*(body.vy + 0.5*body.kvy1);
body.kvy2 = h*yforce(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, body);
body.kx3 = h*(body.vx+ 0.5*body.kv2);
body.kv3 = h*xforce(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2, body);
body.ky3 = h*(body.vy+ 0.5*body.kvy2);
body.kvy3 = h*yforce(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2,body);
body.kx4 = h*(body.vx+body.kv3);
body.kv4 = h*xforce(body.x+ body.kx3, body.y + body.ky3, body);
body.ky4 = h*(body.vy + body.kvy3);
body.kvy4 = h*yforce(body.x + body.kx3, body.y + body.ky3, body);
}
void bereken(){
set();
zon.create(1, 1, -2, 1, 2*pow(10.0,30));
aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx + -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
ofstream file;
file.open("3body.txt");
for(int i =0; i <=30000; i++){
corr(maan);
corr(zon);
corr(aarde);
zon.update();
aarde.update();
maan.update();
file << i*h <<" "<< zon.x << " "<< zon.y << " "<< zon.vx<< " "<< zon.vy <<" "<< aarde.x << " " << aarde.y <<" "<< aarde.vx <<" " << aarde.vy <<" "<< maan.x<<" "<<maan.y<<"n";
}
file.close();
}
int main()
{
bereken();
system("pause");
return 0;
}
我认为问题出在函数xforce（）和yforce（。我不知道这里到底出了什么问题。我将在这里包含一个适用于3个实体的代码的工作版本，这样比较两者会更容易。
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#define h 1000.0
#define G 6.67384*pow(10.0,-11)
using namespace std;
class particle{
public:
double kx1,kx2,kx3,kx4, kv1, kv2, kv3, kv4;
double ky1, ky2, ky3, ky4, kvy1, kvy2, kvy3, kvy4;
double x,y,vx,vy,m;
double dist(particle aap){
double dx = x - aap.x;
double dy = y - aap.y;
return sqrt(pow(dx,2.0)+pow(dy,2.0));
}
double g(double x1, double y1,particle aap, particle bever){
return G*aap.m*(aap.x-x1)/pow(dist(aap),3.0) + G*bever.m*(bever.x-x1)/pow(dist(bever),3.0);
}
double p(double x1, double y1, particle aap, particle bever){
return G*aap.m*(aap.y-y1)/pow(dist(aap),3.0) + G*bever.m*(bever.y-y1)/pow(dist(bever),3.0);
}
void update(){ //zet het object 1 stap vooruit
x = x + (1/6.0)*(kx1+2*kx2+2*kx3+kx4);
vx = vx + (1/6.0)*(kv1+2*kv2+2*kv3+kv4);
y = y + (1/6.0)*(ky1+2*ky2+2*ky3+ky4);
vy = vy + (1/6.0)*(kvy1+2*kvy2+2*kvy3+kvy4);
}
void create(double x1, double y1, double vx1, double vy1, double m1){
x = x1;
y = y1;
vx = vx1;
vy = vy1;
m =m1;
}
bool operator ==(particle &other){
if(x == other.x && y == other.y && vx == other.vx && vy == other.vy){
return true;
}
}
};
particle zon, maan, aarde;
void set(){
zon.create(1, 1, -2, 1, 2*pow(10.0,30));
aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx + -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
}
double xforce(double x1, double y1, particle& aap){ //kracht in de x-richting
particle bodies[] = {zon, aarde, maan};
double fx;
for (int i = 0; i < 3; i++){
if (bodies[i].x == aap.x && bodies[i].y == aap.y && bodies[i].vx == aap.vx && bodies[i].vy == aap.vy ){;}
else{
fx += aap.g(x1,y1,bodies[i],aap);
}
}
return fx;
}
double yforce(double x1, double y1, particle& aap){ //kracht in de y-richting
particle bodies[] = {zon, aarde, maan};
double fy;
for (int i = 0; i <= 3; i++){
if (bodies[i].x == aap.x && bodies[i].y == aap.y && bodies[i].vx == aap.vx && bodies[i].vy == aap.vy) {;}
else{
fy += aap.p(x1,y1,bodies[i],aap);
}
}
return fy;
}
void corr(particle& body, particle aap, particle bever){
body.kx1 = h*body.vx;
body.kv1 = h*body.g(body.x, body.y, bever ,aap);
body.ky1 = h*body.vy;
body.kvy1 = h*body.p(body.x, body.y, bever, aap);
body.kx2 = h*(body.vx + 0.5*body.kv1);
body.kv2 = h*body.g(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, bever,aap);
body.ky2 = h*(body.vy + 0.5*body.kvy1);
body.kvy2 = h*body.p(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, bever,aap);
body.kx3 = h*(body.vx+ 0.5*body.kv2);
body.kv3 = h*body.g(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2, bever,aap);
body.ky3 = h*(body.vy+ 0.5*body.kvy2);
body.kvy3 = h*body.p(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2,bever,aap);
body.kx4 = h*(body.vx+body.kv3);
body.kv4 = h*body.g(body.x+ body.kx3, body.y + body.ky3, bever,aap);
body.ky4 = h*(body.vy + body.kvy3);
body.kvy4 = h*body.p(body.x + body.kx3, body.y + body.ky3, bever,aap);
}
void bereken(){
set();
zon.create(1, 1, -2, 1, 2*pow(10.0,30));
aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx + -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
ofstream file;
file.open("3body.txt");
for(int i =0; i <=30000; i++){
corr(maan, aarde, zon);
corr(zon, maan , aarde);
corr(aarde, maan , zon);
zon.update();
aarde.update();
maan.update();
file << i*h <<" "<< zon.x << " "<< zon.y << " "<< zon.vx<< " "<< zon.vy <<" "<< aarde.x << " " << aarde.y <<" "<< aarde.vx <<" " << aarde.vy <<" "<< maan.x<<" "<<maan.y<<"n";
}
file.close();
}
int main()
{
bereken();
system("pause");
return 0;
}
在这段代码中，总力是在函数g（）和p（）中计算的，而不是在xforce和yforce（）中。
您没有初始化变量。
例如：
double fy;
for (int i = 0; i <= 3; i++){
if (bodies[i].x == aap.x && bodies[i].y == aap.y && bodies[i].vx == aap.vx && bodies[i].vy == aap.vy)
{;}
else {
fy += aap.p(x1,y1,bodies[i]);
}
你应该做
double fy = 0;
相关文章：
- 如何使用Google Mock来模拟gettimeofday(）
- G锁定铸造到基础上会释放模拟行为
- 有什么好的方法可以让系统调用代理允许在单元测试中进行模拟
- 落砂模拟碰撞检测C++和SFML
- 在gtest.中使用fff.h模拟系统API
- 谷歌模拟和覆盖关键字
- 用C#中的并集模拟C++嵌套结构
- 在同一模拟中使用静脉和静脉_ inet内容时出现运行时错误
- 在模拟器中使用并集来模拟CPU寄存器有多合适
- 我写了一个C++程序来模拟Enigma机器.我没有得到输出
- 如何模拟不同边数的骰子滚动?
- 模拟持久按键
- 使用SIR模型的疾病爆发模拟
- 在 c++ 中模拟输入并在 JAVA 中读取它?
- 转发变量参数列表以模拟 std::thread
- 如何在谷歌模拟中匹配 C 样式数组
- 如何使用兰德随机化模拟点击
- 模拟GPS数据，以便使用Qt与Traccar一起使用
- QKeyPress - 在Qt中模拟按键
- 如何使用不同的谷歌模拟运行相同的谷歌测试用例?
