模拟粒子运动

Question

描述

该代码进行粒子模拟。在这个代码中，一些粒子以一定的分布长度采取步骤。如果粒子超过一条曲线，代码就会把它们放回去。

码

   #include <iostream>
#include <vector>
#include <fstream>
#include <random>
#include <cmath>
using namespace std;
double randnum (double aa, double bb)  //defining a function to create random numbers
{
  static std::default_random_engine generator;
  std::uniform_real_distribution<double> distribution (aa,bb);
  return distribution(generator);
}

int main() {
    double dr;
    double a = 1.0;               //minimum value for generated random number with power-law distribution
    int N = 100000;                  //number of time steps
    double u;                    //uniform random number used to create random number with desired distribution
    int Particles = 1000;       //number of particles
    int H = 500;                    //potential parameter
    int h = 200;                    //potential parameter
    double C = 0.004;                  //potential parameter
    double A = 0.002;              //Potential Parameter
    double Delta = 1/7;            //Assymetric Parameter
    double b = 25.0;               //maximum value for generated random number with power-law distribution
    double  phi;   //first angle
    std::vector<double>  x(Particles); 
    std::vector<double>  pre_x(Particles);  //x_copmpnent of i'th particle position in previous step
    std::vector<double>  pre_y(Particles);
    std::vector<double>  J(N); 
    std::vector<double>  y(Particles); 
    for(int i = 0; i < Particles; i++)  //Initializing initial angle
        {
        pre_x[i] = 0;
        pre_y[i] = 200;
        } 

    ofstream myfile;
    myfile.open ("J.txt");
    if(alph < 1.01 && alph > 0.99)
        B = 1  //renormalazation constant for distribution
        for(int i = 0; i < N; i++)  //time steps
            {
            J[i] = 0;
            for (int j = 0; j < Particles; j++)  //Particles
                {
                 u = randnum(0,1);
                 dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
                 phi = randnum(0,M_PIl);
                 x[j] = pre_x[j] + cos(phi) * dr;
                 y[j] = pre_y[j] + sin(phi) * dr;
                 while( (sin(A * x[j]) + Delta * sin(C * x[j])/2) * h + H < y[j] || y[j] < 0)
                     {
                     u = randnum(0,1);
                     dr = pow( pow( a, 1-alph ) + u * (1-alph)/B, 1/(1-alph));
                     phi = randnum(0,M_PIl);
                     x[j] = pre_x[j] + cos(phi) * dr;
                     y[j] = pre_y[j] + sin(phi) * dr;
                     }
                 pre_x[j] = x[j];
                 pre_y[j] = y[j];
                 J[i] = J[i] + cos(phi);
                 }
                     myfile<<J[i]/Particles<<endl; //Outputs array to txtFile
            }
            myfile.close();
    }

关注点：

• 代码看起来很慢
• 通用代码质量

它是缓慢的，换句话说，我认为有可能使它更快。我需要"N“和”粒子“都是大值。在这种情况下，如何使代码运行得更快？有什么变化可以改善运行时间吗？任何其他反馈也是受欢迎的。

Answer 1

以下是一些可以帮助您改进代码的东西。

不要滥用using namespace std

将using namespace std放在每个程序的顶端是您最好避免的坏习惯。一旦删除了该语句，就必须在适当的地方添加名称空间。例如，fstream将改为std::fstream，endl将更改为std:endl，等等。

将代码分解为较小的函数

在当前代码中，几乎所有事情都是用main()完成的。与其把所有东西都放在一个长函数中，如果每一个离散的步骤都是它自己的功能，它将更容易阅读和维护。

消除了“幻数”

这段代码中到处都是“魔术数字”，即未命名的常量，如0.004、200、25.0等。一般来说，最好避免这种情况，并给这些常量取有意义的名称。这样，如果任何事情都需要更改，您就不必对所有" 200“实例的代码进行搜索，然后尝试确定这个特定的200是一个需要更改的对象，还是一个恰好具有相同值的其他常量。

考虑改进名称

我认为randnum不是一个糟糕的函数名，但大多数变量名称都是单个字母名称，根本无法描述它们要表示的内容。更多描述性变量名称将帮助代码读者了解它正在做什么。

将结果写入std::cout

与其有一个硬编码的文件名，不如将数据写入std::cout。如果用户随后想保存副本，那么这样做只是命令行输出重定向的简单问题。

使用<cmath>而不是<math.h>

使用新样式的<cmath>而不是C样式的<math.h>有两个原因。首先，它是更惯用的现代C++，也是因为它使用了名称空间。与此相一致，您还应该使用std::cos，而不仅仅是cos。log、pow和sin也是如此。

不使用std::endl，当'\n'做

时

使用std::endl会发出一个\n并刷新流。除非您确实需要流刷新，否则可以通过简单地发射'\n'来提高代码的性能，而不是使用可能计算成本更高的std::endl。

只使用必要的#includes

#include <functional>行不是必需的，可以安全地移除。

在实际的

中使用constexpr

在main中，几乎所有的变量实际上都被用作常量，所以如果使用符合C++11的编译器，那么至少将它们声明为const (最好是constexpr )是有意义的(如果没有，则确实应该更新编译器)。

使用对象

代码的目的是模拟称为粒子的物体。它是用C++写的。为什么代码不使用名为C++的Particle对象？下面是一个我们可以使用的Particle类：

class Particle {
public:
    Particle() : x_{0}, y_{200}, phi_{0} {}
    Particle &operator++();
    double phi() const { return phi_; }
public:
    static double get_dr(); 
    void next(const Particle &other);
    bool bad() const;
    double x_;
    double y_;
    double phi_;
};

现在main可以如下所示：

int main() {
    constexpr int N{10000};
    std::vector<double> J(N, 0);
    constexpr int numParticles{1000};
    std::vector<Particle> particles(numParticles);

    for (auto &step : J) {
        for (auto &p : particles) {
            ++p;
            step += std::cos(p.phi());
        }
        std::cout << step / numParticles << '\n';
    }
}

这不是很好吗？

简化了数学

让我们看看上面的类函数是如何实现的。首先是operator++：

Particle &Particle::operator++() {
    Particle testp;
    for (testp.next(*this); testp.bad(); testp.next(*this)) 
    { /* do nothing */ }
    return *this = testp;
}

这在原始代码中具有内循环的效果。next函数如下：

void Particle::next(const Particle &other) {
    double dr = get_dr();
    phi_ = randnum(0, M_PIl);
    x_ = static_cast<int>(other.x_) + std::cos(phi_) * dr;
    y_ = static_cast<int>(other.y_) + std::sin(phi_) * dr;
}

get_dr函数如下：

double Particle::get_dr() {
    constexpr double a = 1.0;
    constexpr double b = 25.0;
    constexpr int alph = 0;
    constexpr double B = (alph < 1.01 && alph > 0.99) ?
        1 / (std::log(b) * std::pow(b, 1 - alph) - std::log(a) * std::pow(a, 1 - alph))
       : (1 - alph) / (std::pow(b, (1 - alph)) - std::pow(a, (1 - alph)));
    return std::pow(std::pow(a, 1 - alph) + randnum(0, 1) * (1 - alph) / B, 1 / (1 - alph));
}

最后，bad()的简单版本如下：

bool Particle::bad() const {
    constexpr int H = 500;
    constexpr int h = 200;
    constexpr double C = 0.004;
    constexpr double A = 0.002;
    constexpr double Delta = 1 / 7;
    return 
        (std::sin(A * x_) + Delta * std::sin(C * x_) / 2) * h + H < y_;
}

然而，sin的最大值是1，所以这个不等式的左手边的上界可以表示为常数。由于与常量比较通常比调用std::sin要快得多，所以我们可以如下所示使用这些知识：

bool Particle::bad() const {
    constexpr int H = 500;
    constexpr int h = 200;
    constexpr double C = 0.004;
    constexpr double A = 0.002;
    constexpr double Delta = 1 / 7;
    constexpr double maxval = (1 + Delta * 1 / 2) * h + H;
    return maxval < y_ || 
        (std::sin(A * x_) + Delta * std::sin(C * x_) / 2) * h + H < y_;
}

短路表达式求值意味着我们只调用计算量大的sin函数，即与maxval的比较失败。在我的机器上的测试表明，这使程序大约快了20%。您可以挤出更多的时间，注意到只有当第一个测试失败时才需要计算x。通常，您希望安排一些事情，以便首先进行更简单的测试。

如果结果方程简化，那么考虑使用不同的坐标方案(例如，极坐标和矩形坐标)可能是有利的。我怀疑他们可能会。

Answer 2

您的代码很难读懂。这主要是因为大多数变量都有神秘的名称。

在物理模拟中，可以用名字来命名自然常数，比如c或G，但这不是你要做的。读者如何区分h和H？为什么有些变量是大写的，而另一些则是小写的？每个变量意味着什么？

性能和std::endl不能结合在一起。

将int变量与0.99进行比较是无稽之谈。