C++模拟牛顿力学（2D）

简介

如何用计算机来模拟真实世界呢？计算机最大的功能是计算，而物理学的种种公式就把现实世界中的物理规律以数学的语言描绘了出来，从而使我们可以通过计算大致模拟现实世界的物体运动。因此不难想到把物理学定律（这里用的是牛顿力学）用计算机语言描述出来，模拟现实世界。这个项目就用C++实现了此功能，为了方便，目前实现的是二维宇宙。
项目完成后，只要输入各项参数，我们就可以用它模拟现实中的运动了。下面是用该项目模拟某飞船围绕某恒星运动的过程（注意，运动轨迹并不是预设的动画，而是通过计算得出的）。

代码

完整代码如下：

namespace NewtonianPhysics
{
	using Scalar = long double;//标量
	const Scalar G = 6.67259e-11l;//万有引力常数
	namespace _2D
	{
		class Vector//矢量
		{
		public:
			static const Vector Zero;
			Scalar x;
			Scalar y;
			Vector() :x(0.l), y(0.l) {}
			Vector(Scalar inX, Scalar inY) :x(inX), y(inY) {}
			friend Vector operator+(const Vector& a, const Vector& b)
			{
				return Vector(a.x + b.x, a.y + b.y);
			}
			friend Vector operator-(const Vector& a, const Vector& b)
			{
				return Vector(a.x - b.x, a.y - b.y);
			}
			friend Vector operator*(const Vector& vec, Scalar n)
			{
				return Vector(vec.x * n, vec.y * n);
			}
			friend Vector operator*(Scalar n, const Vector& vec)
			{
				return Vector(vec.x * n, vec.y * n);
			}
			Vector operator-()const
			{
				return Vector(-x, -y);
			}
			Vector& operator*=(Scalar n)
			{
				x *= n;
				y *= n;
				return *this;
			}
			Vector& operator+=(const Vector& right)
			{
				x += right.x;
				y += right.y;
				return *this;
			}
			Vector& operator-=(const Vector& right)
			{
				x -= right.x;
				y -= right.y;
				return *this;
			}
			Scalar LengthSq() const
			{
				return (x * x + y * y);
			}
			Scalar Length() const
			{
				return (std::sqrt(LengthSq()));
			}
			void Normalize()
			{
				float length = Length();
				x /= length;
				y /= length;
			}
			static Vector Normalize(const Vector& vec)
			{
				Vector temp = vec;
				temp.Normalize();
				return temp;
			}
		};
		const Vector Vector::Zero;
		class MassPoint
		{
		public:
			Scalar m;//质量/kg
			Vector location;//当前位置
			Vector v;//速度/(m/s)
			Vector f;//受到的合力/(N)
			MassPoint(Scalar mass, Vector _location = { 0,0 }, Vector velocity = { 0,0 })
				:m(mass), location(_location), v(velocity) {}
			void Update(Scalar deltaTime)
			{
				//由F=ma得出加速度并更新速度
				v += f * (1 / m) * deltaTime;
				//由速度得出位移并更新坐标
				location += v * deltaTime;
				f = Vector::Zero;
			}
		};
		class World
		{
		public:
			std::vector<MassPoint> objects;
			void SetGravitation()//设置引力
			{
				for (size_t i = 0; i < objects.size(); ++i)
				{
					for (size_t j = i + 1; j < objects.size(); ++j)
					{
						//F=GMm/r²
						objects[i].f += _2D::Vector::Normalize(objects[j].location - objects[i].location) * 
							(G * objects[i].m * objects[j].m / (objects[i].location - objects[j].location).LengthSq());
						objects[j].f += -objects[i].f;
					}
				}
			}
			void Update(Scalar deltaTime)
			{
				for (MassPoint& m : objects)
					m.Update(deltaTime);
			}
		};
	}
}

标量和矢量

首先，我们定义了标量和矢量的类型。标量就是只有大小的量，如时间、质量等，而矢量是既有大小又有方向的量，如速度、加速度。标量用内置的long double类型即可。
矢量就是数学中的向量，在这里我们用Vector类来表示。根据数学定理，任何一个矢量都可以用平面直角坐标系中的一个坐标来表示。在Vector类里面，我们用两个标量来表示横坐标和纵坐标，并且定义了矢量的加法、减法、与实数的乘法。Length函数的作用是求矢量的长度，而LengthSq则是求矢量长度的平方。因为Length函数需要开方，速度比较慢，所以能用LengthSq就尽量用LengthSq。而Normalize函数的作用则是将矢量变成方向相同、长度为1的单位矢量。

增量时间

根据主流物理学理论，时间是连续的，然而，通过计算机模拟现实世界时，我们只能“一帧一帧”地计算。例如，如果我们每一帧代表0.1s，当前物体速度为10m/s，我们只能默认在这0.1s内物体的速度保持10m/s不变，尽管实际上速度可能变化了。这有点类似于微分，然而微分可以引入无穷小量，而计算机计算每一帧的时间不可能是无穷小，所以计算出来和实际结果一定有误差。并且，每帧代表的时间越短，误差越小。
在上面的例子中，0.1s就叫做增量时间（$\Delta t)$。

质点

牛顿力学的核心就是将物体抽象成质点（有质量无大小的点）进行研究。这里的MassPoint类就表示质点。在MassPoint类的成员变量中，m表示质量，v表示当前的速度，f表示受到的合力，location是当前的坐标。
Update函数根据物理定律对质点的状态进行更新，它的参数就是增量时间。在这个函数中，我们先通过牛顿第二定律求出加速度，然后乘以增量时间得出速度，再乘以时间便得出了位移。

这样，我们便实现了牛顿力学的基本功能。

World

在模拟程序中，所有在同一个“宇宙”中（即可以相互作用）的质点都被放到一个World中，这样就可以方便地管理它们了。此外，World还提供了一键计算引力的方法，利用万有引力公式计算出各个质点之间的引力。如果想忽略引力，不调用该函数即可。

实例

下面是太阳系中太阳、水星、金星、地球的运动模拟，每帧代表地球上的1天，每秒10帧。输入的数据全部是真实的，所以模拟结果和现实很相近（如地球公转一周在程序中大约37秒）。
当然，大家也可以调节这些数据，或者加入更多的元素，实现更加复杂的运动过程（例如三体运动）。

using namespace NewtonianPhysics;
using namespace _2D;
void Print(const World & w)
{
	cleardevice();
	for (const auto& p : w.objects)
	{
		setfillcolor(YELLOW);
		setlinecolor(RED);
		setlinestyle(PS_SOLID, 1);
		Scalar r = 2e11;
		auto x = (p.location.x + r) * 800 / (2 * r), y = (r - p.location.y) * 800 / (2 * r);
		fillcircle(x, y, 5);
	}
}
int main()
{
	initgraph(800, 800); 
	World w; 
	w.objects.push_back(MassPoint(1.9891e30, { 0,0 }, { 0,0 }));//太阳
	//行星均为远日点数据
	w.objects.push_back(MassPoint(3.3011e23, { 6.8817e10,0 }, { 0,3.9563e4 }));//水星
	w.objects.push_back(MassPoint(4.8675e24, { 1.082e11,0 }, { 0,3.5e4 }));//金星
	w.objects.push_back(MassPoint(5.97237e24, { 1.521e11,0 }, { 0,2.93e4 }));//地球
	while (1)
	{
		auto s = GetTickCount64();
		Print(w);
		w.SetGravitation();
		w.Update(86400);
		while (GetTickCount64() - s <= 100);
	}
	closegraph();
	return 0;
}

这一期就先讲到这里，如果喜欢别忘了点个赞！