一、说明

        关于混沌如何实现？能否用计算机模拟？本文从简单的物理道具：双臂摆的物理方程，引进混沌理念。进而进行复杂的自然状态中。本文只是研究题目的引出，日后如果需要进一步加深，不妨提供一个踮脚的石头。

二、钟摆物理方程

        任何物理专业的学生都研究过低级钟摆，可能比他们想要的要多。与弹簧上的质量一样，摆锤是简单谐波运动的典型例子。那里的关键词很简单。您可能知道，它们非常容易进行数学建模。钟摆的周期可以用下面的公式表示，其中 g 是重力，L 是钟摆的长度。

        角度位移更有趣，但仍然不太深。方程本身更复杂，但实际上它只是一个余弦波。下面的公式表示单个钟摆的角位移，其中 θ 是角位移，A 是最大角位移（在大多数情况下这只是释放角），ω 是角频率，t 是时间。

        说完这些，让我们来看看更有趣的，尽管实际上是不可能的（你会看到），数学和物理。

        很有可能，如果你使用过普通的钟摆，你至少对双摆有点熟悉，这让人想起臭名昭著的三体问题。三体问题代表了这样一种思想，即在所有物体的力相互影响的系统中，对三个物体的运动的解析解是不可能求解的，因为每个物体的力同时与其他物体相互作用。

        简而言之，双摆是无法分析求解的。就上面提到的摆方程而言，您真正无法求解的是 θ 方程，即两个摆锤的角位移。你也无法求解周期，但这主要是因为双摆并不完全有规律的周期。下面是我们讨论过的双摆系统的直观表示，其中枢轴点锚定在墙上，角度用 θ1 和 θ2 标记，悬挂质量标记为 m1 和 m2。

        重申一下——在我们目前的物理学中，我们无法分析地预测双摆的运动。但是为什么？嗯，运动不是简单的随机的，而是混乱的。我们如何证明混沌和随机是不一样的？这是我们稍后将寻求回答的问题。

三、无法确定解的混沌方程

        你可能看过一些关于混沌的标题，我敢打赌缩略图看起来像这样：

        这幅画，以及许多其他蝴蝶翅膀式的图画，最常与混沌理论联系在一起。混沌理论中通常不被认为是混沌理论一部分的东西是双摆、三摆和任何其他具有不止一只臂的钟摆。让我们探索混沌运动的一些特性，并尝试在计算机程序的帮助下尽可能多地了解它，我们稍后将深入研究。

        我们可以研究的一个特性是，对起始条件进行微小的改变对最终结果的影响。最终，这种变化可以忽略不计，还是会产生一个完全来自原来的结果的答案？下面是 t=50 时双摆的计算机模拟，其中 θ1 和 θ2 的起始角分别为 80° 和 0°。（角度与上图相同）

        红线拖曳在钟摆的底部质量 （m2） 后面，并在绘制后保持 2 秒。

        现在，让我们运行完全相同的模拟，图像为 t=50，所有质量和臂长完全相同。但是，θ1 和 θ2 的角度分别为 79° 和 0°。θ1 的起始值发生了 1° 的变化，让我们看看它如何影响结果。

        不出所料（或者令人惊讶的是，我不知道你期望什么），仅仅1°的变化，系统就发生了巨大的变化，以至于系统的最终位置和它所追踪的线与第一张图像完全不同。这与任何外部变量无关，因为在我们的模拟中，没有摩擦、空气阻力或任何其他障碍。这个例子表明，当一个系统表现出混沌运动时，就像在双摆中一样，即使是最微小的条件变化也会导致结果的巨大变化。

        因此，我们已经确定了混沌运动的一个有趣的方面——启动条件的微小变化会导致结果的巨大变化。让我们试着再确定几个。如果我们再次进行同样的双摆实验，会发生什么？从实验上讲，这几乎是不可能的;这将涉及在测量重力场、质量、空气阻力和实验中所有其他变量等事物时具有无限的精度，然后使用您神奇的无限精度数据再次设置实验，完全相同。然而，我们有一个很好的模拟可以使用，这使得这个挑战消失了。

        让我们运行与之前类似的双摆模拟，但更改了一些变量，因此图像不相同：t=50，θ1=80°，θ2=15°，并且摆锤的物理性质相同（长度，质量）。结果如下：

        剧透警报，第二次完全一样......第三次、第四次、第五次和第六次。为了证明一个观点，让我们也尝试一下三重钟摆的练习。我们将使用相同的起始条件，尽管现在我们有第三个角度 θ3，因此我们将使用 θ1=80°、θ2=15°、θ3=0°。 不过，三重摆的计算量要大得多，因此图像将在 t=10 时拍摄。

        在三重摆锤模拟中，红线被拖到第三个质量的后面。

        尽管它比双摆模拟更复杂，但三摆在相同的起始条件下表现出相同的重复其模式的趋势这一事实很重要。在相同条件下获得相同结果的想法证明了混沌运动和混沌系统是确定性的，因此不是随机的。同样，由于这很难在实验中证明或在实践中使用，因此不会对我们对现实世界的理解产生太大影响。然而，我们现在确实知道混沌系统不是随机的，可以通过正确的知识进行预测。此外，由于我们证明了微小的变化可以产生很大的影响，我们知道需要非常高的精度才能有效地预测混沌系统的结果。

        我们现在已经确定了混沌运动的两个重要特性：第一个是微小的变化可以对系统的结果产生巨大的影响，第二个是混沌运动不是随机的，而是确定性的。让我们试着再找一个——正如我们在文章开头所讨论的，在双摆中你真正无法求解的变量是角位移。然而，我们也说过，求解周期是徒劳的，因为双钟摆不会重复——让我们来检验一下这个理论。

        首先，让我们在以下条件下运行双摆实验：θ1=30°，θ2=30°，质量和长度与之前相同。这一次，让我们使用 VPython 的绘图功能来记录底部质量高度随时间变化的图表：

        如果你问我，这看起来很周期性。那么，这是否意味着上述系统不是一个混沌系统呢？如果运动不是混沌的，这是否意味着当双摆的运动能量小于变得混沌所需的阈值能量时，我们可以得出一个解析解？我将把它留给读者作为练习。不过，要做到这一点，你可能需要我在本文中编写和使用的程序的帮助，所以让我们来看看它。这个程序是用 GlowScript 用 VPython 编写的，所以请确保你在 WebVPython 系统上运行它，而不仅仅是任何 python 编辑器。

四、计算机模拟

        首先，让我们定义将在程序中使用的物理常量。这些包括质量、手臂长度和重力。此外，我们将定义系统的起始条件，其中 theta1 和 theta2 是起始角，theta1dot 和 theta2dot 是臂的起始角速度。最后，我们还将开始时间设置为 0，并定义 dt，这是我们稍后在循环的每次迭代中向前迈进的时间量。

#lengths of the strings
L1 = 1.5 #top string
L2 = 1 #middle string

#masses
M1 = 2 #top moving ball
M2 = 1 #middle moving ball

#g
g = 9.8

#starting angles and velocities
theta1=30*pi/180 #release angle top ball
theta2= 30*pi/180 #release angle middle ball
theta1dot = 0
theta2dot = 0

t = 0
dt = 0.001

接下来，我们定义系统的所有部分，枢轴点，质量 1 （m1），臂 1 （stick1），第二组质量和棒 （m2， stick2） 也是如此。

pivot = sphere(pos=vector(0,L1,0), radius=0.05)

m1 = sphere(pos=pivot.pos-vector(0,L1,0),radius=0.05,color=color.white)
stick1 = cylinder(pos=pivot.pos,axis=m1.pos-pivot.pos,radius=0.015,color=color.yellow)

m2 = sphere(pos=m1.pos-vector(0,L2,0),radius=0.05,color=color.white)
stick2 = cylinder(pos=m1.pos, axis=m2.pos-m1.pos, radius=0.015,color=color.yellow)

        现在我们定位质量块和摇杆，使它们反映起始角度，就好像您运行代码之前它只是两根垂直杆一样。为此，我们只需使用一些三角函数将质量放在棍子将要到达的末端，然后将棍子设置为每个质量之间的距离。在这里，我们还使用 VPython 的“attach_trail”方法使底部质量生成轨迹。

m1.pos = pivot.pos+vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos+vector(L2*sin(theta2),-L2*cos(theta2),0)

stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos

attach_trail(m2, retain=200, color=color.red)

接下来，我们开始程序的主循环;在这个模拟中，与许多物理模拟一样，我们将使用一个 while 循环，其退出条件为 t=50。我们还将使用 rate（） 方法来设置程序的运行速度。

while t <= 50:  
    rate(1000)

        现在我们必须深入研究真正的数学。我们将使用一些微积分（不是真的，但我们确实处理导数）、三角学和代数来首先计算角加速度，然后用它来增加角速度，然后将该角速度应用于质量和杆的位置，以根据计算出的加速度移动它们。

        为了计算顶部质量的加速度，让我们用这个公式来表示角加速度，或θ1的二阶导数：

您还可以使用以下代码表示相同的方程：

theta1ddot = (-g*(2*M1 + M2)*sin(theta1) -M2*g*sin(theta1 -2*theta2)-2*sin(theta1 - theta2)*M2*(L2*theta2dot**2 + L1*cos(theta1 - theta2)*theta1dot**2))/(L1*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))

对于第二个质量的加速度，我们必须使用不同的公式，因为它代表系统的不同部分。

该等式在程序中写为：

theta2ddot = (2*sin(theta1 - theta2)*((M1 + M2)*L1*theta1dot**2 + g*(M1 + M2)*cos(theta1) + L2*M2*cos(theta1 - theta2)*theta2dot**2))/(L2*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))

        接下来，让我们递增角速度的变量 'theta1dot' 和 'theta2dot'。为此，我们只需使用下面的公式，该公式将角速度增加角加速度乘以dt。虽然图像只说明了 theta1dot 的方程，但 θ1 和 θ2 的方程是相同的。

        让我们用代码来表示这一点：

theta1dot = theta1dot + theta1ddot*dt
theta2dot = theta2dot + theta2ddot*dt

        对于我们计算双摆运动的最后一部分，必须将角度 θ1 和 θ2 增加角速度。这由下面的公式表示，其中我们将角度 θ 设置为等于前一个 θ 加上角速度乘以 dt。同样，θ1 和 θ2 的计算完全相同。

        在代码中，我们将等式写成：

theta1 = theta1 + theta1dot*dt
theta2 = theta2 + theta2dot*dt

        现在我们必须更新质量和杆在程序中的位置，以便我们可以看到程序的结果。首先，我们将更新群众的立场。我们将采用直观的方法来更新这些职位;由于系统中永远不会改变的一件事是杆的长度，因此中间质量 （M1） 将等于 L1 与枢轴点的距离。从逻辑上讲，我们可以为质量 m1 的位置构造以下方程：

        相同的公式可用于反映 m2 位置的公式，仅使用 θ2 和 L2，并且它基于 m1 的位置，而不是枢轴。请参阅下面的等式。

        这两个公式在我们的程序中编写得非常简单，因为 VPython 的一个特点是它使向量计算变得非常简单。所有对象的位置都已存储为向量，因此您需要做的就是使用 vector（） 方法定义一个新的向量，并在计算中使用它。

m1.pos = pivot.pos + vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos + vector(L2*sin(theta2),-L2*cos(theta2),0)

        现在我们必须迈出最后一步，通过更新杆的位置来实现我们的计算结果。从技术上讲，这并不是完全必要的，因为群众仍然会展示并遵循他们应该遵循的模式。不过，为了视觉效果，让我们编写这段代码。数学不是太有趣，所以代码如下。首先，我们设置第一个摇杆的轴，使其填充枢轴点和质量 m1 之间的距离。接下来，我们将第二根棍子的位置设置为等于质量 m1。最后，我们将第二根棍子的轴设置为等于第一根和第二根质量之间的距离。最后一行是按 dt 递增 t，使程序向前推进。

stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos

t += dt

以下是易于复制和粘贴的最终代码：

Web VPython 3.2

#lengths of the strings
L1 = 1.5 #top string
L2 = 1 #middle string

#masses
M1 = 2 #top moving ball
M2 = 1 #middle moving ball

#g
g = 9.8

t = 0
dt = 0.001

#starting angles and velocities
theta1=80.5*pi/180 #release angle top ball
theta2= 0*pi/180 #release angle middle ball
theta1dot = 0
theta2dot = 0

pivot = sphere(pos=vector(0,L1,0), radius=0.05)

m1 = sphere(pos=pivot.pos-vector(0,L1,0),radius=0.05,color=color.white)
stick1 = cylinder(pos=pivot.pos,axis=m1.pos-pivot.pos,radius=0.015,color=color.yellow)

m2 = sphere(pos=m1.pos-vector(0,L2,0),radius=0.05,color=color.white)
stick2 = cylinder(pos=m1.pos, axis=m2.pos-m1.pos, radius=0.015,color=color.yellow)


m1.pos = pivot.pos+vector(L1*sin(theta1),-L1*cos(theta1),0)
m2.pos = m1.pos+vector(L2*sin(theta2),-L2*cos(theta2),0)

stick1.axis = m1.pos - pivot.pos
stick2.pos = m1.pos
stick2.axis = m2.pos - m1.pos


attach_trail(m2, retain=200, color=color.red)

eChart = gdisplay(x=500, y=0, width=600, height=400, title="", xtitle="", ytitle="", foreground=color.black, background=color.white)
ePlot = gdots(color=color.red)

while t < 50:  
    rate(1000)
    
    #angular acceleration
    theta1ddot = (-g*(2*M1 + M2)*sin(theta1) -M2*g*sin(theta1 -2*theta2)-2*sin(theta1 - theta2)*M2*(L2*theta2dot**2 + L1*cos(theta1 - theta2)*theta1dot**2))/(L1*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))
    theta2ddot = (2*sin(theta1 - theta2)*((M1 + M2)*L1*theta1dot**2 + g*(M1 + M2)*cos(theta1) + L2*M2*cos(theta1 - theta2)*theta2dot**2))/(L2*(2*M1 + M2 - M2*cos(2*theta1 -2*theta2)))
    
    #angular velocity
    theta1dot = theta1dot + theta1ddot*dt
    theta2dot = theta2dot + theta2ddot*dt
    
    #angular position
    theta1 = theta1 + theta1dot*dt
    theta2 = theta2 + theta2dot*dt
    
    m1.pos = pivot.pos + vector(L1*sin(theta1),-L1*cos(theta1),0)
    m2.pos = m1.pos + vector(L2*sin(theta2),-L2*cos(theta2),0)
    stick1.axis = m1.pos - pivot.pos
    stick2.pos = m1.pos
    stick2.axis = m2.pos - m1.pos

    t = t+dt

        谢谢你走到这一步，我希望你喜欢我对双摆和三摆的混沌运动的分析。最后，我既不是专业的物理学家，也不是程序员，所以如果你对我如何改进我的工作有任何建议，请不要犹豫，告诉我。

        此外，本文中使用的所有图像都是由我使用 Python 和 Microsoft Word 创建的，两者都具有令人难以置信的功能。

        最后，我想以一个耐嚼的想法结束，这将在我的下一个项目中发挥作用。如果你拿一个已经经历过混沌运动的多臂摆，如果你要无缝地添加另一个质量非常高的臂（可能是系统中已经最大的 50 倍），起始 θ 为 0°，它的运动会有什么变化？