作为培训的一部分,我在C ++中实现了n-body类,以模拟物体的引力相互作用并更加熟悉C ++提供的功能,例如面向对象的编程。

此实现使用直接微分方程的积分(Verlet积分),导致时间复杂度为\ $ \ mathcal {O}(n ^ 2)\ $,其中\ $ n \ $是粒子数。

请在实施过程中尽可能地努力,并给我建设性的反馈。

我特别希望在以下领域提供建议:


代码样式(可读性,命名约定)
类设计
Efficieny(如何避免不必要的复杂性)
重新发明轮子(STL是否提供应在我的代码中使用的功能?)
内存使用情况

main.cpp

#include "nbody.h"

int main(int argc, char* argv[]) {
    Nbody nbody(16, 0.001, 1);
    nbody.timeIntegration();
    return 0;
}


nbody.h

#ifndef NBODY_H
#define NBODY_H

// Parameters
const int DIM = 2;          // dimensions
const double EPS = 1e-4;    // smoothing parameter

// Function prototypes
inline double sqr(double);

struct Particle{
    double m;           // mass
    double x[DIM];      // position
    double v[DIM];      // velocity
    double F[DIM];      // force 
    double F_old[DIM];  // force past time step
};

// Nbody class 
class Nbody {
    private:
        int step = 0;
        double t = 0;
        const int n;         // number of particles
        const double dt;           // step size
        const double t_max;        // max simulation time
        Particle *p = new Particle[n]; // allocate memory 
        void init_data();
    public:
        ~Nbody();
        Nbody(int n_, double dt_, double t_max_);
        inline void print_parameter() const;
        inline void print_data() const;
        inline void write_data(int step) const;
        void timeIntegration();
        void comp_force();
        void force(Particle*, Particle*);
        void comp_position();
        void comp_velocity();
        void update_position(Particle*);
        void update_velocity(Particle*);
};

#endif


nbody.cpp

#include <iostream>
#include <fstream> 
#include <cmath>
#include <random>

#include "nbody.h"

// Class methods
Nbody::Nbody(int n_, double dt_, double t_max_) : n(n_), dt(dt_), t_max(t_max_) {
    init_data();
}

Nbody::~Nbody() {
    delete[] p; 
    p = 0; 
}

void Nbody::timeIntegration() {
    comp_force();
    for(; t<t_max; t+=dt, step+=1) {
        comp_position();
        comp_force();
        comp_velocity();
        if (step % 10 == 0) {
            write_data(step);
            //print_data();
        }
    }
}

void Nbody::update_velocity(Particle *p) {
    double a = dt * 0.5 / p->m;
    for (int d=0; d<DIM; d++) {
        p->v[d] += a * (p->F[d] + p->F_old[d]);
    }
}

void Nbody::update_position(Particle *p) {
    double a = dt * 0.5 / p->m;
    for (int d=0; d<DIM; d++) {
        p->x[d] += dt * (p->v[d] + a * p->F[d]);
        p->F_old[d] = p->F[d];
    }
}

void Nbody::comp_velocity() {
    for (int i=0; i<n; i++) {
        update_velocity(&p[i]);
    }
}

void Nbody::comp_position() {
    for (int i=0; i<n; i++) {
        update_position(&p[i]);
    }
}

void Nbody::comp_force() {
    for (int i=0; i<n; i++) {
        for (int d=0; d<DIM; d++) {
            p[i].F[d] = 0;
        }
    }
    for (int i=0; i<n; i++) {
        for (int j=i+1; j<n; j++) {
            force(&p[i], &p[j]);
        }
    }
}

void Nbody::force(Particle *i, Particle *j) {
    double r=EPS; // smoothing
    for (int d=0; d<DIM; d++) {
        r += sqr(j->x[d] - i->x[d]);
    }
    double f = i->m * j->m / (sqrt(r) * r);
    for (int d=0; d<DIM; d++) {
        i->F[d] += f * (j->x[d] - i->x[d]);
        j->F[d] -= f * (j->x[d] - i->x[d]);
    }
}

void Nbody::write_data(int step) const {
    std::ofstream results;
    std::string file_name = "data_" + std::to_string(step) + ".log";
    results.open(file_name);
    if (results.fail()) { // or (!results) ?
        std::cerr << "Error\n" << std::endl;
    } else {
        for (int i=0; i<n; i++) {
            results << t << " ";
            results << p[i].m << " ";
            for (int d=0; d<DIM; d++) {
                results << p[i].x[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << p[i].v[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << p[i].F[d] << " ";
            }
            results << std::endl;
        }
        results.close();
    }
}

void Nbody::print_data() const {
    std::cout.setf(std::ios_base::scientific);
    std::cout.precision(5);
    for (int i=0; i<n; i++) {
        std::cout << t << " ";
        std::cout << p[i].m << " ";
        for (int d=0; d<DIM; d++) {
            std::cout << p[i].x[d] << " ";
        }
        for (int d=0; d<DIM; d++) {
            std::cout << p[i].v[d] << " ";
        }
        for (int d=0; d<DIM; d++) {
            std::cout << p[i].F[d] << " ";
        }
        std::cout << std::endl;
    }
}

void Nbody::init_data() {
    std::random_device rd;          
    std::mt19937 generator(rd()); 
    std::uniform_real_distribution<double> distribution_x(0.0,1.0);
    std::uniform_real_distribution<double> distribution_v(-1.0,1.0);
    for (int i=0; i<n; i++) {
        p[i].m = 1./n;
        for (int d=0; d<DIM; d++) {
            p[i].x[d] = distribution_x(generator);
            p[i].v[d] = distribution_v(generator);
            p[i].F[d] = 0.0;
            p[i].F_old[d] = 0.0;
        }
    }
}

inline void Nbody::print_parameter() const {
    std::cout << n << " " << dt << " " << t_max << std::endl;
}

// Other Functions

inline double sqr(double x) {
    return x*x;
}



评论

@Samuel如果您想提供更多信息,我认为对N体模拟是什么做一个解释会很棒。这样,人们就可以理解您的代码应该达到的目的。

您是否有要定位的特定C ++版本(C ++ 11,C ++ 14,C ++ 17)?

@pacmaninbw我不确定要使用哪个版本。有什么道理?编写C ++程序时应使用哪个版本?

您可以添加注释,但此时不要更改代码。此时您可能应该学习和使用的版本是C ++ 17。

@Samuel在考虑了我们的审核反馈并可能将其纳入您的代码之后,非常欢迎您在后续问题中发布它。我对这段代码的进一步发展很感兴趣。

#1 楼

太好了!
使用using namespace std;并不会犯初学者的基本错误! main()函数只有3行代码。
nbody类中不变的函数声明包括const,这将有助于以后的优化。
该代码利用C ++随机数生成而不是C srand()rand()函数。
由于Nbody是作为一个类实现的,因此更改main()非常容易,因此它可以接受用户输入的ndtt_max的值。
缺少标题
#include <string>中缺少nbody.cpp;在大多数情况下,这是编译代码所必需的。
已过时的
现在使用inline函数声明只是对编译器的建议。优化编译器可以并且将通过内联基于代码的代码来做更好的优化工作。
Nbody构造函数的主体使用过时的初始化形式,而不是像下面的代码中那样使用()
Nbody::Nbody(int n_, double dt_, double t_max_) : n(n_), dt(dt_), t_max(t_max_) {
    init_data();
}

使用花括号{}
Nbody::Nbody(int n_, double dt_, double t_max_)
: n{n_}, dt{dt_}, t_max{t_max_}
{
    init_data();
}

将初始化放在单独的行上可以更容易地找到。
首选STL容器类
优先使用STL容器类,例如std::vectorstd::array C样式数组。 std::array<type, size>类是固定大小的数组。 std::vector<type>是可变大小的阵列。 STL容器类提供了迭代器,因此不需要指针。使用std::vector<Particle> p;可能会减少构造函数的参数数量。绝对可以消除n类中对变量Nbody的需要,因为p.size()运行后始终包含粒子数。同样,在Nbody::init_data()运行之后,可以使用迭代器访问Nbody::init_data()中的粒子,并允许代码使用带范围的for循环,例如
void Nbody::write_data(int step) const {
    std::ofstream results;
    std::string file_name = "data_" + std::to_string(step) + ".log";
    results.open(file_name);
    if (results.fail()) { // or (!results) ?
        std::cerr << "Error\n" << std::endl;
    } else {
        for (auto particle : p) {
            results << t << " ";
            results << particle.m << " ";
            for (int d=0; d<DIM; d++) {
                results << particle.x[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << particle.v[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << particle.F[d] << " ";
            }
            results << std::endl;
        }
        results.close();
    }
}

使p成为STL容器类的另一个好处是,可以将类p的析构函数用作默认构造函数,并且无需在类声明中分配粒子数组。
变量名
仅通过阅读代码并不清楚变量Nbodyn_ndt_dtt_max_t_maxxFv是什么。例如,我假设p表示增量时间,但目前尚不清楚这是真的。数组dt可以重命名为p,如果我对particles的理解比对deltaTime更为合适。
是的,有些变量名有注释,但是如果我必须维护代码,我宁愿工作带有自记录代码而不是依赖注释。
示例
void Nbody::write_data(int step) const {
    std::ofstream results;
    std::string file_name = "data_" + std::to_string(step) + ".log";
    results.open(file_name);
    if (results.fail()) { // or (!results) ?
        std::cerr << "Error\n" << std::endl;
    } else {
        for (auto particle : particles) {
            results << t << " ";
            results << particle.mass << " ";
            for (int d=0; d<DIM; d++) {
                results << particle.position[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << particle.velocity[d] << " ";
            }
            for (int d=0; d<DIM; d++) {
                results << particle.Force[d] << " ";
            }
            results << std::endl;
        }
        results.close();
    }
}

样式
有些(不是全部),开发人员更喜欢在私有声明之前查看类的公共声明。一堂课这是因为查找类的公共接口变得更加容易。
除非计划要有多个构造函数,否则不需要dt函数,将代码移入构造函数可能更好。
如果功能void init_data()print_parameter()是调试功能,最好将它们放在print_data()#ifdef DEBUG中。
在当前实现中,不需要#endif中的return 0;。如果添加了错误处理代码并且有一个main(),则最好保留它。最好使用在return 1;return EXIT_SUCCESS;)中定义的EXIT_FAILUREcstdlib
建议
最好允许用户通过输入来命名结果进入的输出文件通过用户界面或作为命令行参数的一部分。如果用户未指定名称,则该名称可以默认为当前文件名。
最好只有一个输出文件。

评论


\ $ \ begingroup \ $
在我看来,对于这种领域,变量名具有很强的自我描述性。我认为假设n体交互程序的维护者具有物理学背景也是安全的。
\ $ \ endgroup \ $
– mcocdawc
19-10-23在21:41

\ $ \ begingroup \ $
你永远都不知道谁会继续维护你的代码。可能是一个完全不同的人。正确命名变量的成本可以忽略不计,但潜在的收益却很大。
\ $ \ endgroup \ $
– Hakaishin
19-10-24 9:00



\ $ \ begingroup \ $
@Hakaishin的代价是,任何具有领域知识的人都必须维护代码,每次他们阅读代码时,都必须从域的通用表示法$ dt $进行翻译。海事组织这使它的维护性较差。
\ $ \ endgroup \ $
– Pete Kirkham
19-10-24在11:15

\ $ \ begingroup \ $
@Hakaishin物理方程式使用数学符号(变量),而不是出于某种原因使用其名称和句子。使所有内容更具可读性。通过使用长名称,实现数学方程式(例如,数值方法)的代码变得更难解析。
\ $ \ endgroup \ $
–弗拉基米尔F
19-10-24在17:44

\ $ \ begingroup \ $
我已经看到了很多带有简洁的,看起来像数学的变量名的代码,这些变量名完全是伪造的(或者来自一篇晦涩的论文),而那些代码库却很糟糕。但是d代表变化,t代表时间,m代表质量,F代表力,等等,确实是这种基本的,众所周知的符号。
\ $ \ endgroup \ $
–hegel5000
19-10-24在19:30

#2 楼

首先

作为初学者,您做得很好。我已经编程了10年了,很长一段时间以来我的代码都比您编写的代码可读性低得多。就是说:

什么需要修复

我不了解n体问题的所有细节,但是我知道它的作用。我不是数值精度方面的专家,所以我不会对您执行的算术发表评论。从设计的角度来看,这是我要注意的几件事。

该类实际上无法测试

在构造输入数据随机化和使用一种方法进行输入之间在绝大多数工作中,为此类编写有意义的自动化测试非常困难。部分原因是该类做的太多了。

Public接口不能反映其用法

public接口比客户端使用的要广泛得多。据我所知,客户唯一需要做的就是构造这些对象之一,并立即对其调用timeIntegration(),然后以某种方式记录结果。

您将使用非标准方法来传达标准概念。

您将提供“ print_data”和“ write_data”方法。对于此类无需依赖<iostream><fstream>,这将使以自动化(阅读:单元测试)的方式进行测试非常困难。您应该为粒子类提供一个<<运算符,并允许客户端决定如何处理结果。

无法获取此类的原始数据

此外,由于print_data()write_data()方法似乎是从此类中获取数据的唯一方法,因此使用除了简单的命令提示符程序外,此类的其他类均受限制。一种以非打印形式获取内部数据的方法将很有帮助。

做什么

此类的更好设计是使用带有必要参数的公共构造函数,该参数立即调用计算集成所需的一切,然后使用一种方法来获取已处理的数据。没有其他东西可以公开了。这样,客户端很难正确使用此类。拥有自己唯一数据的具有吸气剂的类应该在OOP设计中提出一个危险信号,因此所有这些重新思考实际上导致了更大的认识,即...这不应该是class

我最大的考虑是根本不要上这个课。它拥有的所有数据都不会在有用的公共接口上保持不变。有关Wikipedia上类设计的不变量的更多信息。没有理由将引入的状态在整个生命周期中都归此类所有,并且有很多机会使用此类来生成完全无效的数据。相反,它应该具有由一个高级函数组成的接口。

n体计算器的公共接口应该包含两到三件事:


设置结构。除“热”数据外,这将包括所有必要的部分以正确运行计算。这将由客户端初始化。如果struct数据无效(即分母为零的东西),则函数应以某种返回码退出(如果环境允许,则返回异常)。这可能应该由const l值引用
通过std::vector<Particle>(可能是const l值)引用来获取,这是n身体计算器的输入数据,需要运行一个时间步。这可能是设置结构的一部分,但在我看来,它与设置结构中的其他概念截然不同。

此函数应保证在适当位置修改std::vector<Particle>或返回转换后的std::vector<Particle>。我的个人偏爱是后者,但是,取决于所使用的C ++版本,这可能会妨碍良好的性能。本质上,此功能所做的全部工作就是转换粒子状态列表。它可以(并且应该)使用其他辅助函数来完成其工作,并且这些函数很可能会在较大粒子框架的其他部分中重用。除了传入的粒子集之外,所有函数都应该是无状态的。

此多重折叠的附加值:


使用此函数更明显正确连接。参见最小惊讶原则。 Wiki文章。测试一组无状态函数比测试大型纠缠类要容易得多。
随着此代码库的扩展,这将允许更高程度的重用基本操作。

其他建议

名称

我会为Particle结构成员建议更好的名称。如果在较大的程序中正确使用它们,它们将很可能成为基本数据类型。输入质量,位置,速度和力没有错。的确,当您将位置称为x时,人们可能会知道您的意思,但是当您键入位置时,他们一定会知道您的意思。

强类型

我将对粒子成员使用强类型。乔纳森·博卡拉(Jonathan Bocarra)上有一些关于cppfluent的出色博客文章(例如CppFluent Strong类型)。可以将它们视为双精度型,其优点是使在函数调用中切换参数更加困难,并使代码更具表达性。

摆脱掉Globals

Globals是一件坏事,应该避免。不管是否摆脱了面向对象的方法,都应该将它们合并到某种设置结构中。

比您更多地使用STL

许多求和for循环可以使用std::accumulate();您应该使用std::vector而不是原始的c样式数组。您应该使用基于范围的for循环,而不能使用std::vector或STL算法。

评论


\ $ \ begingroup \ $
您完全误解了单一责任原则。
\ $ \ endgroup \ $
– gnasher729
19-10-24在10:28

\ $ \ begingroup \ $
你说得对!是的我将编辑我的答案
\ $ \ endgroup \ $
– Marc Olberding
19-10-24在13:13

\ $ \ begingroup \ $
谢谢!我认为很棒的评论。
\ $ \ endgroup \ $
–塞缪尔
19-10-25在20:25

#3 楼

除了其他答案之外:


DIMNbody.stepNbody.n使用无符号整数类型,因为这都不是负数;
由于C ++ 11而不使用constexpr仅对constDIM都使用EPS
摆脱argc中未使用的argvmain参数;
考虑更多使用const。例如,在f中的Nbody::force()可以是const,在a中的Nbody::update_position可以是const,依此类推。


评论


\ $ \ begingroup \ $
评论不用于扩展讨论;此对话已移至聊天。
\ $ \ endgroup \ $
– Mathieu Guindon♦
19-10-25在16:27

#4 楼

您的代码以C / C ++混合样式编写。例如,您的析构函数有一个delete(我找不到对应的new所在的位置),从根本上不需要。使用std::vector来存储类似数组的数据。

还可以像void Nbody::update_position(Particle *p)一样进行很多参数传递。请改用引用,如果只读取粒子,请使用const Particle &p。否则,对我来说,它看起来像n体代码。它是平方的,而不是更复杂/有效的东西,但这可能没问题。

哦,我发现了new:您在类定义中有
Particle *p = new Particle[n];
,但是n未初始化。这可能是未定义的行为,绝对是非常危险的,而且很可能是完全错误的。

不要使用new分配数组!使用std::vector,如下所示:

std::vector<Particle> the_particles;
public:
  Particles(int n) : the_particles(vector<Particle>(n)) {}
}```


评论


\ $ \ begingroup \ $
的确如此。作为编写答案的一部分,我将代码重构为使用std :: vector ,感觉好多了。
\ $ \ endgroup \ $
–罗兰·伊利格(Roland Illig)
19-10-23在20:13

\ $ \ begingroup \ $
我在类定义中使用new。我使用错误的方式吗?参数传递到底有什么问题?更新粒子的速度或位置时,是否不更改粒子?那么使用常量引用会不会出错?您能否更具体地说明comp_velocity()和update_velocity()的样子?
\ $ \ endgroup \ $
–塞缪尔
19-10-23在21:47

\ $ \ begingroup \ $
请您详细解释一下代码吗?我对STL不太熟悉。我是否将std :: vector the_particles定义为私有?粒子到底是什么?它在哪里定义?如何在其余代码中使用它?
\ $ \ endgroup \ $
–塞缪尔
19-10-25在7:28

\ $ \ begingroup \ $
我的编译器抛出一个错误:错误:'int'Particle(int n)之前的预期unqualified-id:p(vector (n)){};。看到我做错了吗?
\ $ \ endgroup \ $
–塞缪尔
19-10-25在7:33

\ $ \ begingroup \ $
1.它使用C ++ 17习惯用法。确保使用适当的编译器标志! 2.是,将向量设为私有。 3.我使用了自己的命名。要查找代码中的所有内容,需要太多的滚动。我确定你知道我的意思。
\ $ \ endgroup \ $
– Victor Eijkhout
19-10-25在15:48



#5 楼

除了其他答案:

init_data函数不属于Nbody类。在N体问题的定义中,找不到“随机”一词,使用随机输入数据仅与您的特定情况有关,因此应将此代码移至main.cpp中。

Nbody的构造函数,参数名称中不需要结尾的下划线。以下代码看起来更简洁,并且等效于您当前的代码:

 Nbody::Nbody(int n, double dt, double t_max)
: n(n), dt(dt), t_max(t_max) {
    init_data();  // should be removed, as I said above
}
 


用于调试不仅可以使用timeIntegration方法,而且还可以使用仅执行一个步骤的简单step方法,这将是一个很好的选择。这使您可以编写更好的单元测试。这也使另一个构造函数参数t_max变得不必要。

仍然在timeIntegration中,而不是step+=1,您应该编写++step。编写step++是等效的,但这将告诉每个读者您不太了解C ++。在C ++中,++通常位于变量之前,在其他语言(如Java或C或Go)中,通常位于变量之后。有关更多详细信息,请参见此Stack Overflow答案。

timeIntegrationupdate_velocity的代码进行比较,发现您使用的编程风格不一致。您应该自己决定是否使用camelCase或snake_case标识符。然后,始终使用该样式。另一件事是,您在运算符*/周围放置了空格,但不在+周围放置了空格。我本来希望如此,因为*/+更紧密地绑定操作数。通常的样式是始终用空格包围二进制运算符。因此t < t_max; t += dt; step++

您的Nbody类无法解决棘手的情况,在这些棘手的情况下,粒子之间的距离非常近,以至于dt对于实际仿真而言变得太大。这是必须记录的东西。

我喜欢将updated_velocityupdate_position分成两个单独的方法。这使它们易于阅读。 (另外,从实现的角度出发,这是必要的,因为必须先更新所有粒子的速度,然后才能更新任何粒子的位置,否则结果取决于粒子的顺序。)

comp中的缩写comp_position不明确。它可能意味着比较或计算。

Nbody::force中,您不应命名参数ij,因为按照惯例,这些变量名是为整数保留的。我宁愿选择p和q。如果将Nbody::p重命名为ps,因为它无论如何都是复数形式,那么就不再存在命名冲突。您只需删除参数即可。

方法write_data应该称为step,因为它涉及的是所有参数,而不仅仅是一个参数。

在API级别,我会不要将Nbody::stepprint_parameter放在构造函数中,而是将print_parameters作为参数传递给dt方法,而将t_max作为参数传递给dt方法。在step中有一​​个t_max常量,看起来很可疑。对于0.001的timeIntegration,它可能具有适当的nbody.h值,但是如果我想使用EPS模拟该怎么办?我认为这不应该是一个全局常数。甚至不应该是光速,因为取决于精确的实验,光速是如此之多。

dt中,您写的0.0001没有尾随0。当然,它可以保存一个按键,但是我认为这是不值得的。只需编写规范的dt = 1.0e-9即可,就像您在同一函数中的其他几个地方所做的一样。

写入Nbody::init_data文件的数据非常不准确。典型的1.类型提供16到17位的精度,但是您只写出了6位,这是C ++的默认值。自2017年以来,C ++终于支持准确打印浮点数了。

评论


\ $ \ begingroup \ $
C和Java都像C ++一样具有前递增和后递增(我不知道Go,但如果有所不同,将会感到惊讶)。可以说,当值将被丢弃时,我们通常更喜欢预增量(尽管您忽略了在for循环的增量表达式中编写的位置-为什么?),并且对于非固定类型可能更有效,但是-increment在C ++中确实有它的位置,应在需要时使用。也许该段可以改写?
\ $ \ endgroup \ $
– Toby Speight
19-10-24在7:02

\ $ \ begingroup \ $
@TobySpeight谢谢,做到了。
\ $ \ endgroup \ $
–罗兰·伊利格(Roland Illig)
19-10-24在17:37

\ $ \ begingroup \ $
感谢您的好评!在n-body类中也删除尾随下划线是否正确?还是这是一个不好的风格?您将dt添加到step方法是什么意思?那你要把EPS放在哪里?
\ $ \ endgroup \ $
–塞缪尔
19-10-25在19:51

\ $ \ begingroup \ $
是的,可以在所有地方删除下划线。关于step方法:我将定义Nbody :: step(double dt,double eps)。这样可以完全控制呼叫者。
\ $ \ endgroup \ $
–罗兰·伊利格(Roland Illig)
19-10-26在6:34

\ $ \ begingroup \ $
您说过init_data()不属于Nbody类。您能在这点上详细说明一下吗?那我该如何初始化粒子呢?将init_data()添加到粒子结构是否也有意义?
\ $ \ endgroup \ $
–塞缪尔
19-10-28的16:16

#6 楼

使用向量数学库

找到一个合适的库来实现坐标向量,因此您不必将它们实现为双精度数组。理想情况下,您的struct Particle应该看起来像:

struct Particle {
    double m;   // mass
    vec3 x;     // position
    vec3 v;     // velocity
    vec3 F;     // force
    vec3 F_old; // force past time step
};


,合适的库将提供函数和运算符重载,以使使用这些类型变得非常容易。您应该能够编写如下内容:

void Nbody::update_position(Particle *p) {
    double a = dt * 0.5 / p->m;
    p->x += dt * (p->v + a * p->F);
    p->F_old = p->F;
}


有很多可用的库。我本人也喜欢GLM。有关可能的库的讨论,请参见https://stackoverflow.com/questions/1380371/what-are-the-most-widely-used-c-vector-matrix-math-linear-algebra-libraries-a。 />使函数操纵Particle的成员函数


您有很多函数主要操纵粒子的状态,但它们不是Particle本身的一部分。例如,struct Particle是除时间步之外的东西,它仅操纵update_position()的成员变量。如果将其设置为dt的成员函数,则它将变得更加简洁:

struct Particle {
    ...
    void update_position(double dt);
};

void Particle::update_position(double dt) {
    double a = dt * 0.5 / m;
    x += dt * (v + a * F);
    F_old = F;
}


您可以这样称呼:

void Nbody::comp_position() {
    for (auto &p: particles) {
        p.update_position(dt);
    }
}


Particle甚至Particle都可以做同样的事情。

评论


\ $ \ begingroup \ $
我尝试首先结合您答案的第二部分,但出现以下错误:nbody.cpp:在成员函数'void Particle :: update_position(double)'中:nbody.cpp:47:33:error:'在此范围内未声明p'const double a = dt * 0.5 / p-> m; ^ nbody.cpp:在成员函数“ void Nbody :: comp_position()”中:nbody.cpp:55:19:错误:未在此范围内为“(自动&p:粒子){
\ $ \ endgroup \ $
–塞缪尔
19-10-25在7:51

\ $ \ begingroup \ $
是建议将其更好地用于(自动&p:粒子)还是(粒子&p:粒子)?什么时候不应该使用自动?
\ $ \ endgroup \ $
–塞缪尔
19-10-25在15:36

\ $ \ begingroup \ $
我建议您在上下文中已经明确类型或不关心类型时使用auto。在这里,很明显p是质点,因为它是质点的元素。在这些情况下,使用auto可以避免不必要的重复类型,并可以避免错误。
\ $ \ endgroup \ $
– G. Sliepen
19-10-25在16:14

\ $ \ begingroup \ $
请您概述一下如何用方法计算力?我遇到了困难,因为在力计算的情况下,我有一个嵌套循环,因此我不确定在这种情况下该如何做。
\ $ \ endgroup \ $
–塞缪尔
19-10-25在16:24

\ $ \ begingroup \ $
将void force(Particle&other)写为struct Particle的成员函数后,将其中的force(&p [i],&p [j])更改为p [i] .force(p [j])嵌套循环。在这里,您可能应该继续使用int i和int j作为迭代器变量,因为这是确保每个对仅被访问一次的最简单方法。
\ $ \ endgroup \ $
– G. Sliepen
19-10-25在18:28

#7 楼

由于这里是新手,所以我无法评论,但是Roland Illig断言它应该是++step而不是step++,并且它表明您不了解C ++是不正确的。

在C ++中, ++确定如何计算表达式的顺序。因此,在++step中,变量在执行任何操作之前先递增,而在step++中,变量在执行值之前执行。仅将step++++step作为一行代码基本上是等效的,但是在这样的示例中,区别是显而易见的:

int step = 0;
std::cout << ++step << std::endl; // would print 1
std::cout << step << std::endl; // would print 1


while

int step = 0;
std::cout << step++ << std::endl; // would print 0
std::cout << step << std::endl; // would print 1


只是要澄清一下,因为您应该理解区别,而不是出于风格/声誉原因而不是一个选择另一个!

评论


\ $ \ begingroup \ $
我认为Rolan Illig表示++ i通常比i ++快。
\ $ \ endgroup \ $
–eanmos
19-10-24在6:52

\ $ \ begingroup \ $
@eanmos,这是可能的,但是在您提供的链接中指出,只有当要递增的类型不是原始类型时,这种性能影响通常才存在。但是在OP的帖子中,step的类型为int,编译器应正确对其进行优化。
\ $ \ endgroup \ $
– ashiswin
19-10-24在7:03

\ $ \ begingroup \ $
是的,但是默认情况下编写++ i并仅在需要使用后增量时才使用它是一个好习惯。
\ $ \ endgroup \ $
–eanmos
19-10-24在7:11

\ $ \ begingroup \ $
@eanmos足够公平。更多的断言是:“在C ++中,++出现在变量之前,在其他语言中,例如Java或C或Go,则在变量之后。”但是,是的,约定是使用++ i。
\ $ \ endgroup \ $
– ashiswin
19-10-24在7:12



\ $ \ begingroup \ $
@Samuel,简而言之。 i ++必须创建旧值的副本,“返回”它,然后递增i。同时++ i不必创建副本并“返回”它,它会使i递增,而只是“返回”它。对于内置类型(例如int),编译器可能能够优化您的代码。但是我可以是C ++类的实例,因此i ++和++ i可以调用operator ++函数之一。在那种情况下,编译器可能无法优化,然后在使用i ++时复制对象i可能真的很慢。因此,在一般情况下,只需使用++ i。
\ $ \ endgroup \ $
–eanmos
19-10-24在10:42



#8 楼

我将专注于已经由另一个答案解决的一件事,但我认为值得更多关注:单一责任原则。

您的NBody类具有合并为一个的多个功能,建议分开。据我所知:


它代表一组N粒子
它提供了进行物理模拟的算法
它提供了实现打印模拟结果

我认为有足够的材料将它们分成三个单独的实体,为将来的更改留出了更大的灵活性。

此外,其中一些NBody类中的方法实际上仅作用于给定的Particle,因此可以将它们重构为Particle结构的方法。

另一个建议是看看模板方法模式,它可以是模拟框架的有用起点,可以在必要时提供适当的灵活性来更改集成方法。

评论


\ $ \ begingroup \ $
你是说我分开写不同的课吗?粒子类,模拟类和打印结果的类?
\ $ \ endgroup \ $
–塞缪尔
19-10-24在17:08

#9 楼

除了G. Sliepen的想法之外,您还可以使用STL的std::valarray<double>。这样一来,您可以将

for (int d = 0; d < DIM; ++d) {
    p->x[d] += dt * (p->v[d] + a * p->F[d]);
    p->F_old[d] = p->F[d];
}

替换为

p->F_old = p->F;
p->x += dt * (p->v + a * p->F);


布置数组结构而不是结构数组。如果粒子多于尺寸,则可以让您对所有x坐标,然后对所有y坐标和所有z坐标执行更宽的矢量运算,而不仅限于坐标系的宽度。也就是说,每个p可能只有两个或三个并行计算,但是如果您有多个std::array<std::valarray<double>, DIM>,其中q坐标在x[0]中,则x坐标在x[1]中,z坐标在x[2]中,速度在v[0]中,等等。,可能看起来像:

for (size_t i = 0; i < x.size(); ++i) {
  F_old[i] = F[i];
  x[i] += dt * (v[i] + a * F[i]);
}


,并且能够使用向量寄存器的全宽。但是,如果计算不是那么容易分离的话,这将无法正常工作。