数学中最重要的难题之一，求解微分方程，人类支配自然的第一步

欢迎大家来到IT世界,在知识的湖畔探索吧!

自然的语言究竟是什么？

思考一下这个问题。如果你问一个物理学家，我敢打赌他们的回答会是微分方程（differential equation）的某种变体。自然界中的一切，从原子的振动，到飓风的形成，再到行星和恒星在浩瀚太空中的运动，都是由微分方程控制的。

实际上，流体的混沌流动是由纳维-斯托克斯方程（Navier-Stokes）微分方程支配的，甚至牛顿的第二定律F=ma也是一个微分方程，因为加速度是位置的二阶导数。

这个想法宏大而美丽，我们可以用一个相对简单的数学概念来描述周围的许多事物。然而，有一个不小的问题：大多数微分方程根本无法求解。自然界的所有秘密都被它们的数学复杂性完全封锁住了，或者说真的如此吗？虽然大多数方程不可解，但可以通过迭代找到近似解。这个主题非常广泛且复杂。不过现在我们先聚焦一个简单的例子。

一个在弹簧上来回运动的方块，这是一个经典的简谐振子（Harmonic Oscillator）例子。事实上，它足够简单，以至于描述该系统的微分方程确实有一个确切的解。这其实就是牛顿的第二定律——力等于质量乘以加速度，再结合弹簧力定律——力等于负的弹簧常数乘以位移。

通过重新排列，可以解出加速度。当解这个微分方程时，会得到一个关于方块位移的正弦波解，随着时间变化而变化。

由于位移是正弦波形式的，所以速度也是。如果我们知道所有的参数，比如振幅、相位因子（即启动时钟的位置）、方块的质量和弹簧的刚度，那么就可以精确预测在任何时刻方块的位置。

不过假设你不知道如何解牛顿方程，那你该如何近似模拟这个系统的运动呢？通过使用牛顿第二定律和胡克定律，我们发现系统的加速度为负的弹簧常数除以质量再乘以位移。

如果将连续时间分割成离散的时间步长ΔT，就可以通过迭代计算速度v和位移x。

在每一个时间步长上，力会增加或减少方块的速度，

速度也会相应地增量位置,

每个时间步长，位移、速度和力都会发生变化，因此每次都需要重新计算。

让我们用这个方法来近似位移与时间的正弦曲线，

红线是精确解。你可以看到一开始近似还是挺好的，但随着时间推移，准确性迅速降低。这是任何近似方法的常见问题。然而还有其他的问题。注意到近似的振幅随着时间增加在增大，尽管系统是封闭的，这意味着弹簧在每个周期都会越来越伸展，显然这是不合逻辑的，而且也违反了能量守恒。

这其实是我们所用方法的一个特定缺陷——它不太能很好地守恒量，比如能量。经过多次迭代后，近似结果可能会趋向无限大。

通过选择较小的时间步长，可以缓解这个问题。步长越小，误差积累得越慢。不过，这也有问题，较小的步长会增加计算量，因为需要计算更多的迭代步骤。另一方面，能有一种能处理较大时间步长的技术会更好，这样既提高计算效率，也能在不需要过多小步长的情况下推进系统的演算。这种近似方法有其局限性，但通过调整时间步长和方法，我们可以逐步改善结果，并更好地模拟自然现象。

能否找到一种没有当前方法问题的方式？让我们用更通用的形式写出隐式欧拉方法。

首先，我们实际上是在解一个由两个微分方程组成的系统，但为了简化，现在我们只讨论一个一般性的方程，我们试图求解的方程形式是

带有某些初始条件

之前使用的方法实际上有个名字，叫欧拉法 (Euler’s Method)。由于有初始条件，我们通常称这种问题为初值问题 (Initial Value Problems，简称 IVP)。

再一次，如果我们解析地求解这个方程，那么就将时间t离散化为由小时间步长 Δt的集合。然后，通过欧拉法，如果知道某些初始条件x_i在时间t_i，那么通过微分方程，可以找到x_(i+1)的近似值

之所以能这么做，是因为根据原方程，函数 f(x,t)是该点处的导数，

因此我们本质上是在做一个微小的线性近似。

如前所述，时间步长足够小的情况下，这种方法有效。然而如果时间步长过大，会遇到问题。我们可以对欧拉法做一些小的调整来解决这些问题。通过隐式欧拉法，我们在计算x_(i+1)时，利用的是

而不是

通过这种方式，或许可以减少之前遇到的问题，因为我们在计算中考虑了下一个位置。

注意，这种数值方法将更加数值稳定。我们称这种方法为隐式的，因为x_(i+1)依赖于自身。而之前介绍的欧拉法可以被视为显式方法。

GIF中从上到下显示的是：精确解、显式欧拉法和隐式欧拉法。为什么隐式弹簧慢慢停止了？因为我们在计算中没有包含阻尼项，这就是我们做出的取舍。隐式方法比显式方法数值上更稳定，但换来的代价是每个时间步长会漏出一些能量，导致系统表现得像是有阻尼一样。实际上，如果我们重新排列隐式方法，就会发现它与显式方法的计算是反过来的。我们从x_(i+1)计算x_i，因此在某种意义上，隐式方法就是一个反向的欧拉法。从这个角度理解，就很容易明白为什么会有阻尼的现象。

显式欧拉法随着时间推移振幅会增大，而隐式方法则相反，振幅随着时间减少。因此，这种近似方法并不适合弹簧模拟。

那么我们如何结合显式和隐式方法的优点，同时弱化它们的缺点呢？最简单和最直接的方式可能是取它们的平均值，虽然不完全是，但想法很类似。我们这么做的目的是希望显式方法的能量增加可以抵消隐式方法的能量减少。为此，我们计算两个值K_1和K_2，

注意我们在计算K_2时使用的是显式欧拉法中的结果。然后我们将这两个值的平均值加到x_i 上，这就得到了更精确的结果。

这一方法能在保持数值稳定性的同时，也较好地守恒能量。

这就是我们的方法了，它比之前的方法更加数值稳定，同时能量守恒也更好。看一看，几乎完全吻合，

不过我们可以看出这仍然是一个近似解，因为随着时间的推移，准确度还是会下降，但这无疑是目前最好的方法。实际上，我们刚刚使用的这个技巧有一个专门的名称，叫做二阶龙格-库塔方法 (Runge-Kutta 2nd Order，RK2)。这个方法被称为“二阶”是因为我们平均了两个项K_1和K_2。事实上，显式欧拉法等同于一阶龙格-库塔方法，这就是为什么显式欧拉法中的公式也会出现在我们对 RK2 的计算中。

同样需要注意的是，龙格-库塔方法既可以是显式的，也可以是隐式的。现在，实际应用中最常用的是四阶龙格-库塔方法，简称 RK4，它也是基于相同的概念，只是加了几个额外的权重或 K 值。

在这里，我们不仅关心t_i和 t(i+1)​，还关心它们之间的中点。和之前一样，我们计算K_1，并将其加到x_i上，然后将其外推到中间点。在该点计算K_2​，并再次用该斜率找到另一个中间点，接着计算K_3​。最后，用K_3斜率完成一个完整的时间步长，并在该点计算K_4。

然后平均这四个斜率值，再乘以 Δt，得到 x(i+1)。

这里是完整的代数公式，

注意K_2和K_3的权重比K_1和K_4更高。这种权重分配来源于 RK4 的推导过程，其中涉及到泰勒级数。我更多的是想建立一种直观理解，说明为什么以及如何使用龙格-库塔方法。

这就是我们最终构建的近似器，RK4 广泛应用于科学和工程领域的初值问题。其他龙格-库塔方法也被广泛使用，例如Matlab的ODE45 就使用了四阶和五阶龙格-库塔方法的组合。你会注意到，虽然 RK2 方法已经足够好了，但 RK4 方法更加精确，随着时间推移，两者的准确度差距会越来越明显。具体使用哪种近似方法取决于应用场景，对于某些系统，显式或隐式欧拉法就足够了，而 RK4 通常是首选，因为它收敛速度快且总体上更稳定。

求解初值问题的近似解是数值分析领域中一个极为广泛的课题。但仅使用我概述的这些原理，你已经可以创造许多很酷的东西了。