微分方程（Blanchard Differential Equations 4th）中文版Section5.6

周期强迫非线性系统与混沌

到目前为止，本章仅讨论了自主非线性系统。另一方面，在第四章中，我们遇到了许多重要的二阶非自主方程。特别是，周期强迫的谐振子方程显式依赖于时间，并表现出相当有趣的行为。例如，在第4.2节中，我们发现周期性强迫影响了振荡解的振幅，甚至在共振强迫的情况下可能导致系统“爆炸”。

此外，在第4.5节中，我们研究了一个特定的具有周期强迫的非线性系统（塔科马峡湾桥），并发现它也表现出了令人惊讶的行为。在本节中，我们将继续研究周期性强迫的非线性系统。

如同洛伦兹系统一样，这些系统通常展示混沌行为。在这个上下文中，“混沌”一词意味着解展示了无限多种不同的定性行为。此外，这些不同的解紧密地聚集在一起，以至于初始条件的任何变化都会对解的长期行为产生剧烈影响。由于这些原因，找到混沌系统的解析解是不可能的，因为存在太多不同类型的解。同样，进行完整的定性分析也是不可能的。相反，我们提供一些对这些系统令人震惊行为的窥视，以便我们在其他系统中识别“混沌”。需要注意的是，这些系统是数学研究的活跃领域。如果我们的介绍似乎不完整，那是因为完整的故事尚未被发现。

周期性强迫的杜芬方程

我们从一个“完全理解”的非线性系统的例子开始，即我们可以绘制出其相图并描述每个解的行为。然后，我们在这个系统中添加一个周期性强迫项，并发展出用于深入了解结果的非自主系统的技术。

我们考虑的非线性系统是形式为

$$\frac{dy}{dt}=v$$

$$\frac{dv}{dt}=y-y^3$$

的杜芬方程。
这是一个哈密顿系统，其能量为 $H(y,v)=\frac{v^2}{2}-\frac{y^2}{2}+\frac{y^4}{4}$。因此，我们可以通过绘制 $H$ 的等高线来绘制这个系统的相图（见图5.49）。该系统有两个中心平衡点，位于 $(\pm 1,0)$，以及一个鞍点平衡点，位于 $(0,0)$。当 $t\to -\infty$ 时，趋向于鞍点的解在 $t\to +\infty$ 时也会趋向于它。所有其他解都是周期性的，要么绕着一个中心循环，要么绕着所有三个平衡点循环。

我们通过添加一个正弦强迫项来修改这个系统，得到

$$\frac{dy}{dt}=v$$

$$\frac{dv}{dt}=y-y^3+ϵ\sin t$$

注意，强迫的幅度由 $ϵ$ 给出，周期为 $2\pi$（参见习题以了解此系统作为二阶方程的解释）。
回归映射

正如我们在第4.1节中指出的那样，处理系统的通常工具——矢量场、方向场和相图——在非自治系统中并不适用，因为矢量场随时间变化。我们实际上需要一个三维的图像，包含 $y$、$v$ 和 $t$ 轴。我们将解视为在三维 $yvt$ 空间中沿曲线移动。尝试绘制这样的图像的结果如图5.50所示。这个图包含仅三条解曲线，但仍然相当难以可视化。此外，为了理解解的长期行为，我们必须将 $t$ 轴延伸得非常远。我们需要找到一种更好的方式来描绘解的行为。

为了表示这个系统的解，我们使用了由首位认识和理解混沌系统含义的数学家亨利·庞加莱（Henri Poincaré）提出的一个思想。目标是用二维图像替代图5.50中的三维图像。固定一个初始点 $y(0)=y_0$，$v(0)=v_0$。在图5.51中，这个点在平面 $t=0$ 上表示为一个点。如果我们绘制平面 $t=T$（在某个稍后时间 $T$）在 $yvt$ 空间中，我们可以询问在时间 $t=0$ 时从 $(y_0,v_0)$ 开始的解在哪里与这个平面相交。通过沿着解曲线走到时间 $t=T$，我们可以确定这个交点的位置（见图5.51）。

现在，我们利用强迫项是周期性的事实。取 $T=2\pi$。当 $t=0$ 时开始的解在平面 $t=2\pi$ 上交于某个点 $(y_1,v_1)=(y(2\pi ),v(2\pi ))$。跟踪这个解的另一种方法是扩展 $t$ 轴。然而，由于强迫是周期性的，微分方程也是周期性地具有周期 $2\pi$。因此，当 $t=2\pi$ 时从 $(y_1,v_1)$ 开始的曲线是从平面 $t=0$ 上的 $(y_1,v_1)$ 向前平移 $2\pi$ 时间单位的曲线。所以我们可以在不增加 $t$ 轴的情况下再跟踪解 $2\pi$ 时间单位。经过 $2\pi$ 更多的时间单位后，解再次与平面 $t=2\pi$ 相交，这次是在一个点 $(y_2,v_2)$，即 $(y(4\pi ),v(4\pi ))$。我们可以通过将点 $(y_2,v_2)$ 移到平面 $t=0$ 上，然后再跟踪解 $2\pi$ 时间单位来重复这个过程（见图5.52）。

当 $t$ 很大时，这个图像仍然非常复杂，因为从 $t=0$ 到 $t=2\pi$ 的解有很多条。但解的行为可以从一个更简单的图像中恢复。假设我们将这个图像旋转，使得 $t$ 轴直接进入纸面。这样我们只看到解穿过平面 $t=0$ 的点，即 $(y_0,v_0)=(y(0),v(0))$，$(y_1,v_1)=(y(2\pi ),v(2\pi ))$，$(y_2,v_2)=(y(4\pi ),v(4\pi ))$，……这一系列点足以提供解的长期行为的一个良好概念。

因此，我们将三维相空间中的解曲线图像替换为一个函数，该函数将平面 $t=0$ 上的给定初始值 $(y_0,v_0)$ 映射到它的“首次返回”点 $(y_1,v_1)$，也在平面 $t=0$ 上。这个函数称为回归映射或庞加莱回归映射；我们可以将平面 $t=0$ 上的点看作是通过回归映射“映射”到下一个与该平面交点的点。通过反复应用这个函数，我们可以看到解依次穿过平面的点。

无强迫系统的回归映射

为了说明这个过程，我们首先考虑一个简单的情况。令 $ϵ=0$，则我们有无强迫系统：

$$\frac{dy}{dt}=v$$

$$\frac{dv}{dt}=y-y^3$$

这是一个自治系统，但我们使用上述思想来描述回归映射在这种特殊情况下的图像。如果我们从一个接近平衡点 $(1,0)$ 的初始条件开始，那么在 $yvt$ 空间中，解会围绕这个平衡点螺旋运动（见图5.53）。这个解的回归映射选择在 $2\pi$ 时间间隔上的点，从而生成一系列围绕 $(1,0)$ 的点。

如果我们选择一个离 $(1,0)$ 较远的初始条件，解以不同的速率围绕 $(1,0)$ 螺旋运动，结果回归映射生成的点序列在 $(1,0)$ 周围进展得较慢。系统的平衡点给出的解总是返回到平面上的同一点。这些点被称为庞加莱回归映射的固定点。最后，从鞍点平衡点开始的解会给出一系列点，通过庞加莱回归映射从原点出发并返回到原点（见图5.54）。

在图5.55中，我们展示了将回归映射应用于多个不同初始条件的结果。在这种情况下，庞加莱回归映射的图像看起来非常像相空间图。接近平衡点 $(\pm 1,0)$ 的初始条件产生的解会沿着闭合循环保持在 $(\pm 1,0)$ 周围，因此在庞加莱映射的图像中，这些解生成了一系列围绕 $(\pm 1,0)$ 的点。由于解绕原点螺旋运动的速率随着半径的增加而减小，因此这些循环上的点的分布在不同的半径处是不同的。连接鞍点平衡点自身的轨道产生了从原点出发并返回到原点的点序列。

就像在相空间中一样，观察这些计算出的图像是非常有启发性的。回归映射中点序列出现的顺序提供的信息远比静态图像要多。然而，我们仍然可以利用完成的庞加莱映射图像来对解进行结论性分析。

回归映射对于强迫非线性系统

我们现在回到强迫非线性系统：

$$\frac{dy}{dt}=v$$
$$\frac{dv}{dt}=y-y^3+ϵ\sin t$$

为了数值计算解，我们需要选择一个 $ϵ$ 的值。为了说明问题，我们取 $ϵ=0.06$；其他值的 $ϵ$ 可以在练习中考虑。我们的目标是使用庞加莱回归映射图像来解释这个系统解的行为。

我们固定一个接近平衡点 $(1,0)$ 的初始条件，计算在 $yvt$-空间中的解，然后得到庞加莱回归映射图像。得到的图像如图5.56所示。该解的 $y(t)$ 图形见图5.57。我们可以利用这个图像预测解的行为：至少在时间是 $2\pi$ 的倍数时，解的位置保持相对接近 $(1,0)$。
以原点附近的鞍点为起点的解则更为有趣。该解的庞加莱回归映射如图5.58所示。与将原点连接到自身的简单曲线不同，我们看到在原点两侧出现了一簇点。

图5.58 使我们能够预测，起始于原点附近的解表现得非常不稳定。$y$-坐标会同时取正值和负值。此外，由于这些点似乎没有遵循任何特定的模式，我们可以预测这个解的 $y(t)$ 图形会无规则地振荡。这正是我们在图5.59中观察到的情况。

我们可以通过重新查看未受迫系统的相位图（见图5.60）来理解为什么这个解会这样表现。一个从原点附近开始的解，例如在 $(0,0)$ 右侧，靠近离开原点的分界线的地方，会在远离原点一段时间后回到原点附近。由于强迫项非常小，因此在解远离 $(0,0)$ 时对解的影响很小，因为该区域的向量场幅度相对较大。因此，我们预计解会在靠近入射分界线时返回到这一附近。当解接近原点时，未受迫系统的向量场较小，因此强迫项变得更加重要。当解接近 $(0,0)$ 时，如果强迫项推动解向上（即，$0.06\sin t>0$），那么解会向 $v>0$ 区域移动，从而沿着离开分界线返回到右半平面。

另一方面，如果当解接近原点时强迫项为负，则可能会有足够的“推动力”使解越过入射分界线。在这种情况下，解会进入左半平面。一旦进入左半平面，解就会绕左侧平衡点做一个循环。因此，对于受迫系统的每一个循环，解在返回到原点附近时必须“做出选择”。解可以进入右半平面或左半平面。解的方向取决于解相对于原点的位置以及当解靠近原点时强迫项的符号。因此，这个选择依赖于时机。
因此，我们无法预测当初始条件略有变化时会发生什么。当解接近 $(0,0)$ 时移动较慢，所以一个稍微改变的初始条件将解推得更接近原点，会显著影响解在 $(0,0)$ 附近停留的时间。这反过来会影响解返回到原点附近的时间，因此可能会影响解进入 $y-v$ 平面的哪个侧面。因此，初始条件的微小差异会在长期行为上产生剧烈的影响。这一点在图5.61中得到了证明，该图展示了两个初始条件非常接近的解的 $y(t)$ 图。解在一段时间内保持接近，但在之后的某个时刻，它们之间的差距足够大，从而在接近原点时做出了不同的决策。此后，解变得截然不同。

现实检验

这种行为相当令人不安。我们知道系统是“确定性的”。解的行为完全由系统微分方程的右侧决定。然而，当我们观察这些解时，它们的行为似乎没有任何特定的模式。我们甚至给解赋予了拟人化的特征，比如说“解决定了前进的方向。”解并不会“思考”。它们不需要思考，因为它们的行为由微分方程的右侧决定。发生的情况是，初始条件的微小变化会导致系统长期行为的巨大变化。

虽然这种行为相比于我们在之前章节研究的微分方程看起来不寻常，但在自然界中并不罕见。物理系统，如湍流中的水流、地球上的天气模式，甚至抛硬币，都以这种方式表现。这些系统是确定性的，因为它们遵循严格的物理规律。这并不意味着它们是可预测的。初始条件的微小变化可以对它们的行为产生巨大的影响。

依赖这些系统的数值模拟也存在风险。我们知道每种数值方法都只提供解的近似值。在模拟的每一步中都会存在小的误差。对于像上述这样的系统，数值近似中的一个小误差会导致一个与预期解略有不同的近似解。但相邻的解可能具有截然不同的长期行为。因此，数值模拟可能会给出与期望解相差甚远的结果。

这也是为什么长期（超过五天）天气预报通常不太准确的原因之一。对天气系统的知识不完全以及数值模拟中的误差会导致预测结果可能远离实际情况。

周期性强迫的摆

作为另一个周期性强迫非线性系统的例子，我们回到描述摆运动的系统。我们可以把摆看作是放在桌子上并且周期性地被摇动的摆。事实证明，上述观察到的许多行为也会在这个系统中发生。

方程

质量为 1，臂长为 1 的周期性强迫摆的方程为：

$$\frac{d\theta }{dt}=v$$

$$\frac{dv}{dt}=-g\sin \theta +ϵ\sin t$$

其中 $g$ 是重力常数。强迫项 $ϵ\sin t$ 模拟了一个外部力，这个外部力周期性地使摆向顺时针和逆时针方向施加幅度为 $ϵ$，周期为 $2\pi$ 的力。为了方便起见，我们假设时间和距离的单位已经选择，使得 $g=1$，因此我们的系统为：

$$\frac{d\theta }{dt}=v$$

$$\frac{dv}{dt}=-\sin \theta +ϵ\sin t$$

回归映射

我们按照上述方式构造强迫摆系统的回归图。强迫项的周期再次为 $2\pi$，因此我们在 $\theta vt$ 空间中跟踪从 $t=0$ 平面开始的解，并标记它们穿越 $t=2\pi$ 平面的位置。

对于以下示例，我们固定 $ϵ=0.01$；其他 $ϵ$ 的值在练习中考虑。在图 5.62 中，我们展示了初始条件接近 $(0,0)$ 的解的庞加莱回归图。得到的厚实的回路对应于一个振幅变化的摆动解。图 5.63 中相同解的 $\theta (t)$ 图显示了这种振荡。

对于强迫摆系统，初始条件接近 $(0,0)$ 对应于以小角度和小速度启动摆。在这种情况下，未强迫（无阻尼）摆以小的恒定振幅永远振荡。强迫项的加入意味着，就像强迫谐振子一样，强迫有时会沿着运动方向施加力，使摆摆动得更高，有时则会反向施加力，使摆摆动得更低。与谐振子不同的是，摆动的周期依赖于振幅。因此，强迫以一种长期上非常复杂的方式向系统添加和减去能量。（从图像中这不是很明显，因为强迫很小。）

靠近鞍点的解

图 5.64 显示了周期性强迫摆方程（$ϵ=0.01$）的庞加莱回归图，初始条件接近 $(-\pi ,0)$。这些点形成了一个“云”，没有特定的结构。此外，$\theta$ 坐标变得非常大。这意味着摆的摆臂在一个方向上已经完全旋转了多次。

就像在强迫杜芬系统中一样，每当解接近一个鞍点时，它必须“决定”转向哪个方向。如果解保持在上半平面中，则解的 $\theta$ 坐标在返回到另一个鞍点附近之前增加了 $2\pi$ 的倍数。如果它选择进入下半平面，则 $\theta$ 坐标在做出另一个选择之前减少了 $2\pi$。如果我们在 $t\theta$ 平面上绘制解的 $\theta$ 坐标，我们会看到它的运动方式非常混乱。

特别是，如果摆反复“决定”以相同的方向旋转，解的 $\theta$ 坐标可能会变得非常大（正或负）（参见练习 5–8）。

解在接近鞍点时的转向决策依赖于强迫项的符号。此时，解的行为对初始条件非常敏感。如果我们从接近 $(-\pi ,0)$ 的两个几乎相同的初始条件开始解，它们最终会分开并变得截然不同（参见图 5.65）。

从上述图像中，我们可以推断出强迫摆系统中存在一些有趣的行为。初始条件接近 $(-\pi ,0)$ 的强迫摆系统的解对应于几乎垂直但速度非常小的初始摆放位置。当然，摆会向下摆动。在摆动过程中，强迫项对摆动的影响非常小。当摆摆到几乎垂直的位置时，它再次减速，强迫项的影响变得更加明显。当摆接近摆动的最高点时，强迫项要么“推过顶部”，使其在相同方向上再转一圈，要么“拉回”，使其按原路返回。摆的运动方向取决于强迫项的符号，这又依赖于摆到达摆动顶部的时间。由于摆在摆动顶部附近的速度非常慢，因此初始条件的微小变化可以大大改变时机，从而改变摆的方向。

我们强调，上述物理论证并不是要替代庞加莱回归图及解的分析。尽管物理论证有其合理性，但它并不能说明 $ϵ=0.01$ 的强迫项是否足够大以导致这种行为。

总结

本节的总结是，即使是我们“理解”的系统，如摆，当添加周期性强迫等附加项时，也可以变得非常复杂。周期性强迫摆系统初看起来并不那么复杂，但从庞加莱回归图中我们看到其解表现得非常不可预测。初始位置的微小变化往往对解的行为产生剧烈影响。

如果这种我们现在称之为混沌的行为可以在像周期性强迫摆这样简单的系统中观察到，那么在自然界中发现它也就不足为奇了。这并不是一个新发现。亨利·庞加莱在 1880 年代首次考虑了自然界中存在“混沌”的可能性。他当时研究了一颗小行星（小行星）在恒星（太阳）和一颗大型行星（木星）的引力影响下的运动。庞加莱开发了回归图来研究该系统的行为。令人惊讶的是，庞加莱在那时没有像我们今天那样使用数值模拟结果的优势。尽管如此，他可以看到回归图会以非常复杂的方式表现出来，他明智地表示自己不会尝试绘制它。