用Python的SciPy和Matplotlib，5分钟搞定单摆与双摆的混沌动画（附完整代码）

用Python打造单摆与双摆的混沌之美从基础模拟到炫酷动画你是否曾被物理课本上那些复杂的力学公式劝退却又对单摆和双摆的优雅运动轨迹充满好奇今天我们将用Python的SciPy和Matplotlib这两个强大的工具带你绕过繁琐的数学推导直接动手实现单摆和双摆的模拟与可视化。即使你完全不懂拉格朗日力学也能在短短几分钟内看到这些迷人的物理现象在你的屏幕上活灵活现。1. 准备工作与环境配置在开始之前我们需要确保你的Python环境已经安装了必要的库。如果你使用的是Anaconda这些库通常已经预装。如果没有可以通过以下命令快速安装pip install numpy scipy matplotlib这三个库将构成我们整个项目的基础NumPy提供高效的数值计算功能SciPy包含我们需要的微分方程求解器Matplotlib用于数据可视化和动画制作为了验证安装是否成功可以运行以下测试代码import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt print(所有库已正确安装)2. 单摆模拟从简单开始单摆是最基础的物理系统之一也是理解更复杂系统如双摆的绝佳起点。让我们先建立一个简单的单摆模型。2.1 单摆的物理模型单摆的运动可以用以下微分方程描述θ (g/L) * sin(θ) 0其中θ是摆角弧度g是重力加速度9.8 m/s²L是摆长米这个方程看起来简单但由于sin(θ)的非线性特性解析求解相当复杂。幸运的是我们可以用数值方法轻松解决。2.2 使用SciPy求解单摆运动下面是一个完整的单摆模拟代码from math import sin import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt # 物理参数 g 9.8 # 重力加速度 (m/s^2) L 1.0 # 摆长 (m) # 微分方程定义 def pendulum_eq(y, t, L): theta, omega y # 解包当前状态角度和角速度 dydt [omega, - (g/L) * sin(theta)] return dydt # 初始条件初始角度1弧度(约57度)初始角速度0 y0 [1.0, 0.0] # 时间点0到10秒步长0.01秒 t np.linspace(0, 10, 1000) # 求解微分方程 sol odeint(pendulum_eq, y0, t, args(L,)) # 提取结果 theta sol[:, 0] # 角度随时间变化 omega sol[:, 1] # 角速度随时间变化 # 绘制角度随时间变化图 plt.figure(figsize(10, 5)) plt.plot(t, theta, b, label角度(rad)) plt.plot(t, omega, g, label角速度(rad/s)) plt.legend(locbest) plt.xlabel(时间 (s)) plt.ylabel(角度/角速度) plt.title(单摆运动) plt.grid() plt.show()这段代码会生成一个展示单摆角度和角速度随时间变化的图表。你可以尝试修改初始角度y0中的第一个值观察不同初始条件下单摆的运动有何不同。2.3 单摆周期与非线性效应对于小角度摆动θ 0.5弧度单摆的周期近似为T ≈ 2π√(L/g)但当摆动角度增大时这个近似就不再准确。我们可以用SciPy计算不同初始角度下的精确周期from scipy.optimize import fsolve from scipy.special import ellipk def pendulum_period(L, theta0): 计算单摆的精确周期 return 4 * np.sqrt(L/g) * ellipk(np.sin(theta0/2)**2) # 计算不同初始角度下的周期 angles np.linspace(0, np.pi/2, 20) periods [pendulum_period(L, angle) for angle in angles] # 绘制周期与初始角度的关系 plt.figure(figsize(10, 5)) plt.plot(angles, periods, ro-) plt.xlabel(初始角度 (rad)) plt.ylabel(周期 (s)) plt.title(单摆周期与初始角度的关系) plt.grid() plt.show()3. 双摆模拟探索混沌世界双摆系统由两个单摆连接而成是一个经典的混沌系统。微小的初始条件变化可能导致完全不同的运动轨迹这就是所谓的蝴蝶效应。3.1 双摆的物理模型双摆的运动方程比单摆复杂得多。我们不需要自己推导这些方程这需要拉格朗日力学知识可以直接使用现成的公式class DoublePendulum: def __init__(self, m11.0, m21.0, l11.0, l21.0): self.m1, self.m2 m1, m2 # 两个小球的质量 self.l1, self.l2 l1, l2 # 两根杆的长度 self.g 9.8 # 重力加速度 def equations(self, state, t): 双摆的微分方程 theta1, z1, theta2, z2 state c, s np.cos(theta1-theta2), np.sin(theta1-theta2) # 方程组矩阵形式M * [theta1, theta2] F M np.array([ [(self.m1self.m2)*self.l1, self.m2*self.l2*c], [self.l1*c, self.l2] ]) F np.array([ -self.m2*self.l2*z2**2*s - (self.m1self.m2)*self.g*np.sin(theta1), self.l1*z1**2*s - self.g*np.sin(theta2) ]) # 解线性方程组得到角加速度 acc np.linalg.solve(M, F) return [z1, acc[0], z2, acc[1]]3.2 模拟双摆运动现在我们可以用这个类来模拟双摆的运动# 创建双摆实例 pendulum DoublePendulum(m11.0, m21.0, l11.0, l21.0) # 初始条件[theta1, omega1, theta2, omega2] y0 [np.pi/2, 0, np.pi/2, 0] # 时间点 t np.linspace(0, 20, 1000) # 求解微分方程 sol odeint(pendulum.equations, y0, t) # 提取结果 theta1, theta2 sol[:, 0], sol[:, 2] # 计算小球坐标 x1 pendulum.l1 * np.sin(theta1) y1 -pendulum.l1 * np.cos(theta1) x2 x1 pendulum.l2 * np.sin(theta2) y2 y1 - pendulum.l2 * np.cos(theta2) # 绘制轨迹 plt.figure(figsize(10, 10)) plt.plot(x1, y1, label上球轨迹) plt.plot(x2, y2, label下球轨迹) plt.legend() plt.axis(equal) plt.title(双摆运动轨迹) plt.grid() plt.show()尝试修改初始条件y0中的值观察双摆运动轨迹如何变化。你会发现即使初始角度只有微小差异长期运动轨迹也可能完全不同——这正是混沌系统的特征。4. 制作炫酷动画静态图像难以展现双摆运动的魅力让我们用Matplotlib的动画功能创建一个动态可视化。4.1 创建单摆动画首先是一个简单的单摆动画from matplotlib.animation import FuncAnimation # 重新计算单摆运动 t np.linspace(0, 10, 200) sol odeint(pendulum_eq, y0, t, args(L,)) theta sol[:, 0] x L * np.sin(theta) y -L * np.cos(theta) # 创建动画 fig, ax plt.subplots(figsize(8, 8)) ax.set_xlim(-1.2*L, 1.2*L) ax.set_ylim(-1.2*L, 1.2*L) ax.grid() line, ax.plot([], [], o-, lw2) trace, ax.plot([], [], b-, lw1, alpha0.5) time_text ax.text(0.02, 0.95, , transformax.transAxes) def init(): line.set_data([], []) trace.set_data([], []) time_text.set_text() return line, trace, time_text def animate(i): thisx [0, x[i]] thisy [0, y[i]] line.set_data(thisx, thisy) trace.set_data(x[:i], y[:i]) time_text.set_text(f时间 {t[i]:.1f}s) return line, trace, time_text ani FuncAnimation(fig, animate, frameslen(t), init_funcinit, blitTrue, interval20) plt.show()4.2 创建双摆动画双摆动画更加复杂但也更加迷人# 重新计算双摆运动 t np.linspace(0, 20, 500) sol odeint(pendulum.equations, y0, t) theta1, theta2 sol[:, 0], sol[:, 2] x1 pendulum.l1 * np.sin(theta1) y1 -pendulum.l1 * np.cos(theta1) x2 x1 pendulum.l2 * np.sin(theta2) y2 y1 - pendulum.l2 * np.cos(theta2) # 创建动画 fig, ax plt.subplots(figsize(10, 10)) ax.set_xlim(-2.2, 2.2) ax.set_ylim(-2.2, 2.2) ax.grid() line, ax.plot([], [], o-, lw2) trace1, ax.plot([], [], b-, lw1, alpha0.5) trace2, ax.plot([], [], r-, lw1, alpha0.5) time_text ax.text(0.02, 0.95, , transformax.transAxes) def init(): line.set_data([], []) trace1.set_data([], []) trace2.set_data([], []) time_text.set_text() return line, trace1, trace2, time_text def animate(i): thisx [0, x1[i], x2[i]] thisy [0, y1[i], y2[i]] line.set_data(thisx, thisy) trace1.set_data(x1[:i], y1[:i]) trace2.set_data(x2[:i], y2[:i]) time_text.set_text(f时间 {t[i]:.1f}s) return line, trace1, trace2, time_text ani FuncAnimation(fig, animate, frameslen(t), init_funcinit, blitTrue, interval20) plt.show()4.3 保存动画如果你想保存动画为GIF或视频文件可以取消注释以下代码# 保存为GIF需要安装pillow # ani.save(pendulum.gif, writerpillow, fps30) # 保存为MP4需要安装ffmpeg # ani.save(pendulum.mp4, writerffmpeg, fps30)5. 探索与扩展现在你已经掌握了单摆和双摆模拟的基本方法可以尝试以下扩展参数调整改变质量、摆长等参数观察系统行为变化能量分析计算并绘制系统的动能、势能和总能量相图绘制角度与角速度的关系图观察相空间轨迹三摆系统挑战更复杂的多摆系统交互式控制使用ipywidgets创建可交互的模拟界面例如这里是一个简单的能量计算函数def calculate_energies(theta1, theta2, omega1, omega2, m1, m2, l1, l2, g9.8): # 计算坐标和速度 x1 l1 * np.sin(theta1) y1 -l1 * np.cos(theta1) x2 x1 l2 * np.sin(theta2) y2 y1 - l2 * np.cos(theta2) # 速度分量 vx1 l1 * omega1 * np.cos(theta1) vy1 l1 * omega1 * np.sin(theta1) vx2 vx1 l2 * omega2 * np.cos(theta2) vy2 vy1 l2 * omega2 * np.sin(theta2) # 动能 T 0.5 * m1 * (vx1**2 vy1**2) 0.5 * m2 * (vx2**2 vy2**2) # 势能 V m1 * g * (y1 l1) m2 * g * (y2 l1 l2) return T, V, T V # 计算双摆的能量 T, V, E calculate_energies(theta1, theta2, sol[:,1], sol[:,3], pendulum.m1, pendulum.m2, pendulum.l1, pendulum.l2) # 绘制能量随时间变化 plt.figure(figsize(10, 5)) plt.plot(t, T, r, label动能) plt.plot(t, V, b, label势能) plt.plot(t, E, g, label总能量) plt.legend() plt.xlabel(时间 (s)) plt.ylabel(能量 (J)) plt.title(双摆能量变化) plt.grid() plt.show()通过这个项目你不仅学会了如何用Python模拟物理系统还掌握了科学计算和可视化的核心技能。这些技术可以应用于从工程分析到金融建模的广泛领域。最重要的是你看到了即使是最复杂的物理现象也可以通过计算工具变得触手可及。