恭喜您已经发现了这份神奇的计算材料学课件！这份课件是我在熬夜吃掉不计其数的披萨和咖啡后创作出来的，配套的教材是由单斌、陈征征、陈蓉合著的《计算材料学--从算法原理到代码实现》。学习资料合集您可以在这个网址找到：www.materialssimulation.com/book，您也可以跟着up主无人问津晦涩难懂的B站视频一起进行学习。希望它能帮您在计算材料学的道路上摔得不那么痛。

就像您尊重那些一边烘焙披萨一边写代码的大厨一样，当您使用这份课件时，请：

• 记得告诉大家这份课件是斯坦福大厨写的，并且他在华中科技大学微纳中心工作
• 别用它去赚大钱，这个课件是用来学习的，不是用来买披萨的
• 保持开放共享的精神

如果你有关于计算材料学的想法，或者你只是想和我讨论最好吃的披萨口味，欢迎通过邮件 bshan@mail.hust.edu.cn 联系我。

目录（对应《计算材料学》章节6.1.2）

点击购书

• Lyapunov不稳定性简介
• 三体问题Python求解
• 分子动力学中的Lyapunov不稳定性

##1. Lyapunov不稳定性简介

李亚普诺夫不稳定性是理解动态系统特性中一个非常关键的概念，尤其是在处理复杂系统如天体物理中的三体问题以及化学物理中的分子动力学模拟时。这种不稳定性表现为系统对初始条件的极度敏感性，即使是微乎其微的初始状态差异也能导致长时间演进后结果的巨大不同。这里我们将探讨李亚普诺夫不稳定性的影响，特别是在三体问题和分子动力学模拟中的具体体现。

三体问题中的李亚普诺夫不稳定性

三体问题是指在牛顿引力作用下，三个质量体互相作用的运动问题。尽管这个问题的表述相对简单，但其解决却极富挑战，因为系统显示出高度的混沌特征和敏感性。这意味着初始条件的轻微变化可以导致完全不同的运动轨迹，从而使得长期的精确预测变得不可能。这种不稳定性不仅理论上的有趣，而且对于理解天体物理系统的稳定性和演化具有重要意义。

分子动力学模拟中的李亚普诺夫不稳定性

在化学物理领域，分子动力学模拟是研究原子和分子交互作用的一种强有力的工具。通过模拟原子和分子的运动，科学家可以研究复杂系统的微观结构和动态过程。然而，由于李亚普诺夫不稳定性，模拟中的轨迹同样对初始条件极为敏感。这一特性意味着，尽管无法用分子动力学模拟来精确预测系统状态的长期演变，但可以通过对不同时间点的系统状态进行采样，来获得系统的统计特性，如温度、压力和能量分布等。

应对策略和科学意义

尽管李亚普诺夫不稳定性带来预测上的困难，但它也提供了探索系统动态行为的机会。在天体物理中，研究这种不稳定性有助于我们理解行星系统的形成和演化过程中可能出现的各种复杂交互作用。在分子动力学模拟中，通过运用合适的统计方法和增强采样技术，可以有效地探究和解释物质的性质，尽管存在不稳定性。此外，通过改进计算方法和算法，如引入并行计算和多尺度建模策略，可以进一步提高模拟的准确性和效率，从而在更宏观的尺度上洞察微观世界的复杂性。

总之，虽然李亚普诺夫不稳定性在理论和应用上都提出了挑战，它也极大地丰富了我们对动态系统的理解，使我们能够更全面地掌握自然界的复杂性和多样性。通过持续的研究和技术创新，我们可以更好地利用这一现象，深入探索从分子到宇宙的无数奥秘。

##2. 三体问题Python求解

三体问题，这一古老而复杂的天体物理难题，源自牛顿力学时代，至今仍吸引着众多科学家的研究兴趣。在中国科幻大师刘慈欣的科幻小说《三体》中，三体问题不仅仅是一个物理学问题，而是整个故事的核心元素。在这部作品中，三体世界由于其星系中存在三颗恒星的动态不稳定关系，使得其文明面临极端的生存挑战。这种不可预测的恒星运动导致三体世界的环境极端恶劣，时而酷热，时而严寒，给三体文明带来了一系列的生存危机。

在现实中，三体问题描述了三个质量体在相互引力作用下的运动规律。这个问题的复杂性在于，与两体问题相比，三体问题无法找到一个简洁的封闭形式解，因为系统的动态非常敏感，微小的初始条件变化都可能导致完全不同的轨迹演化。从数学角度看，三体问题是典型的混沌系统，具有高度敏感的初始条件依赖性和预测的不确定性。

本项目计划使用Python来实现三体问题的数值模拟，通过编写程序模拟三个天体在彼此引力作用下的动态演变。通过这样的实现，我们不仅可以更直观地理解三体问题的复杂性，还可以探索各种可能的动态系统行为，比如稳定的轨道、不规则运动甚至混沌状态。此外，借助《三体》的背景故事，我们将使这一古老科学问题更加生动有趣，增加读者的阅读兴趣，使科学探索与文学创作相结合，开启一次跨界的科学冒险之旅。

2.1 计算引力

此段代码定义了一个名为 gravitational_force 的函数，用于计算两个质量体之间的引力作用力向量。该函数接收三个参数：m1 和 m2 分别代表两个物体的质量，而 r_vec 是表示两个质量之间位置向量的 numpy 数组。

函数内部的处理步骤如下：

1. 首先，定义了引力常数 G 并简化设置为 1，以便于计算和理解。
2. 使用 numpy 的 norm 函数计算两质点之间的距离 r，这个距离是向量 r_vec 的模。
3. 接下来，根据牛顿万有引力定律计算力的大小 force_magnitude，计算公式为 G * m1 * m2 / r**2，即两质量乘积与距离平方的比值乘以引力常数。
4. 然后，计算力的方向 force_direction，它是单位向量，由原位置向量 r_vec 除以其模 r 得到。
5. 最后，函数返回负的力向量，即 -force_magnitude * force_direction。这是因为引力是吸引力，作用方向应当是指向两物体之间的直线上。

通过这个函数，可以方便地根据给定的质量和位置计算出两个物体之间的引力向量，这在模拟天体运动和其他需要考虑引力作用的物理系统中非常有用。

2.2 模拟三体运动的轨道

这段代码定义了一个函数 simulate_three_body，用于执行三体问题的数值模拟。这个函数通过求解牛顿运动方程来模拟三个质量体在彼此的引力作用下的动态演化。以下是代码的详细解释：

函数参数

• masses: 一个列表，包含了三个物体的质量。
• positions: 一个二维数组，初始时刻三个物体的位置坐标。
• velocities: 一个二维数组，初始时刻三个物体的速度向量。
• t_span: 一个元组，定义了时间模拟的起始和结束时间。
• num_points: 模拟中要计算的时间点数量，默认为500。

模拟过程

1. 内部函数 three_body_gravity：

   • 输入参数 t (当前时间) 和 y (当前所有物体的位置和速度的扁平化数组)。
   • 将 y 数组重新整形为位置和速度两部分，分别对应所有物体的当前位置和速度。
   • 初始化一个力的数组 forces，初始值为零，用于存放计算结果。
   • 通过双重循环计算每一对物体之间的引力作用，并更新相应的力数组。引力的计算使用之前定义的 gravitational_force 函数，根据牛顿万有引力定律和两物体之间的相对位置计算得出。
   • 根据质量调整计算得到的力，更新速度。
   • 将速度和加速度（由力除以质量得到）的扁平化数组合并后返回。

2. 初始条件设置：

   • 将初始位置和速度的数组扁平化合并为一个一维数组 y0。

3. 时间数组设置：

   • 使用 np.linspace 根据 t_span 和 num_points 生成模拟的时间点数组 t_eval。

4. 求解微分方程：

   • 使用 solve_ivp（初始值问题求解器）来求解定义好的微分方程。该求解器将 three_body_gravity 作为微分方程系统，从 t_span[0] 积分到 t_span[1]，使用 y0 作为初始条件。
   • 指定求解器的相对和绝对误差容忍度 rtol 和 atol 以控制计算精度。

返回值

• 返回 solve_ivp 的结果，其中包含了整个模拟期间每一个计算时间点的位置和速度。

这段代码为理解和研究天体物理中的三体问题提供了一个强有力的工具，通过数值方法能够模拟和观察三体系统随时间的动态变化，尽管系统可能表现出高度的混沌和不可预测性。

2.3 可视化轨迹

这段代码定义了一个名为 animate_simulation 的函数，用于将三体问题的数值模拟结果动画化展示。此函数不仅可以在屏幕上实时显示动画，还提供了选项保存动画为GIF文件。以下是对该函数的详细说明：

函数参数

• solution: 来自 solve_ivp 的模拟结果，包含了时间点和每个时间点对应的所有物体的位置。
• save_animation: 布尔类型参数，默认为 False。当设置为 True 时，会将动画保存为GIF文件。

动画创建流程

1. 创建绘图环境：

   • 使用 plt.subplots 创建一个图形和一个轴（fig 和 ax）。
   • 设置图形的大小为10x8英寸。
   • 为每个物体初始化一条线（使用列表推导），这些线将用于绘制物体的轨迹。

2. 设置图形属性：

   • 设置坐标轴的显示范围为x和y均为-6到6。
   • 设置图形的长宽比为1，确保x轴和y轴的单位比例相同。
   • 开启网格，以便于观察运动轨迹。

3. 初始化动画：

   • init 函数用于初始化动画，它将所有线条的数据设置为空，仅在动画开始时调用一次。

4. 定义动画更新方式：

   • animate 函数定义了每一帧动画如何更新，它根据 solution.y 中的数据更新每条线的数据，以此来展示三个物体的运动轨迹。
   • 参数 i 表示当前帧的索引，函数会根据这个索引从 solution.y 中提取每个物体对应的位置数据。

5. 生成动画对象：

   • 使用 FuncAnimation 创建一个动画对象，该对象会循环调用 animate 函数来更新图形，frames 参数设置为时间点的总数，以控制动画的总帧数。

6. 动画输出控制：

   • 如果 save_animation 为 True，则将动画保存为GIF文件。文件名为 "three_body_simulation.gif"，使用 pillow 作为写入器，帧率设置为每秒30帧。
   • 如果 save_animation 为 False，则直接在Notebook中渲染动画并显示。

2.4 执行模拟

这段示例代码展示了如何使用之前定义的 simulate_three_body 和 animate_simulation 函数来运行并可视化一个具体的三体问题模拟。代码详细地定义了初始条件，包括质量、位置、速度，并调整了初始速度以确保系统的总动量为零。接着执行模拟并进行动画展示。以下是详细解释和步骤：

定义初始条件

1. 质量：定义一个数组 masses，包含三个物体的质量，这里每个物体的质量都设为1.0。
2. 位置：定义一个二维数组 positions，包含每个物体的初始位置。本示例中，三个物体分别位于 (1.0, 0)、(-1.0, 0) 和 (0, 1.0)。
3. 速度：定义一个二维数组 velocities，包含每个物体的初始速度。为了展示动态效果，三个物体分别具有向正y轴、负y轴和负x轴的初速度。

确保总动量为零

• 计算总质心速度 center_of_mass_velocity，这是所有物体速度的质量加权平均。
• 通过从每个物体的速度中减去质心速度，调整初始速度以确保系统的总动量为零。这一步骤是为了避免整个系统在空间中漂移，使得模拟更加真实反映相对运动。

执行模拟

• 使用之前定义的 simulate_three_body 函数进行模拟。设置时间跨度 t_span 为从0到150，定义在此期间内计算1000个时间点的状态，这样可以获得足够平滑的动画效果。

可视化动画

• 调用 animate_simulation 函数生成模拟的动画。此函数根据模拟结果创建动态展示，使得观察者可以清楚地看到三体间相互作用的动态过程。
• 如果希望保存动画，可以将 save_animation 参数设置为 True。

显示动画

• 最后，生成的动画以HTML视频的形式输出，可以直接在支持HTML的环境中查看，如Jupyter Notebook。

改变初始条件,body 3的y方向坐标从1.0变成0.99,重新执行模拟。

##3. 分子动力学中的Lyapunov不稳定性

李亚普诺夫不稳定性在分子动力学（MD）模拟中扮演着重要角色，其揭示了微小的初始条件变化能如何引发系统行为的显著不同。这一现象虽然在初看似乎是一个挑战，实际上为我们理解和改进模拟提供了重要的视角和工具。

在分子动力学模拟中，李亚普诺夫不稳定性表现为原子和分子的轨迹对初始位置和速度极为敏感，即使是极小的差异也可能导致演化过程中完全不同的动态行为。这种敏感性意味着预测精确的单个轨迹变得非常困难，特别是在长时间的模拟中。然而，这并不妨碍分子动力学模拟作为一种强有力的工具来研究材料的物理性质、生物分子的动态行为以及化学反应过程，因为这些研究更侧重于统计平均的物理量而不是单个实体的精确位置或轨迹。

此外，李亚普诺夫不稳定性还启示我们，通过增强采样技术和改进模拟算法，可以更有效地探索系统的状态空间，尤其是在面对高维复杂系统时。例如，温度重复、距离增强采样等方法可以帮助系统跨越能量障碍，访问那些在常规MD模拟中难以达到的状态，从而获得更全面的系统行为描述。

总之，虽然李亚普诺夫不稳定性在分子动力学模拟中引入了预测上的不确定性，但通过正确的理解和应用，它反而成为了揭示分子系统复杂动态的一个强有力的工具。这种对初始条件敏感的性质促使科学家开发新的模拟策略，更深入地理解和预测材料和生物分子的行为。