微分方程基础

1.初值解问题

• 微分方程描述了未知函数与其导数之间的关系。 求解微分方程就是找到满足关系的函数，通常同时满足一些附加条件。 在本课程中，我们将主要关注一类特定的问题，称为初始值问题。 在典型的初始值问题中，系统的行为由以下形式的常微分方程 (ODE) 描述
  

• f为已知函数，x 代表当前系统的状态， x˙ 为 x 对时间t的导数，通常，x 和 x˙ 是向量。顾名思义，对于一个初始值问题，给定开始时间 t0 ，初始状态 x(t0) = x0，并希望在之后的时间求得x。

• 通用的初始值问题很容易可视化。 在 2D 中，x(t) 扫出一条曲线，描述平面中点 p 的运动。 在任何点 x 处，函数 f 可以被计算以提供 2-vector，因此 f 定义了平面上的向量场（见图 1。） x 处的向量是移动点 p 在移动时必须具有的速度 通过 x（它可能会或可能不会。）将 f 视为从点到点驱动 p，就像洋流一样。 无论我们最初将 p 存放在哪里，那个时候的“当前”都会抓住它。 p 的携带位置取决于我们最初将其丢弃的位置，但一旦丢弃，所有未来的运动都由 f 决定。 p到f扫出的轨迹形成矢量场的积分曲线。 见图 2。

• 我们将 f 写为 x 和 t 的函数，但导数函数可能直接取决于时间，也可能不直接取决于时间。 如果是这样，那么不仅点 p 移动，向量场本身也移动，因此 p 的速度不仅取决于它在哪里，还取决于它何时到达那里。 在这种情况下，导数 x˙ 以两种方式取决于时间：首先，导数向量本身摆动，其次，点p，因为它在轨迹 x(t) 上移动，所以在不同时间看到不同的导数向量。 如果您保持粒子漂浮在起伏的矢量场中的图像，这种双重时间依赖性不应该导致混乱。

2.数值解

• 标准的导论微分方程课程侧重于符号解，其中未知函数的函数形式将被猜测。 例如，微分方程 x˙ = -kx，其中 x˙ 表示 x 的时间导数，由 x = e-kt 满足。
• 相反，我们将只关注数值解，其中我们从初始值 x(t0) 开始采取离散时间步长。 迈出一步，我们使用导数函数 f 来计算 x 的近似变化，Δx，在时间间隔 Δt 内，然后将 x 增加 Δx 以获得新值。 在计算数值解时，导函数 f 被视为一个黑匣子：我们提供 x 和 t 的数值，然后接收 x˙ 的数值。 数值方法通过在每个时间步执行一个或多个这些导数评估来操作。

2.1 欧拉法 Euler’s Method

• 最简单的数值方法称为欧拉法。 让我们对 x 的初始值由 x0 = x(t0) 表示，我们在稍后的时间 t0 + h 对 x 的估计由 x(t0 + h) 表示，其中 h 是一个步长参数。 欧拉方法通过在导数方向上迈出一步来简单地计算 x(t0 + h)，
  

• 您可以使用 2D 矢量场的图片来可视化 Euler 方法。 p 不是真正的积分曲线，而是遵循多边形路径，其每条边都是通过在开始时评估向量 f 并按 h 缩放来确定的。 见图 3。
  

• 欧拉的方法虽然简单，但并不准确。 考虑积分曲线是同心圆的二维函数 f 的情况。 一个由 f 支配的点 p 应该永远在它开始的任何一个圆上运行。 相反，随着欧拉的每一步，p 将沿直线移动到半径更大的圆，因此它的路径将遵循向外的螺旋。 缩小步长会减慢这种向外漂移的速度，但永远不会消除它。

• 此外，欧拉方法可能不稳定。考虑一个一维函数 f = -kx，它应该使点 p 以指数方式衰减到零。对于足够小的步长，我们得到了合理的行为，但是当 h > 1/k 时，我们有 |Δx| > |x|，所以解在零附近振荡。超过h = 2/k，振荡发散，系统爆炸。参见图 4。
  

• 最后，欧拉的方法甚至没有效率。大多数数值求解方法几乎将所有时间都花在执行导数评估上，因此每一步的计算成本取决于每一步的评估次数。虽然 Euler 的方法每一步只需要一次评估，但方法的真正效率取决于它允许您采取的步骤的大小（同时保持准确性和稳定性）以及每一步的成本。更复杂的方法，即使有些方法每步需要多达四到五次评估，也可以大大优于 Euler 方法，因为它们每一步的更高成本被它们允许的更大步长所抵消。要了解我们如何改进 Euler 方法，我们需要更仔细地查看该方法产生的错误。
  理解发生了什么的关键是泰勒级数：假设 x(t) 是平滑的，我们可以将其在步骤结束时的值表示为包含开始时的值和导数的无限和：
  

• 如你所见，我们通过截断级数得到欧拉更新公式，丢弃除右侧前两项以外的所有项。这意味着欧拉方法只有在除第一个之外的所有导数都为零时才是正确的，即如果 x(t) 是线性的。误差项，即欧拉步长与完整的、未截断的泰勒级数之间的差异，主要由前项 (h2/2)x¨(t0) 决定。因此，我们可以将误差描述为 O(h2)（读作“Order h squared”。）假设我们将步长减半；也就是说，我们采取大小为 h/2 的步长。尽管这只会产生大约四分之一的误差，但我们必须在任何给定的时间间隔内采取两倍的步数。这意味着我们在 t0 到 t1 的区间内累积的误差线性地取决于 h。理论上，使用欧拉方法，我们可以通过选择一个合适的小 h，在 t0 到 t1 的区间内以尽可能小的误差数值计算 x。在实践中，可能需要很多时间步长，具体取决于误差和函数f。

2.2中点法 The Midpoint Method

• 如果我们能够评估 x¨ 以及 x˙，我们可以通过在截断的泰勒级数中保留一个附加项来实现 Oh3) 精度而不是 O(h2)：
  

• 回想一下，时间导数 x˙ 由函数 f (x(t), t) 给出。为简单起见，我们将假设导函数 f 仅通过 x 间接依赖于时间，因此 x˙ = f (x(t))。然后链式法则给出
  

• 为了避免计算 f˙ ，这通常是复杂和昂贵的，我们可以只用 f 来近似二阶项，并将近似值代入方程 1，留下 O(h***3) 误差.为此，我们执行另一个泰勒展开，这次是 f 的函数，

• 我们首先通过选择将 x¨ 引入这个表达式

• 以便

• 其中 x0 = x(t0)。我们现在可以将两边都乘以 h（将 O(h2) 项变为 O(h*3)）并重新排列，得到

• 将右侧代入等式1给出更新公式

• 这个公式首先计算一个欧拉步，然后在步的中点执行二阶导数计算，使用中点计算来更新 x。因此命名为中点法。中点法在 O(h3) 范围内是正确的，但需要对 f 进行两次评估。有关该方法的图示，请参见图 5。
  

• 我们不必因 O(h3) 的错误而停止。通过多次评估 f，我们可以消除越来越高阶的导数。最流行的方法是 4 阶 Runge-Kutta 方法，每步误差为 O(h*****5)。 （中点法可以称为二阶龙格-库塔法。）四阶龙格-库塔法我们不会推导，但计算x(t0 + h)的公式如下：
  

2.3自适应步长 Adaptive Stepsizes

• 无论采用何种底层方法，一个主要问题在于确定一个好的步长。理想情况下，我们希望选择尽可能大的 h，但不要大到给我们带来不合理的误差，或者更糟的是，导致不稳定。如果我们选择一个固定的步长，我们只能在 x(t) 的“最差”部分允许的范围内尽可能快地进行。我们想做的是随着时间的推移而改变 h。只要我们可以使 h 变大而不会产生太多错误，我们就应该这样做。当必须减少 h 以避免过多的错误时，我们也想这样做。这就是自适应步长的想法：在求解 ODE 的过程中改变 h。

• 在这里，我们将介绍 Euler 方法的自适应步长。基本思路如下。假设我们有一个给定的步长 h，我们想知道我们可以考虑改变多少。
  假设我们计算 x(t0 + h) 的两个估计值。我们通过从 t0 到 t0 + h 采取大小为 h 的欧拉步长来计算估计值 xa。我们还通过从 t0 到 t0 + h 采取两个大小为 h/2 的欧拉步长来计算估计值 xb。 xa 和 xb 都与 x(t0 + h) 的真实值相差 O(h2)。这意味着 xa 和 xb 彼此相差 O(h2)。因此，我们可以写出当前误差 e 的度量是
  

• 这使我们可以方便地估计采用大小为 h 的欧拉步长时的误差。

• 假设我们愿意每步有多达 10-4 的误差，而当前的误差仅为 10-8。由于误差随着 h2 上升，我们可以将步长增加到
  

• 相反，如果我们当前有 10-3 的误差，并且只能容忍 10-4 的误差，我们将不得不将步长减小到
  

• 自适应步长是一种强烈推荐的技术

2.4 实现

• 我们要求解的 ODE 可能代表许多事物——例如，质量和弹簧的集合、一些刚体或可变形物体。我们希望实现 ODE 求解器和它们运行的​​模型，以将它们与另一个的内部细节隔离开来。这将使得可以轻松更改求解器，并使求解器代码可重用。幸运的是，这种模块化并不难实现，因为所有求解器都可以用一组小的、刻板的操作来表示。据推测，ODE 管理的对象系统将体现在某种结构中。该方法是编写对该结构进行操作的特定于类型的代码以执行标准操作，然后根据这些通用操作实现求解器。

• 从求解器的角度来看，它所运行的系统是一个黑盒函数 f (x, t)。求解器需要能够根据需要在 x 和 t 的任何值下评估 f，然后安装采取时间步长时更新 x 和 t。为了支持这些操作，代表被求解的 ODE 的对象必须能够处理来自求解器的这些请求：

  • 返回dim (x)。 由于 x 和 x˙ 可能是向量，求解器必须知道它们的长度，才能分配存储空间、执行向量算术运算等。

  • 获取/设置 x 和 t。 求解器必须能够在步骤结束时安装新值。此外，多步骤方法必须在执行导数评估的过程中将 x 和 t 设置为中间值。

  • 在当前 x 和 t 处评估 f。

• 在面向对象的语言中，这些操作自然会被实现为以特定于类型的方式处理的通用函数。 在非面向对象的语言中，通用函数可以通过在结构槽中安装指向特定类型函数的指针来伪造，或者只是通过将函数指针作为参数传递给求解器。 稍后我们将详细考虑如何为特定模型（例如粒子和弹簧系统）实现这些操作。