1. 从单摆到双摆一个迷人的混沌世界如果你玩过单摆或者至少见过钟摆那你对规则、可预测的周期性运动一定不陌生。但当你把两个单摆用铰链连接起来事情就变得完全不一样了。这个看似简单的物理系统——双摆是展示混沌动力学最经典、最直观的例子之一。没有复杂的数学公式仅仅通过观察它的运动你就能亲眼目睹“确定性混沌”一个完全由确定性方程描述的系统却对初始条件有着极端的敏感性微小的扰动会导致运动轨迹在短时间内变得截然不同、无法长期预测。用Python来模拟和可视化双摆是一个绝佳的练手项目。它完美融合了物理建模、数值计算和图形可视化。你不仅能重温牛顿力学和拉格朗日方程还能亲手用代码“创造”出一个混沌系统并把它画出来看着它如何从有序走向无序。这对于学习Python的科学计算栈如NumPy, SciPy和绘图库如Matplotlib是极好的实践。更重要的是整个过程充满了探索的乐趣你会不断调整参数观察系统如何响应就像在做一个数字物理实验。2. 理论基石推导双摆的运动方程在动手写代码之前我们必须搞清楚双摆到底遵循什么样的物理规律。直接使用牛顿力学分析铰链处的约束力会非常繁琐。这里拉格朗日力学是我们的得力工具它通过能量来推导运动方程完美规避了约束力的复杂计算。2.1 建立模型与坐标系我们考虑一个理想化的平面双摆模型两个摆杆都是刚性的质量分别为m1和m2长度分别为L1和L2。第一个摆的顶端固定在原点(0, 0)其与竖直向下方向的夹角为θ1。第二个摆的顶端也就是第一个摆的末端与第一个摆铰接其与竖直向下方向的夹角为θ2。忽略所有形式的摩擦空气阻力、铰链摩擦。根据这个模型我们可以写出两个摆锤的位置坐标摆锤1的位置x1 L1 * sin(θ1),y1 -L1 * cos(θ1)这里y轴向下为正方便绘图摆锤2的位置x2 x1 L2 * sin(θ2),y2 y1 - L2 * cos(θ2)2.2 计算系统的动能与势能系统的动能T是两个摆锤动能之和。摆锤的速度是其位置对时间的导数。v1^2 (dx1/dt)^2 (dy1/dt)^2 (L1 * θ1)^2因为θ1 dθ1/dtv2^2的计算稍复杂因为摆锤2的位置依赖于θ1和θ2v2^2 (L1 * θ1 * cos(θ1) L2 * θ2 * cos(θ2))^2 (L1 * θ1 * sin(θ1) L2 * θ2 * sin(θ2))^2化简后可得v2^2 L1^2 * θ1^2 L2^2 * θ2^2 2 * L1 * L2 * θ1 * θ2 * cos(θ1 - θ2)因此总动能为T 0.5 * m1 * v1^2 0.5 * m2 * v2^2系统的势能V则以悬挂点为参考点势能为0V m1 * g * y1 m2 * g * y2 -m1 * g * L1 * cos(θ1) - m2 * g * (L1 * cos(θ1) L2 * cos(θ2))这里g是重力加速度。2.3 应用拉格朗日方程拉格朗日量L定义为动能与势能之差L T - V。 对于每个广义坐标θ1和θ2拉格朗日方程为d/dt (∂L/∂θi) - ∂L/∂θi 0, 其中i 1, 2。将L T - V代入并执行这些微分运算这个过程比较冗长通常借助符号计算软件如SymPy来完成我们可以得到两个耦合的二阶微分方程它们描述了θ1和θ2随时间t的演化(m1 m2) * L1 * θ1 m2 * L2 * θ2 * cos(θ1 - θ2) m2 * L2 * θ2^2 * sin(θ1 - θ2) (m1 m2) * g * sin(θ1) 0 m2 * L2 * θ2 m2 * L1 * θ1 * cos(θ1 - θ2) - m2 * L1 * θ1^2 * sin(θ1 - θ2) m2 * g * sin(θ2) 0其中θi表示角加速度。这两个方程就是我们的“游戏规则”代码的任务就是数值求解它们。提示在实际编程中我们很少直接处理这两个二阶方程。标准的做法是将它们转化为一个一阶微分方程组4个方程因为大多数数值积分器如scipy.integrate.solve_ivp都是为求解一阶方程组设计的。定义状态向量u [θ1, θ2, ω1, ω2]其中ω1 θ1,ω2 θ2。那么我们的任务就是写出du/dt [ω1, ω2, α1, α2]的表达式其中α1和α2需要从上面两个二阶方程中联立解出。3. 代码实现从方程到动画理论准备就绪现在进入实战环节。我们将使用NumPy进行数值计算SciPy进行微分方程积分Matplotlib进行可视化并用它的动画模块让双摆动起来。3.1 环境准备与核心函数定义首先确保你的Python环境安装了必要的库。打开终端或命令提示符运行pip install numpy scipy matplotlib接下来我们开始编写核心代码。创建一个新的Python文件比如double_pendulum.py。import numpy as np from scipy.integrate import solve_ivp import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation # 定义双摆参数 m1, m2 1.0, 1.0 # 质量 (kg) L1, L2 1.0, 1.0 # 长度 (m) g 9.81 # 重力加速度 (m/s^2) # 初始条件: [theta1, theta2, omega1, omega2] # 尝试微小的差异观察混沌 initial_state_1 [np.pi / 2, np.pi / 2, 0.0, 0.0] # 两个摆都水平释放 initial_state_2 [np.pi / 2 0.001, np.pi / 2, 0.0, 0.0] # 仅第一个角度有千分之一弧度的差异 # 模拟总时间和时间步 t_max 30.0 # 模拟总时长 (秒) fps 30 # 动画帧率 t_eval np.linspace(0, t_max, int(t_max * fps)) # 用于动画渲染的时间点 def derivs(t, state): 计算状态向量 state [theta1, theta2, omega1, omega2] 的导数。 返回 dstate/dt [omega1, omega2, alpha1, alpha2]。 theta1, theta2, omega1, omega2 state # 一些中间计算避免重复 delta theta1 - theta2 sin1, cos1 np.sin(theta1), np.cos(theta1) sin2, cos2 np.sin(theta2), np.cos(theta2) sind, cosd np.sin(delta), np.cos(delta) # 方程组系数的分母部分 denom1 (m1 m2) * L1 - m2 * L1 * cosd**2 denom2 L2 - (L1 * cosd**2 * L2) / ((m1 m2) * L1 / m2) # 简化形式实际需联立求解 # 更稳健的解法构建线性方程组 A * [alpha1, alpha2] b # 根据推导出的二阶方程形式 # a11 * alpha1 a12 * alpha2 b1 # a21 * alpha1 a22 * alpha2 b2 a11 (m1 m2) * L1 a12 m2 * L2 * cosd b1 -m2 * L2 * omega2**2 * sind - (m1 m2) * g * sin1 a21 L1 * cosd a22 L2 b2 L1 * omega1**2 * sind - g * sin2 # 解线性方程组 A np.array([[a11, a12], [a21, a22]]) B np.array([b1, b2]) # 使用线性代数求解比直接写显式公式更清晰且不易出错 alpha np.linalg.solve(A, B) return [omega1, omega2, alpha[0], alpha[1]]这段代码的核心是derivs函数。我在这里选择了一种更清晰、数值稳定性更好的方法将运动方程整理成关于角加速度α1,α2的线性方程组A * α B然后用np.linalg.solve求解。这种方法避免了手动推导并输入那两个复杂的显式表达式减少了出错的可能尤其是在参数如质量变化时更具通用性。3.2 数值积分与轨迹计算有了导数函数我们就可以用solve_ivp来模拟双摆的运动了。# 对两组略微不同的初始条件进行积分以对比混沌效应 print(正在积分第一个双摆轨迹...) sol1 solve_ivp(derivs, [0, t_max], initial_state_1, t_evalt_eval, methodRK45, rtol1e-9, atol1e-9) print(正在积分第二个双摆轨迹...) sol2 solve_ivp(derivs, [0, t_max], initial_state_2, t_evalt_eval, methodRK45, rtol1e-9, atol1e-9) # 解包结果 theta1_1, theta2_1 sol1.y[0, :], sol1.y[1, :] theta1_2, theta2_2 sol2.y[0, :], sol2.y[1, :] # 计算两个摆锤的笛卡尔坐标轨迹 x1_1 L1 * np.sin(theta1_1) y1_1 -L1 * np.cos(theta1_1) x2_1 x1_1 L2 * np.sin(theta2_1) y2_1 y1_1 - L2 * np.cos(theta2_1) x1_2 L1 * np.sin(theta1_2) y1_2 -L1 * np.cos(theta1_2) x2_2 x1_2 L2 * np.sin(theta2_2) y2_2 y1_2 - L2 * np.cos(theta2_2)这里我使用了高精度的RK45方法即Runge-Kutta 4/5阶并设置了很低的容差rtol,atol。对于混沌系统数值积分误差本身也会被系统放大因此使用高精度积分器并记录尽可能多的时间点数据是必要的。t_eval参数确保了我们在均匀的时间点上获取解这对于生成流畅的动画至关重要。3.3 创建静态轨迹图与相位图在制作动画前先绘制静态图来直观感受双摆的复杂运动和混沌特性。# 创建画布和子图 fig plt.figure(figsize(14, 10)) # 1. 双摆轨迹对比图 ax1 plt.subplot(2, 2, 1) ax1.plot(x2_1, y2_1, b-, alpha0.6, lw0.8, label轨迹1 (θ1π/2)) ax1.plot(x2_2, y2_2, r-, alpha0.6, lw0.8, label轨迹2 (θ1π/20.001)) ax1.set_xlabel(x (m)) ax1.set_ylabel(y (m)) ax1.set_title(第二个摆锤的轨迹对比 (混沌敏感性)) ax1.legend() ax1.grid(True, alpha0.3) ax1.set_aspect(equal, adjustablebox) ax1.set_xlim(-2.5, 2.5) ax1.set_ylim(-2.5, 1) # 2. 角度随时间变化图 ax2 plt.subplot(2, 2, 2) ax2.plot(sol1.t, theta1_1, b-, labelθ1 (轨迹1)) ax2.plot(sol1.t, theta2_1, b--, labelθ2 (轨迹1)) ax2.plot(sol2.t, theta1_2, r-, labelθ1 (轨迹2)) ax2.plot(sol2.t, theta2_2, r--, labelθ2 (轨迹2)) ax2.set_xlabel(时间 (s)) ax2.set_ylabel(角度 (rad)) ax2.set_title(摆角随时间变化) ax2.legend() ax2.grid(True, alpha0.3) # 3. 相位空间图 (θ1, ω1) ax3 plt.subplot(2, 2, 3) ax3.plot(theta1_1, sol1.y[2, :], b-, alpha0.7, lw0.5) ax3.set_xlabel(θ1 (rad)) ax3.set_ylabel(ω1 (rad/s)) ax3.set_title(第一个摆的相位空间 (θ1, ω1)) ax3.grid(True, alpha0.3) # 4. 两个轨迹在相空间中的分离 ax4 plt.subplot(2, 2, 4) # 计算两个系统状态向量的欧氏距离这里简化用角度差 theta_diff np.sqrt((theta1_1 - theta1_2)**2 (theta2_1 - theta2_2)**2) ax4.semilogy(sol1.t, theta_diff, k-) ax4.set_xlabel(时间 (s)) ax4.set_ylabel(角度差 (log scale)) ax4.set_title(两个轨迹的分离 (对数坐标)) ax4.grid(True, alpha0.3) ax4.set_ylim(bottom1e-6) # 设置y轴下限避免log(0)的问题 plt.tight_layout() plt.show()轨迹图展示了第二个摆锤末端划过的路径。两条轨迹初始仅相差0.001弧度但在短时间内就分道扬镳这正是“蝴蝶效应”的直观体现。角度-时间图显示了两个摆角如何随时间演变。它们看起来毫无规律且两组初始条件对应的曲线很快就不再重合。相位图绘制了第一个摆的角位置θ1和角速度ω1的关系。在混沌系统中相位轨迹不会形成闭合的环或稳定的点而是会在一个区域内无限地、永不重复地缠绕。分离图以对数坐标展示两个系统状态差异随时间指数增长的过程。在初期这条线大致是直线其斜率对应于系统的李雅普诺夫指数正指数是混沌的标志。3.4 制作动态动画静态图展示了结果但动画才能生动呈现双摆运动的魅力。我们将使用Matplotlib的FuncAnimation。# 为动画创建新的图形 fig_anim, ax_anim plt.subplots(figsize(8, 8)) ax_anim.set_xlim(-(L1 L2 0.5), L1 L2 0.5) ax_anim.set_ylim(-(L1 L2 0.5), L1 L2 0.5) ax_anim.set_aspect(equal) ax_anim.grid(True, alpha0.3) ax_anim.set_title(双摆混沌运动模拟) # 初始化动画元素 line, ax_anim.plot([], [], o-, lw3, colorblue, markersize10) # 摆杆和摆锤 trace, ax_anim.plot([], [], -, lw1, colorblue, alpha0.5) # 第二个摆锤的轨迹 time_text ax_anim.text(0.02, 0.95, , transformax_anim.transAxes, fontsize12) # 轨迹数据容器 trace_x, trace_y [], [] def init(): 初始化动画函数 line.set_data([], []) trace.set_data([], []) time_text.set_text() return line, trace, time_text def animate(i): 每一帧的更新函数 # 当前帧对应的索引 idx i * 2 # 如果帧数太多可以跳帧以加速动画生成 if idx len(sol1.t): idx len(sol1.t) - 1 # 计算当前帧的摆杆坐标 x_coords [0, x1_1[idx], x2_1[idx]] y_coords [0, y1_1[idx], y2_1[idx]] # 更新摆杆线条 line.set_data(x_coords, y_coords) # 更新轨迹只保留最近N个点避免线条过长影响性能 trace_x.append(x2_1[idx]) trace_y.append(y2_1[idx]) keep_points 500 # 保留的轨迹点数 if len(trace_x) keep_points: trace_x.pop(0) trace_y.pop(0) trace.set_data(trace_x, trace_y) # 更新时间显示 time_text.set_text(f时间: {sol1.t[idx]:.2f} s) return line, trace, time_text # 创建动画对象 # frames 控制总帧数interval 控制帧间延迟(ms)blit 优化渲染 ani FuncAnimation(fig_anim, animate, frameslen(t_eval)//2, init_funcinit, interval1000/fps, blitTrue, repeatTrue) # 保存动画为GIF需要安装pillow # print(正在保存动画这可能需要一些时间...) # ani.save(double_pendulum_chaos.gif, writerpillow, fpsfps, dpi100) plt.show()运行这段代码你将看到一个动态的双摆。观察它如何从简单的周期性摆动如果初始角度很小迅速演变成疯狂的、不可预测的旋转和甩动。轨迹线会逐渐填满一个环形区域这个区域被称为“混沌吸引子”的投影。4. 探索、优化与深度分析一个能动的双摆模拟已经完成了但作为项目我们可以走得更远。4.1 参数空间探索寻找不同的运动模式双摆的行为严重依赖于初始条件。我们可以写一个简单的循环来批量模拟并分类不同的运动形态。def simulate_and_classify(initial_state, t_span[0, 50], max_speed_threshold20.0): 模拟并粗略分类运动模式 sol solve_ivp(derivs, t_span, initial_state, methodRK45, max_step0.01) theta1, theta2, omega1, omega2 sol.y # 简单分类逻辑 avg_energy np.mean(0.5*m1*L1**2*omega1**2 0.5*m2*(L1**2*omega1**2 L2**2*omega2**2 2*L1*L2*omega1*omega2*np.cos(theta1-theta2))) - (m1m2)*g*L1*np.cos(theta1) - m2*g*L2*np.cos(theta2) max_speed np.max(np.abs(omega1)) if max_speed max_speed_threshold: return 混沌高速旋转 elif np.std(theta1) 0.5 and np.std(theta2) 0.5: # 角度波动小 return 近似规则摆动 else: return 混沌复杂摆动 # 测试不同初始角度 initial_angles [(0.1, 0.1), (np.pi/2, 0), (np.pi-0.1, 0), (np.pi/2, np.pi/2), (3.0, 1.0)] for i, (th1, th2) in enumerate(initial_angles): state [th1, th2, 0.0, 0.0] mode simulate_and_classify(state) print(f初始条件 {i1}: (θ1{th1:.2f}, θ2{th2:.2f}) - 运动模式: {mode})这个简单的分类器可以帮助你快速了解哪些初始条件会导致相对温和的摆动哪些会立刻引发狂暴的混沌。你会发现当初始角度接近倒立位置π时系统极其不稳定极易进入混沌。4.2 性能优化与长时间模拟上面的代码对于短时间模拟是没问题的但如果你想模拟几分钟甚至更长时间或者需要极高的时间分辨率性能就会成为瓶颈。这里有几个优化方向使用更高效的积分器solve_ivp的DOP853或Radau方法对于某些刚性或高精度问题可能更快更稳定。向量化与预计算在derivs函数中我们对每个时间点单独计算三角函数。如果使用自定义的固定时间步长积分循环如4阶Runge-Kutta可以预先计算所有时间步的三角函数向量但会牺牲一些灵活性。使用 Numba 加速这是提升性能的大杀器。用numba.jit(nopythonTrue)装饰你的derivs函数可以将其编译成机器码速度提升数十倍甚至上百倍。# 示例使用Numba加速需要先 pip install numba import numba numba.jit(nopythonTrue) def derivs_numba(t, state, m1, m2, L1, L2, g): theta1, theta2, omega1, omega2 state delta theta1 - theta2 sin1, cos1 np.sin(theta1), np.cos(theta1) sin2, cos2 np.sin(theta2), np.cos(theta2) sind, cosd np.sin(delta), np.cos(delta) a11 (m1 m2) * L1 a12 m2 * L2 * cosd b1 -m2 * L2 * omega2**2 * sind - (m1 m2) * g * sin1 a21 L1 * cosd a22 L2 b2 L1 * omega1**2 * sind - g * sin2 # 手动解2x2线性方程组 (克莱姆法则)在Numba中比调用np.linalg.solve更高效 det a11 * a22 - a12 * a21 alpha1 (b1 * a22 - a12 * b2) / det alpha2 (a11 * b2 - b1 * a21) / det return np.array([omega1, omega2, alpha1, alpha2]) # 需要将参数打包传递给 solve_ivp使用 args 参数 sol_fast solve_ivp(lambda t, y: derivs_numba(t, y, m1, m2, L1, L2, g), [0, t_max], initial_state_1, t_evalt_eval, methodRK45, rtol1e-9)在我的测试中使用Numba加速后积分速度可以提升50倍以上对于需要大量模拟或长时间模拟的场景至关重要。4.3 能量守恒检验验证模拟的准确性在一个无摩擦的理想系统中总机械能动能势能应该守恒。由于数值积分会引入微小误差检查能量是否漂移是验证我们模拟是否正确、积分精度是否足够的一个好方法。def total_energy(state, m1, m2, L1, L2, g): 计算给定状态下的总机械能 theta1, theta2, omega1, omega2 state # 动能 v1_sq (L1 * omega1)**2 v2_sq L1**2 * omega1**2 L2**2 * omega2**2 2 * L1 * L2 * omega1 * omega2 * np.cos(theta1 - theta2) T 0.5 * m1 * v1_sq 0.5 * m2 * v2_sq # 势能 (以悬挂点为0点) y1 -L1 * np.cos(theta1) y2 y1 - L2 * np.cos(theta2) V m1 * g * y1 m2 * g * y2 return T V # 计算模拟过程中每一时刻的能量 energy_1 np.array([total_energy(sol1.y[:, i], m1, m2, L1, L2, g) for i in range(len(sol1.t))]) initial_energy energy_1[0] energy_error np.abs((energy_1 - initial_energy) / initial_energy) # 绘制能量相对误差 fig_energy, ax_energy plt.subplots(figsize(10, 5)) ax_energy.plot(sol1.t, energy_error * 100, b-) # 转换为百分比 ax_energy.set_yscale(log) # 对数坐标看得更清楚 ax_energy.set_xlabel(时间 (s)) ax_energy.set_ylabel(总能量相对误差 (%)) ax_energy.set_title(数值模拟的能量守恒检验 (对数坐标)) ax_energy.grid(True, alpha0.3) ax_energy.set_ylim(bottom1e-12) # 看看误差的下限 plt.show() print(f最大相对能量误差: {np.max(energy_error)*100:.2e}%) print(f平均相对能量误差: {np.mean(energy_error)*100:.2e}%)一个优秀的模拟其能量误差应该非常小例如小于1e-6%并且不会随时间系统性增长。如果你发现能量漂移很大可能需要检查运动方程是否正确或者降低solve_ivp的容差rtol,atol。5. 从模拟到艺术可视化技巧与扩展思路基础功能实现后我们可以让可视化更具表现力和启发性。5.1 多摆对比与轨迹着色同时模拟多个初始条件略有不同的双摆并用颜色区分它们的轨迹可以非常直观地展示“轨迹指数分离”这一混沌核心特征。# 模拟5个初始角度有细微差异的双摆 num_pendulums 5 initial_theta1 np.linspace(np.pi/2 - 0.02, np.pi/2 0.02, num_pendulums) colors plt.cm.viridis(np.linspace(0, 1, num_pendulums)) fig_multi, ax_multi plt.subplots(figsize(10, 10)) ax_multi.set_xlim(-2.2, 2.2) ax_multi.set_ylim(-2.2, 1.2) ax_multi.set_aspect(equal) ax_multi.set_title(多个双摆的轨迹分离 (混沌敏感性)) for i in range(num_pendulums): state [initial_theta1[i], np.pi/2, 0.0, 0.0] sol solve_ivp(derivs, [0, 15], state, methodRK45, max_step0.01) theta1, theta2 sol.y[0, :], sol.y[1, :] x2 L1 * np.sin(theta1) L2 * np.sin(theta2) y2 -L1 * np.cos(theta1) - L2 * np.cos(theta2) # 用颜色和透明度表示时间形成“彗星尾”效果 points np.array([x2, y2]).T.reshape(-1, 1, 2) segments np.concatenate([points[:-1], points[1:]], axis1) from matplotlib.collections import LineCollection lc LineCollection(segments, cmapviridis, normplt.Normalize(0, sol.t[-1]), alpha0.7) lc.set_array(sol.t) # 根据时间着色 lc.set_linewidth(1) ax_multi.add_collection(lc) ax_multi.plot(0, 0, ko, markersize10) # 绘制固定点 ax_multi.grid(True, alpha0.3) plt.colorbar(lc, axax_multi, label时间 (s)) plt.show()这张图会显示最初几乎重叠的几条轨迹在短短几秒后就像爆炸一样飞散到整个可达区域这就是对初始条件敏感依赖的完美视觉证明。5.2 三维可视化与庞加莱截面对于高阶爱好者可以尝试更高级的可视化。三维相空间将两个角度和一个角速度如θ1,θ2,ω1作为三个坐标轴绘制轨迹可以更完整地观察系统在相空间中的运动。庞加莱截面这是一种研究高维相空间结构的强大工具。对于双摆四维相空间我们可以固定一个变量如θ10然后每当轨迹穿过这个“截面”时就记录下其他三个变量θ2,ω1,ω2的值并画出来。在混沌系统中庞加莱截面上的点会形成复杂的分形结构而不是清晰的曲线或点集。实现这个需要记录轨迹穿越事件对编程是个不错的挑战。5.3 将项目封装成类与应用为了让代码更整洁、可复用可以将其封装成一个DoublePendulum类。class DoublePendulum: def __init__(self, m11.0, m21.0, L11.0, L21.0, g9.81): self.params (m1, m2, L1, L2, g) def equations(self, t, state): # 使用前面定义的 derivs 或 derivs_numba 逻辑 m1, m2, L1, L2, g self.params theta1, theta2, omega1, omega2 state delta theta1 - theta2 cos_delta np.cos(delta) sin_delta np.sin(delta) # 构建线性方程组系数矩阵 A np.array([ [(m1 m2) * L1, m2 * L2 * cos_delta], [L1 * cos_delta, L2] ]) B np.array([ -m2 * L2 * omega2**2 * sin_delta - (m1 m2) * g * np.sin(theta1), L1 * omega1**2 * sin_delta - g * np.sin(theta2) ]) alpha np.linalg.solve(A, B) return [omega1, omega2, alpha[0], alpha[1]] def simulate(self, initial_state, t_span, t_evalNone): sol solve_ivp(self.equations, t_span, initial_state, t_evalt_eval, methodRK45, rtol1e-9) self.solution sol self.t sol.t self.theta1, self.theta2, self.omega1, self.omega2 sol.y self.x1 self.params[2] * np.sin(self.theta1) self.y1 -self.params[2] * np.cos(self.theta1) self.x2 self.x1 self.params[3] * np.sin(self.theta2) self.y2 self.y1 - self.params[3] * np.cos(self.theta2) return sol # 使用类 dp DoublePendulum(m11.2, m20.8, L11.0, L20.8) dp.simulate([np.pi/3, np.pi/4, 0, 0], [0, 20]) # 现在可以方便地访问 dp.x2, dp.y2 等属性进行绘图这个项目还可以扩展到三摆、多摆或者加入驱动和阻尼研究同步、混沌控制等更前沿的课题。你也可以用PyGame或Pyglet制作交互式模拟允许用户用鼠标拖动摆锤改变初始条件。