65.9K
CodeProject 正在变化。 阅读更多。
Home

我的上帝，它充满了星星：使用 STL 实现隐式欧拉法

coding-thoughts
emptyStarIconemptyStarIconemptyStarIconemptyStarIconemptyStarIcon

0/5 (0投票)

2019 年 3 月 14 日

CPOL

6分钟阅读

viewsIcon

3955

使用 STL 实现隐式欧拉方法。

欢迎回到 n 体问题系列博文的又一精彩部分。这一次，我们将实现欧拉法，以使我们在上一篇博文中定义的测试通过。如果您是第一次看，或者需要刷新一下记忆，请参阅下面的先前博文链接。我还偶然发现了一篇非常好的关于 n 体问题主题的演示文稿。即使它相当数学化，它也能为我们提供良好的背景信息。

  1. 我的上帝，它充满了星星：n 体问题系列
  2. 我的上帝，它充满了星星：然后有了 CMake
  3. 我的上帝，它充满了星星：使用 Catch2 测试质点引力

在上一篇博文中，我们实现了两个简单的测试，但在实现求解器时，我意识到我们显然需要更多的测试。这不仅是为了测试求解器，我们还需要某种形式的 Vector2D（稍后，我们可能需要一个 3D 版本），因为我们必须存储和操作质点的位置、速度和加速度。所有这些属性都可以表示为二维空间中的二维向量。为了有用，Vector2D 必须提供一些典型的操作，例如与标量的加法、减法以及乘法和除法。此外，我们需要执行一些否定测试来比较 Vector2D。我们还需要测试除以零的情况，这由 operator/() 断言，但据我所知，在 2.6.1 版本中，Catch2 无法处理这种情况，而我们正在使用该版本。Vector2D 的测试如下所示。

TEST_CASE( "Test Vector2D operators and functions", "[vector]" ) {
    const Vector2D vecA = { 1, 1 };
    const Vector2D vecB = { 3, 3 };

    SECTION( "Comparing vectors" ) {
        const Vector2D expected = vecA;
        REQUIRE(vecA == expected);
    }

    SECTION( "Length of a vector" ) {
        REQUIRE(vecB.length() == Approx(sqrt(3*3 + 3*3)));
    }

    SECTION( "Adding two vectors" ) {
        const Vector2D expected = { 4, 4 };
        CHECK((vecA + vecB) == expected);

        Vector2D result = vecA;
        result += vecB;
        REQUIRE(result == expected);
    }

    SECTION( "Subtracting two vectors" ) {
        const Vector2D expected = { 2, 2 };
        CHECK((vecB - vecA) == expected);

        Vector2D result = vecB;
        result -= vecA;
        REQUIRE(result == expected);
    }

    SECTION( "Multiplying vector with scalar" ) {
        const Vector2D expected = { 2, 2 };
        CHECK((vecA * 2) == expected);

        Vector2D result = vecA;
        result *= 2;
        REQUIRE(result == expected);
    }

    SECTION( "Dividing vector with scalar" ) {
        const Vector2D expected = { 2, 2 };
        CHECK((vecA * 2) == expected);

        Vector2D result = vecA;
        result *= 2;
        REQUIRE(result == expected);
    }
}

TEST_CASE( "Negative test Vector2D operators and functions", "[vector]" ) {
    const Vector2D vecA = { 1, 1 };

    SECTION( "Comparing vectors" ) {
        Vector2D error = { 0, 1 };
        CHECK_FALSE(vecA == error);

        error = { 1, 0};
        CHECK_FALSE(vecA == error);

        error = { 0.999, 0.999 };
        REQUIRE_FALSE(vecA == error);
    }

    SECTION( "Dividing vector with scalar" ) {
        //Vector2D failing = vecA / 0;
        //Catching asserts not possible with CATCH2 https://github.com/catchorg/Catch2/issues/853
    }
}

有了 Vector2D 类，我们就满足了我们的测试并使其通过。运算符重载用于提供数学二维向量和标量之间的算术运算。重要的是要注意，相等性 operator==() 是通过比较两个向量的绝对差值来完成的。如果它们相等，则差值将接近零，这是通过将结果与 std::numeric_limits<double>::min() 进行比较来检查的。这个实用模板代表了特定类型的最小可能有限值，并且是 limits 头文件的一部分。另一种选择是检查两个值之间的关系，但目前这会导致太多后续问题，例如除以零。Bruce Dawson 建议在与零进行比较时，与基于绝对 epsilon 的值进行比较，我们就是这样使用 std::numeric_limits<double>::min() 的。

struct Vector2D {
    double x;
    double y;

    Vector2D() : x(0), y(0) {}
    Vector2D(double x, double y) : x(x), y(y) {}

    bool operator==(const Vector2D &rhs) const;
    bool operator!=(const Vector2D &rhs) const;

    double length() const;

    Vector2D& operator-=(const Vector2D& rhs);
    Vector2D& operator+=(const Vector2D& rhs);
    Vector2D& operator*=(const double& rhs);
    Vector2D& operator/=(const double& rhs);
};

Vector2D operator-(const Vector2D &lhs, const Vector2D &rhs);
Vector2D operator+(const Vector2D &lhs, const Vector2D &rhs);
Vector2D operator*(const Vector2D &lhs, const double &rhs);
Vector2D operator*(const double &lhs, const Vector2D &rhs);
Vector2D operator/(const Vector2D &lhs, const double &rhs);
bool Vector2D::operator==(const Vector2D &rhs) const {
    auto equalZero = std::numeric_limits<double>::min();

    return fabs(x - rhs.x) <= equalZero &&
           fabs(y - rhs.y) <= equalZero;
}

bool Vector2D::operator!=(const Vector2D &rhs) const {
    return !(rhs == *this);
}

double Vector2D::length() const {
    return sqrt(x*x + y*y);
}

Vector2D& Vector2D::operator-=(const Vector2D& rhs) {
    *this = *this - rhs;
    return *this;
}

Vector2D& Vector2D::operator*=(const double& rhs) {
    *this = *this * rhs;
    return *this;
}

Vector2D& Vector2D::operator/=(const double& rhs) {
    *this = *this / rhs;
    return *this;
}

Vector2D &Vector2D::operator+=(const Vector2D &rhs) {
    *this = *this + rhs;
    return *this;
}

Vector2D operator-(const Vector2D &lhs, const Vector2D &rhs) {
    return { lhs.x - rhs.x,
             lhs.y - rhs.y };
}

Vector2D operator+(const Vector2D &lhs, const Vector2D &rhs) {
    return { lhs.x + rhs.x,
             lhs.y + rhs.y };
}

Vector2D operator*(const Vector2D &lhs, const double &rhs) {
    return { lhs.x * rhs,
             lhs.y * rhs };
}

Vector2D operator*(const double &lhs, const Vector2D &rhs) {
    return rhs * lhs;
}

Vector2D operator/(const Vector2D &lhs, const double &rhs) {
    assert(fabs(rhs) > 0);

    return { lhs.x / rhs,
             lhs.y / rhs };
}

我们还需要扩展求解器测试，现在这些测试也运行一些性能测试。我们这样做是因为我们想大致确定我们的算法有多有效。此外，这使我们能够稍后比较不同的算法。根据 Christoph Schäfer 的论文，我们期望计算复杂度大约为 O(N^2) 。即使 PERFORMANCE 子句处于 beta 状态，它也工作得相当好。

double Inverse(double value) { return -value; }

Particle GenerateStandardParticle(double xPosition, double yPosition) {
    ParticleBuilder particleBuilder;

    return particleBuilder
            .position({xPosition, yPosition})
            .mass(1e10)
            .build();
}

TEST_CASE( "Explicit euler algorithm with two point mass", "[euler]" ) {
    const double EPSILON = 0.1;
    Solver solver(EPSILON);

    const std::vector<Particle> particlesX = { GenerateStandardParticle(0.5, 0),
                                               GenerateStandardParticle(-0.5, 0)};

    const std::vector<Particle> particlesY = { GenerateStandardParticle(0, 0.5),
                                               GenerateStandardParticle(0, -0.5)};

    //Solution
    const double acceleration = -0.6674079993; //m/s^2
    const double velocity = -0.06674079993; //m/s
    const double position = 0.48665184; //m

    SECTION( "Two still standing point mass are attracting each other in x-direction" ) {
        std::vector<Particle> result = solver.solve(particlesX);

        Particle &particle = result.front();
        CHECK(particle.getAcceleration().x == Approx(acceleration));
        CHECK(particle.getVelocity().x == Approx(velocity));
        CHECK(particle.getPosition().x == Approx(position));

        particle = result.back();
        CHECK(particle.getAcceleration().x == Approx(Inverse(acceleration)));
        CHECK(particle.getVelocity().x == Approx(Inverse(velocity)));
        REQUIRE(particle.getPosition().x == Approx(Inverse(position)));
    }

    SECTION( "Two still standing point mass are attracting each other in y-direction" ) {
        std::vector<Particle> result = solver.solve(particlesY);

        Particle &particle = result.front();
        CHECK(particle.getAcceleration().y == Approx(acceleration));
        CHECK(particle.getVelocity().y == Approx(velocity));
        CHECK(particle.getPosition().y == Approx(position));

        particle = result.back();
        CHECK(particle.getAcceleration().y == Approx(Inverse(acceleration)));
        CHECK(particle.getVelocity().y == Approx(Inverse(velocity)));
        REQUIRE(particle.getPosition().y == Approx(Inverse(position)));
    }
}

TEST_CASE("Benchmarking euler", "[benchmark]") {
    const double EPSILON = 0.1;
    Solver solver(EPSILON);

    ParticleBuilder particleBuilder;
    std::vector<Particle> particles = particleBuilder.build(100);
    BENCHMARK("Benchmarking with 100 particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(200);
    BENCHMARK("Benchmarking with 200 particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(400);
    BENCHMARK("Benchmarking with 400 particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(800);
    BENCHMARK("Benchmarking with 800 particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(1600);
    BENCHMARK("Benchmarking with 1.6K particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(3200);
    BENCHMARK("Benchmarking with 3.2K particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(6400);
    BENCHMARK("Benchmarking with 6.4K particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(12800);
    BENCHMARK("Benchmarking with 12.8K particles") {
        particles = solver.solve(particles);
    }

    particles = particleBuilder.build(25600);
    BENCHMARK("Benchmarking with 25.6K particles") {
        particles = solver.solve(particles);
    }
}

因为我们需要生成任意数量的质点，所以我们还引入了一个简单的 ParticleBuilder 工具类，它基于改编和简化的建造者模式，并带有流畅接口，可以生成单个粒子或任意数量的粒子。在我们的例子中，我们使用它来生成我们的基准粒子。建造者模式的基本原理是，每个配置质点特定属性的方法都会返回对自身的引用。这样，我们可以将所有（或仅部分，如果需要）方法连接在一起，最后构建质点。质点的实际生成是在定义过程的最后完成的，这被称为惰性初始化。为了生成质点属性的随机值，我们按照 Vorbrodt 关于随机数生成器的文章的建议，使用了 std::random_device，并将其作为 std::mersenne_twister_engine::mt19937 的随机种子来生成我们需要的实际随机数。随机数的范围将在 1 到要生成的质点数量之间均匀分布。

class ParticleBuilder {
public:
    ParticleBuilder() : mMass(0) {}

    ParticleBuilder & position(const Vector2D &position);
    ParticleBuilder & velocity(const Vector2D &velocity);
    ParticleBuilder & acceleration(const Vector2D &acceleration);
    ParticleBuilder & mass(double mass);
    Particle build() const;
    std::vector<Particle> build(size_t numberOfParticles);

private:
    double mMass;
    Vector2D mAcceleration;
    Vector2D mVelocity;
    Vector2D mPosition;
};
ParticleBuilder & ParticleBuilder::position(const Vector2D &position)
{
    mPosition = position;
    return *this;
}

ParticleBuilder & ParticleBuilder::velocity(const Vector2D &velocity)
{
    mVelocity = velocity;
    return *this;
}

ParticleBuilder & ParticleBuilder::acceleration(const Vector2D &acceleration)
{
    mAcceleration = acceleration;
    return *this;
}

ParticleBuilder & ParticleBuilder::mass(double mass)
{
    mMass = mass;
    return *this;
}

Particle ParticleBuilder::build() const
{
    return {mMass, mAcceleration, mVelocity, mPosition};
}

std::vector<Particle> ParticleBuilder::build(size_t numberOfParticles)
{
    std::vector<Particle> particle(numberOfParticles);

    std::mt19937 mt(std::random_device{}());
    std::uniform_real_distribution real_dist(1.0, static_cast<double>(numberOfParticles));

    for(Particle &p : particle) {
        p = mass(real_dist(mt))
                .acceleration({ real_dist(mt), real_dist(mt) })
                .velocity({ real_dist(mt), real_dist(mt) })
                .position({ real_dist(mt), real_dist(mt) })
                .build();
    }

    return particle;
}

Particle 类然后基本上是直观的。唯一需要我们以某种方式定义的重要主题是，我们如何确定两个粒子是否相同。这不能仅仅通过比较它们的引用来完成，因为我们从不修改粒子本身，而是每次操作后都生成新的粒子。所有粒子都被视为不可变的。因此，我们需要以某种方式为每个质点分配一个唯一的 ID，我们通过简单地将一个类 IDCounter 赋给 id 成员并之后递增 IDCounter 来实现这一点。

class Particle {
public:
    Particle();
    Particle(double mass, const Vector2D &acceleration, 
             const Vector2D &velocity, const Vector2D &position);

    bool operator==(const Particle &rhs) const;
    bool operator!=(const Particle &rhs) const;

    double getMass() const;

    const Vector2D &getAcceleration() const;
    void setAcceleration(const Vector2D &acceleration);

    const Vector2D &getVelocity() const;
    void setVelocity(const Vector2D &velocity);

    const Vector2D &getPosition() const;
    void setPosition(const Vector2D &position);

private:
    static size_t IDCounter;

    size_t id;
    double mass;
    Vector2D acceleration;
    Vector2D velocity;
    Vector2D position;
};
size_t Particle::IDCounter = 0;

Particle::Particle()
    : mass(0) {
    id = IDCounter++;
}

Particle::Particle(double mass, const Vector2D &acceleration, 
                   const Vector2D &velocity, const Vector2D &position)
        : mass(mass), acceleration(acceleration), velocity(velocity), position(position) {
    assert(mass > 0);
    id = IDCounter++;
}

bool Particle::operator==(const Particle &rhs) const {
    return id == rhs.id;
}

bool Particle::operator!=(const Particle &rhs) const {
    return !(rhs == *this);
}

const Vector2D &Particle::getAcceleration() const {
    return acceleration;
}

void Particle::setAcceleration(const Vector2D &acceleration) {
    Particle::acceleration = acceleration;
}

const Vector2D &Particle::getVelocity() const {
    return velocity;
}

void Particle::setVelocity(const Vector2D &velocity) {
    Particle::velocity = velocity;
}

const Vector2D &Particle::getPosition() const {
    return position;
}

void Particle::setPosition(const Vector2D &position) {
    Particle::position = position;
}

double Particle::getMass() const {
    return mass;
}

最后但同样重要的是，我们有了欧拉法的实际实现。Solver 类通过必要的时间步长 \epsilon 进行初始化，并且只有一个简单的接口方法，该方法接受一个 std::vector<Particle> 作为参数，并在计算完成后返回结果。为了执行其计算，solve 方法使用了 calculateAcceleration、calculateVelocity 和 calculatePosition 方法，它们基本上都具有相同的接口。计算的本质，即一个质点由所有其他质点的引力产生的加速度的计算，然后在 static 函数 AccumulateAcceleration 中完成。您可能会注意到我们强烈专注于使用 STL 提供的功能，甚至没有真正的 for 循环。之所以这样做，是因为我有一个想法，稍后使用 STL 的执行策略，但随后意识到 libc++ (P0024R2) 目前不支持这一点。

class Solver {
public:
    explicit Solver(double mEpsilon);

    std::vector<Particle> solve(const std::vector<Particle> &particles) const;

private:
    std::vector<Particle> calculateAcceleration(const std::vector<Particle> &particles) const;
    std::vector<Particle> calculateVelocity(const std::vector<Particle> &particles) const;
    std::vector<Particle> calculatePosition(const std::vector<Particle> &particles) const;
    static Particle AccumulateAcceleration
            (const std::vector<Particle> &particles, const Particle &particle);
    static double CalculateEquivalentMass(const Particle &particleA, const Particle &particleB);

    static const double G;
    static const double EPSILON;
    double mEpsilon;
};
const double Solver::G = 6.67408e-11;
const double Solver::EPSILON = 1e-3;

Solver::Solver(double mEpsilon) : mEpsilon(mEpsilon) {}

std::vector<Particle> Solver::solve(const std::vector<Particle> &particles) const {
    std::vector<Particle> solution = calculateAcceleration(particles);
    solution = calculateVelocity(solution);
    solution = calculatePosition(solution);

    return solution;
}

std::vector<Particle> Solver::calculateAcceleration(const std::vector<Particle> &particles) const {
    std::vector<Particle> solution(particles.size());

    std::transform(begin(particles), end(particles), begin(solution), 
                  [&particles](const Particle &particle) {
        return AccumulateAcceleration(particles, particle);
    });

    return solution;
}

std::vector<Particle> Solver::calculateVelocity(const std::vector<Particle> &particles) const {
    std::vector<Particle> solution(particles.size());

    std::transform(begin(particles), end(particles), begin(solution), 
                  [epsilon = mEpsilon](Particle particle) {
        const Vector2D v0 = particle.getVelocity();
        particle.setVelocity(v0 + particle.getAcceleration()*epsilon);

        return particle;
    });

    return solution;
}

std::vector<Particle> Solver::calculatePosition(const std::vector<Particle> &particles) const {
    std::vector<Particle> solution(particles.size());

    std::transform(begin(particles), end(particles), 
                   begin(solution), [epsilon = mEpsilon](Particle particle) {
        const Vector2D v = particle.getVelocity();
        const Vector2D r0 = particle.getPosition();
        particle.setPosition(r0 + v*epsilon + particle.getAcceleration()*epsilon*epsilon);

        return particle;
    });

    return solution;
}

Particle Solver::AccumulateAcceleration
      (const std::vector<Particle> &particles, const Particle &particle) {
    Particle particleA = particle;
    const double e3 = EPSILON*EPSILON*EPSILON;

    std::for_each(begin(particles), end(particles), [&particleA, e3](const Particle &particleB) {
        if(particleA != particleB) {
            const double M = CalculateEquivalentMass(particleA, particleB);
            const Vector2D r = particleB.getPosition() - particleA.getPosition();
            const double rLength = r.length();
            const double r3 = fabs(rLength*rLength*rLength);

            const Vector2D a0 = particleA.getAcceleration();
            particleA.setAcceleration(a0 + G*M*r/(r3 + e3));
        }
    });

    return particleA;
}

double Solver::CalculateEquivalentMass(const Particle &particleA, const Particle &particleB) {
    const double massA = particleA.getMass();
    assert(massA > 0);
    
    const double massB = particleB.getMass();

    return massA*massB/massA;
}

到目前为止，我们已经完成了相当多的工作。我们现在有一个简单但有效的算法，可以计算 n 体问题离散时间步长的解决方案。但我们最初的草案表现如何？首先，所有测试都通过了，太棒了。我们甚至能够计算相当数量的质点。但我们必须提醒自己，对于 25.6K 个质点，仅一个时间步就需要大约 45 秒。这相当长。如果我们只看数字，可能会意识到欧拉法的计算复杂度声明是 O(N^2) ，下图更明显地显示了这一点。这是完全正确的。如果我们想将 calculateAcceleration 方法中的 std::transform 和 AccumulateAcceleration 函数中的 std::for_each 算法替换为简单的 for 循环，我们会很容易地意识到我们有一个 N 个质点的嵌套 for 循环。其余的只是数学，N \cdot N=N^2 。

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
solverTest is a Catch v2.6.1 host application.
Run with -? for options

-------------------------------------------------------------------------------
Benchmarking euler
-------------------------------------------------------------------------------
/mnt/c/Develop/gravity/solverTest/src/solverTest.cpp:63
...............................................................................

benchmark name                                  iters   elapsed ns      average 
-------------------------------------------------------------------------------
Benchmarking with 100 particles                     1       532000       532 µs 
Benchmarking with 200 particles                     1      3663900    3.6639 ms 
Benchmarking with 400 particles                     1     12501000    12.501 ms 
Benchmarking with 800 particles                     1     49732100   49.7321 ms 
Benchmarking with 1.6K particles                    1    182535600   182.536 ms 
Benchmarking with 3.2K particles                    1    711535000   711.535 ms 
Benchmarking with 6.4K particles                    1   2612469200    2.61247 s 
Benchmarking with 12.8K particles                   1  11103358400    11.1034 s 
Benchmarking with 25.6K particles                   1  46714864800    46.7149 s 

===============================================================================
All tests passed (25 assertions in 4 test cases)

eulerPlot

到目前为止，我们只关注了可读性，我认为我们还可以做得更多。但在下一篇博文中，我想重点关注两个特定问题：算法的效率和并行化。如果您想下载此状态下的项目，请随时获取v0.4.0。

您喜欢这篇博文吗？

您有什么想法吗？

请随时评论和分享这篇博文。

© . All rights reserved.