简介

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

代码

完整代码如下：

namespace NewtonianPhysics
{
	using Scalar = long double;//标量
	const Scalar G = 6.67259e-11l;//万有引力常数
	namespace _2D
	{
		class Vector//矢量
		{
		public:
			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*=(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;
			}
		};
		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;
			}
		};
	}
}

标量和矢量

首先，我们定义了标量和矢量的类型。标量就是只有大小的量，如时间、质量等，而矢量是既有大小又有方向的量，如速度、加速度。标量用内置的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函数根据物理定律对质点的状态进行更新，它的参数就是增量时间。在这个函数中，我们先通过牛顿第二定律求出加速度，然后乘以增量时间得出速度，再乘以时间便得出了位移。

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

实例

下面是简介中飞船围绕恒星运动的代码示例（绘图用了easyx）：

using namespace NewtonianPhysics;
void Print(const _2D::MassPoint& p)
{
	setfillcolor(YELLOW);
	setlinecolor(RED);
	setlinestyle(PS_SOLID, 1);
	auto x = (p.location.x + 2e7) * 800 / 4e7, y = (2e7 - p.location.y) * 800 / 4e7;
	fillcircle(x, y, 5);
}
int main()
{
	initgraph(800, 800); 
	_2D::MassPoint a(5e30l, { 0,0 }), b(1e4, { 1e7,0 }, { 0,5e6 });
	while (1)
	{
		cleardevice();
		Print(a);
		Print(b);
		Sleep(100);
		//F=GMm/r²
		b.f = _2D::Vector::Normalize(a.location - b.location) * (G * a.m * b.m / (a.location - b.location).LengthSq());		a.Update(0.05);
		b.Update(0.1);
	}
	closegraph();
	return 0;
}

该代码中，恒星的质量为 $5\times10^{30}kg$ ，且位于 $(0, 0)$ 静止不动；飞船的质量为 $1\times10^4kg$ ，初始位于 $(1\times10^7,0)$ ，并且具有向上 $5\times10^6m/s$ 的初速度。通过万有引力公式实时计算出恒星对飞船的引力，并调用Update函数进行更新。这样便实现了模拟运动。当然，大家也可以调节这些数据，或者加入更多的元素，实现更加复杂的运动过程（例如三体运动）。
这一期就先讲到这里，如果喜欢别忘了点个赞！