在物理模拟项目中，经常需要处理复杂的数学模型和算法。本文将介绍如何通过CMake设置项目，实施测试驱动开发(TDD)，并实现Euler方法来解决物理问题。
CMake是一个跨平台的自动化构建系统，它使用配置文件（称为CMakeLists.txt）来生成标准的构建文件。根据Reddit上的一些建议，对CMake配置进行了一些调整，但仍然不确定当前的项目结构是否能够持续到项目结束。如果错过了之前的帖子，可以在这里找到： 和 。
为了实施TDD，首先需要定义一些测试案例，然后才能实现求解器库。一个简单的场景是两个具有等质量的点，距离为d，没有初始速度和加速度。显然，可以用它进行数百种不同的复杂测试。但想从两个非常简单和具体的测试开始，点质量只在x和y方向上加速和移动。这样，至少可以检查基础假设、公式和实现是否正确。
还需要考虑引力方程的异常，用它来计算加速度。当两个质量达到奇点，换句话说，碰撞时会发生什么？如果两个点在空间中处于同一位置，那么它们之间的距离r为零。如果发生这种情况，加速度F将变得无限大，因此两个点的加速度也会无限大。两个点将相互飞向虚空。
基本上有两种方法可以解决这个问题。第一种解决方案是将两个质量合并为一个。每个黑洞的粉丝可能会同意这个解决方案，但有一个小问题，合并恒星会发生什么，完全不理解，因此对来说很难建模。显然，模型并不完美，有很多缺陷，但合并恒星的模型会使其复杂性爆炸。第二种解决方案相对简单，不允许质量碰撞。这可以通过向引力方程添加一个错误项来实现。此外，还需要添加一个符号项来获得正确的力的方向。
现在有几种选择来实现所谓的显式和隐式方法，例如Euler、Euler-Cromer、Leapfrog、Verlet等。将从最简单的，但不太准确的显式Euler方法开始。
让从一个名为“显式Euler算法与两点质量”的简单测试案例开始，在solverTest可执行文件中。
TEST_CASE(
"Explicit euler algorithm with two point mass",
"[euler]"
) {
Solver solver;
ParticleBuilder particleBuilder;
Particle particleA = particleBuilder
.position()
.mass()
.build();
Particle particleB = particleBuilder
.position()
.mass()
.build();
std::vector<Particle> particles = { particleA, particleB };
const double EPSILON = 1e-3;
SECTION(
"Two still standing point mass are attracting each other"
) {
std::vector<Particle> result = solver.solve(particles, EPSILON);
REQUIRE(false);
}
}
这个测试案例描述了到目前为止粗糙的设计，希望如何使用求解器。需要一个求解器来执行计算，一个粒子向量的粒子，以及一个粒子构建器，它提供了一种流畅的接口，可以轻松地为生成粒子。Euler算法不仅需要粒子来计算，还需要一个时间步长，正如从第一个帖子中记得的，将在源代码中称之为EPSILON（不要与需要解决的奇异性混淆）。可以通过下载v0.2.0来获取这个起始状态。
第一次测试将检查两个粒子在x方向测试中的以下结果：所有3个结果在纯x方向和y方向移动的点质量的情况下都是有效的。在实现了测试并使其编译之后，solverTest文件看起来像这样：
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]"
) {
Solver solver;
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)};
const double EPSILON = 1e-3;
//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, EPSILON);
Particle &particle = result.front();
REQUIRE(particle.acceleration.x == Approx(acceleration));
REQUIRE(particle.velocity.x == Approx(velocity));
REQUIRE(particle.position.x == Approx(position));
particle = result.back();
REQUIRE(particle.acceleration.x == Approx(Inverse(acceleration)));
REQUIRE(particle.velocity.x == Approx(Inverse(velocity)));
REQUIRE(particle.position.x == Approx(Inverse(position)));
}
SECTION(
"Two still standing point mass are attracting each other in y-direction"
) {
std::vector<Particle> result = solver.solve(particlesY, EPSILON);
Particle &particle = result.front();
REQUIRE(particle.acceleration.y == Approx(acceleration));
REQUIRE(particle.velocity.y == Approx(velocity));
REQUIRE(particle.position.y == Approx(position));
particle = result.back();
REQUIRE(particle.acceleration.y == Approx(Inverse(acceleration)));
REQUIRE(particle.velocity.y == Approx(Inverse(velocity)));
REQUIRE(particle.position.y == Approx(Inverse(position)));
}
}
让首先看看SECTION。执行两个测试，一个测试纯x方向上的加速度、速度和位置，另一个测试纯y方向上的加速度、速度和位置。两个测试都使用了Inverse()函数，该函数用于测试第二个粒子的预期结果。Catch2的Approx包装类重载了比较运算符，并提供了3种不同的配置方法。
ParticleBuilder本身提供了一个流畅的接口，可以轻松配置和生成点质量。对于如此简单的案例，使用带有流畅接口的构建器模式可能没有必要，但将在后面的更高级生成中需要它。
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;
private:
double mMass;
Vector2D mAcceleration;
Vector2D mVelocity;
Vector2D mPosition;
};
Vector2D是一个简单的数据结构，可能会在以后扩展更多功能。测试结果如预期那样失败。