import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint ================== AFM系统动力学方程 ================== def afm_system(state, t, params): x1, x2 = state m, k, D, sigma, Z0, mu, f_amp, omega = params F_spring = -k*x1/m F_vdW = D/(m*(Z0 + x1)**2) F_repulsive = (sigma**6 * D)/(30*m*(Z0 + x1)**8) F_excitation = (f_amp/m) * np.cos(omega*t) dx1 = x2 dx2 = F_spring - F_vdW + F_repulsive - (mu/m)*x2 + F_excitation return [dx1, dx2] ================== 参数设置 ================== params = [ 1e-12, # m [kg] 0.1, # k [N/m] 1e-19, # D [J·m] 3e-9, # sigma [m] 4.5e-9, # Z0 [m] 5e-8, # mu [kg/s] 2e-9, # f_amp [N] 2np.pi1e4 # omega [rad/s] ] ================== 分岔图计算 ================== f_amp_range = np.linspace(1e-10, 3e-9, 200) # 力幅扫描范围 max_x1 = [] for f_amp in f_amp_range: params[6] = f_amp # 更新激励幅值 # 瞬态过程（前200周期） t_transient = np.linspace(0, 200*(2np.pi/params[7]), 2000) sol = odeint(afm_system, [1e-9, 0], t_transient, args=(params,)) # 稳态采样（后50周期） t_steady = np.linspace(0, 50(2*np.pi/params[7]), 500) sol = odeint(afm_system, sol[-1], t_steady, args=(params,)) max_x1.append(np.max(sol[:,0])) ================== 庞加莱截面计算 ================== params[6] = 2.5e-9 # 固定激励幅值 poincare_section = [] total_time = 1000*(2*np.pi/params[7]) # 1000个周期 t = np.linspace(0, total_time, 200000) sol = odeint(afm_system, [1e-9, 0], t, args=(params,)) 每个周期取1个点 poincare_idx = np.arange(80000, 100000, int(2*np.pi/params[7]/(t[1]-t[0]))) poincare_section = sol[poincare_idx] ================== 科学可视化设置 ================== plt.figure(figsize=(15, 6), dpi=100) 分岔图 plt.subplot(1, 2, 1) plt.plot(f_amp_range/1e-9, np.array(max_x1)/1e-9, ‘b.’, markersize=3, alpha=0.6) plt.xlabel(r’ f a m p ( 10 − 9 N ) f amp ​ (10 −9 N)‘, fontsize=12) plt.ylabel(r’ m a x ( x 1 ) ( 10 − 9 m ) max(x 1 ​ ) (10 −9 m)’, fontsize=12) plt.title(‘Bifurcation Diagram’, fontsize=14) plt.grid(True, alpha=0.3) plt.gca().xaxis.set_major_formatter(plt.FormatStrFormatter(‘%.1f’)) plt.gca().yaxis.set_major_formatter(plt.FormatStrFormatter(‘%.2f’)) 庞加莱截面 plt.subplot(1, 2, 2) plt.scatter(poincare_section[:,0]/1e-9, poincare_section[:,1]/1e-3, c=np.arange(len(poincare_section)), cmap=‘viridis’, s=10, alpha=0.7, edgecolors=‘none’) plt.xlabel(r’ x 1 ( 10 − 9 m ) x 1 ​ (10 −9 m)‘, fontsize=12) plt.ylabel(r’ x 2 ( 10 − 3 m / s ) x 2 ​ (10 −3 m/s)’, fontsize=12) plt.title(‘Poincaré Section’, fontsize=14) plt.grid(True, alpha=0.3) cbar = plt.colorbar() cbar.set_label(‘Iteration Index’, rotation=270, labelpad=15) plt.tight_layout() plt.savefig(‘afm_analysis.png’, bbox_inches=‘tight’, dpi=300) plt.show() 参数和动力学方程参考这段代码 之前你给我了相图重构和彭加莱登的代码 请结合动力学方程 把你之前给我的绘制这些图像的完整的python代码发给我

<think>好的，我现在需要帮用户整理关于AFM系统动力学分析的Python代码，并结合之前提供的相图重构和庞加莱截面的代码。首先，用户提供的现有代码存在一些语法错误，比如`2np.pi1e4`应该是`2 * np.pi * 1e4`，还有字符串中的引号格式不正确，比如使用了‘’而不是标准的单引号'。我需要先修正这些错误，确保代码能够正确运行。 接下来，用户希望将分岔图、庞加莱截面以及相图重构的代码整合在一起。之前提供的相图重构涉及到时间延迟嵌入，可能需要从解算后的时间序列数据中提取。因此，我需要检查现有的动力学方程求解部分是否正确，参数设置是否合理，并确保所有可视化部分的数据来源一致。 然后，分岔图部分是通过扫描激励幅值f_amp_range，计算每个f_amp对应的最大位移max_x1。这里需要注意循环中的ODE求解是否正确，特别是瞬态过程和稳态采样的时间设置是否合理，避免初始瞬态影响结果。 庞加莱截面部分需要确认取点的时间间隔是否正确，是否在每个周期内取一个点，并且排除了瞬态阶段的数据。可能需要调整索引范围，确保只取稳态后的点。 对于相图重构，用户之前可能使用时间延迟法，需要从解算后的sol中提取位移x1，并生成嵌入向量。这部分代码需要添加到现有代码中，可能需要单独的子图来展示。 在整合过程中，还需要考虑代码的结构，将参数设置、动力学方程、求解过程、可视化分步骤组织，确保代码清晰易读。同时，检查绘图部分的标签、标题、格式是否正确，比如坐标轴的单位转换（如1e-9米转换为纳米）是否恰当。 最后，要验证整个代码的运行是否顺畅，生成的分岔图、庞加莱截面和相图是否能够正确反映AFM系统的动力学特性，如混沌、周期振荡等行为。可能需要调整参数或求解器的设置（如ODE求解的步长、容差）以确保数值稳定性。 总结来说，我需要系统性地修正现有代码的错误，整合相图重构部分，确保所有可视化元素正确呈现，并验证整个流程的准确性和可靠性，最终提供一个完整、可执行的Python脚本，满足用户对AFM系统动力学分析的需求。</think>以下是整合了分岔图、庞加莱截面和相空间重构的完整Python代码（已修正语法错误并优化格式）： ```python import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint # ================== AFM系统动力学方程 ================== def afm_system(state, t, params): x1, x2 = state m, k, D, sigma, Z0, mu, f_amp, omega = params F_spring = -k * x1 / m F_vdW = D / (m * (Z0 + x1)**2) F_repulsive = (sigma**6 * D) / (30 * m * (Z0 + x1)**8) F_excitation = (f_amp / m) * np.cos(omega * t) dx1 = x2 dx2 = F_spring - F_vdW + F_repulsive - (mu/m)*x2 + F_excitation return [dx1, dx2] # ================== 参数设置 ================== params = [ 1e-12, # m [kg] 0.1, # k [N/m] 1e-19, # D [J·m] 3e-9, # sigma [m] 4.5e-9, # Z0 [m] 5e-8, # mu [kg/s] 2e-9, # f_amp [N] 2 * np.pi * 1e4 # omega [rad/s] ] # ================== 分岔图计算 ================== f_amp_range = np.linspace(1e-10, 3e-9, 200) # 力幅扫描范围 max_x1 = [] for f_amp in f_amp_range: params[6] = f_amp # 更新激励幅值 # 瞬态过程（前200周期） t_transient = np.linspace(0, 200*(2*np.pi/params[7]), 2000) sol = odeint(afm_system, [1e-9, 0], t_transient, args=(params,)) # 稳态采样（后50周期） t_steady = np.linspace(0, 50*(2*np.pi/params[7]), 500) sol = odeint(afm_system, sol[-1], t_steady, args=(params,)) max_x1.append(np.max(sol[:, 0])) # ================== 庞加莱截面计算 ================== params[6] = 2.5e-9 # 固定激励幅值 total_time = 1000*(2*np.pi/params[7]) # 1000个周期 t = np.linspace(0, total_time, 200000) sol = odeint(afm_system, [1e-9, 0], t, args=(params,)) # 每个周期取1个点（排除前80000个瞬态点） poincare_idx = np.arange(80000, 100000, int(2*np.pi/params[7]/(t[1]-t[0]))) poincare_section = sol[poincare_idx] # ================== 相空间重构（时间延迟法） ================== def phase_space_reconstruction(signal, tau, m): N = len(signal) - (m-1)*tau return np.array([signal[i:i+(m-1)*tau+1:tau] for i in range(N)]) tau = 50 # 延迟时间（根据数据特性调整） embed_dim = 3 # 嵌入维度 phase_space = phase_space_reconstruction(sol[:, 0], tau, embed_dim) # ================== 科学可视化 ================== plt.figure(figsize=(18, 6), dpi=100) # 分岔图 plt.subplot(1, 3, 1) plt.plot(f_amp_range/1e-9, np.array(max_x1)/1e-9, 'b.', markersize=3, alpha=0.6) plt.xlabel(r'$f_{amp}\ (10^{-9}N)$', fontsize=12) plt.ylabel(r'$\max(x_1)\ (10^{-9}m)$', fontsize=12) plt.title('Bifurcation Diagram', fontsize=14) plt.grid(True, alpha=0.3) # 庞加莱截面 plt.subplot(1, 3, 2) sc = plt.scatter(poincare_section[:,0]/1e-9, poincare_section[:,1]/1e-3, c=np.arange(len(poincare_section)), cmap='viridis', s=10, alpha=0.7, edgecolors='none') plt.xlabel(r'$x_1\ (10^{-9}m)$', fontsize=12) plt.ylabel(r'$x_2\ (10^{-3}m/s)$', fontsize=12) plt.title('Poincaré Section', fontsize=14) plt.colorbar(sc, label='Iteration Index') plt.grid(True, alpha=0.3) # 相空间重构 plt.subplot(1, 3, 3) plt.plot(phase_space[:-tau, 0]/1e-9, phase_space[:-tau, 1]/1e-9, 'g-', lw=0.5, alpha=0.7) plt.xlabel(r'$x(t)\ (10^{-9}m)$', fontsize=12) plt.ylabel(r'$x(t+\tau)\ (10^{-9}m)$', fontsize=12) plt.title('Phase Space Reconstruction', fontsize=14) plt.grid(True, alpha=0.3) plt.tight_layout() plt.savefig('afm_full_analysis.png', bbox_inches='tight', dpi=300) plt.show() ``` 代码说明： 1. **动力学系统**：使用Van der Waals力、弹簧力和周期激励组成的微分方程描述AFM探针运动 2. **分岔图**： - 扫描激励幅值参数$f_{amp}$ - 每个参数值下排除瞬态后记录最大位移 - 展现系统随参数变化的稳定性特征 3. **庞加莱截面**： - 固定参数时，在相空间周期截取状态点 - 点的分布模式反映系统周期/混沌特性 4. **相空间重构**： - 使用时延嵌入法从单变量时间序列重建相空间 - 参数tau需要根据信号自相关函数确定 5. **可视化增强**： - 统一使用物理量纲归一化显示 - 添加颜色映射显示演化过程 - 高分辨率输出（300dpi） 注意： 1. 计算时间较长（约5-10分钟） 2. 可调整tau和embed_dim优化相空间重构效果 3. 分岔图分辨率可通过增加采样点提升 4. 庞加莱截面的稳定性需要足够长的瞬态排除时间